19 #include "libmesh/libmesh_config.h" 20 #include "libmesh/libmesh_logging.h" 21 #include "libmesh/dof_map.h" 22 #include "libmesh/elem.h" 23 #include "libmesh/dyna_io.h" 24 #include "libmesh/mesh_base.h" 25 #include "libmesh/mesh_tools.h" 26 #include "libmesh/int_range.h" 27 #include "libmesh/utility.h" 33 #ifdef LIBMESH_HAVE_GZSTREAM 34 # include "gzstream.h" 69 for (
unsigned int i=1; i != 2; ++i)
82 {0, 4, 1, 7, 8, 5, 3, 6, 2}));
85 {0, 1, 3, 2, 4, 5, 7, 6}));
87 { 0, 8, 1, 11, 20, 9, 3, 10, 2,
88 12, 21, 13, 24, 26, 22, 15, 23, 14,
89 4, 16, 5, 19, 25, 17, 7, 18, 6}));
103 dyna_type(dyna_type_in),
109 libmesh_not_implemented_msg(
"Support for Polygons/Polyhedra not yet implemented");
111 std::iota(nodes.begin(), nodes.end(), 0);
120 std::vector<unsigned int> && nodes_in) :
122 dyna_type(dyna_type_in),
131 bool keep_spline_nodes) :
133 _keep_spline_nodes (keep_spline_nodes)
143 std::unique_ptr<std::istream> in;
147 #ifdef LIBMESH_HAVE_GZSTREAM 148 auto inf = std::make_unique<igzstream>();
150 inf->open(
name.c_str(), std::ios::in);
153 libmesh_error_msg(
"ERROR: need gzstream to handle .gz files!!!");
158 auto inf = std::make_unique<std::ifstream>();
163 inf->open(new_name.c_str(), std::ios::in);
183 libmesh_error_msg_if(!in.good(),
"Can't read input stream");
201 ELEM_SUBBLOCK_HEADER,
215 FileSection section = FILE_HEADER;
220 std::vector<dyna_int_type>
227 std::vector<dyna_int_type>
228 n_coef_vecs_in_subblock, n_coef_comp;
229 unsigned char weight_index = 0;
231 n_elem_blocks_read = 0,
233 n_elem_nodes_read = 0,
234 n_elem_cvids_read = 0,
235 n_coef_headers_read = 0,
236 n_coef_blocks_read = 0,
237 n_coef_comp_read = 0,
238 n_coef_vecs_read = 0;
245 std::vector<Real> spline_weights;
250 std::vector<std::vector<std::vector<dof_id_type>>> elem_global_nodes;
257 std::vector<std::vector<std::vector<Real>>> dense_constraint_vecs;
261 std::vector<std::vector<std::vector<dof_id_type>>> elem_constraint_rows;
272 if (s.find(
"B E X T 2.0") ==
static_cast<std::string::size_type
>(0))
274 libmesh_error_msg_if(section != FILE_HEADER,
275 "Found 'B E X T 2.0' outside file header?");
277 section = PATCH_HEADER;
282 if (s.find(
"$#") ==
static_cast<std::string::size_type
>(0))
286 if (s.find_first_not_of(
" \t") == std::string::npos)
290 std::stringstream stream(s);
294 stream >> n_spline_nodes;
296 stream >> n_coef_vec;
297 stream >> weight_control_flag;
298 libmesh_error_msg_if(stream.fail(),
"Failure to parse patch header\n");
301 spline_weights.resize(n_spline_nodes);
310 const Real default_weight = 1.0;
311 weight_index = cast_int<unsigned char>
318 section = NODE_LINES;
322 std::array<dyna_fp_type, 4> xyzw;
328 if (weight_control_flag)
329 spline_weights[n_nodes_read] = xyzw[3];
343 libmesh_error_msg_if(stream.fail(),
"Failure to parse node line\n");
345 if (n_nodes_read >= n_spline_nodes)
346 section = N_ELEM_SUBBLOCKS;
348 case N_ELEM_SUBBLOCKS:
349 stream >> n_elem_blocks;
350 libmesh_error_msg_if(stream.fail(),
"Failure to parse n_elem_blocks\n");
352 block_elem_type.resize(n_elem_blocks);
353 block_n_elem.resize(n_elem_blocks);
354 block_n_nodes.resize(n_elem_blocks);
355 block_n_coef_vec.resize(n_elem_blocks);
356 block_p.resize(n_elem_blocks);
357 block_dim.resize(n_elem_blocks);
359 elem_global_nodes.resize(n_elem_blocks);
360 elem_constraint_rows.resize(n_elem_blocks);
362 n_elem_blocks_read = 0;
363 section = ELEM_SUBBLOCK_HEADER;
365 case ELEM_SUBBLOCK_HEADER:
366 stream >> block_elem_type[n_elem_blocks_read];
367 stream >> block_n_elem[n_elem_blocks_read];
368 stream >> block_n_nodes[n_elem_blocks_read];
369 stream >> block_n_coef_vec[n_elem_blocks_read];
370 stream >> block_p[n_elem_blocks_read];
372 libmesh_error_msg_if(stream.fail(),
"Failure to parse elem block\n");
374 block_dim[n_elem_blocks_read] = 1;
377 stream >> block_other_p;
380 block_dim[n_elem_blocks_read] = 2;
382 if (block_other_p != block_p[n_elem_blocks_read])
383 libmesh_not_implemented();
385 stream >> block_other_p;
388 block_dim[n_elem_blocks_read] = 3;
390 if (block_other_p != block_p[n_elem_blocks_read])
391 libmesh_not_implemented();
396 auto & block_global_nodes = elem_global_nodes[n_elem_blocks_read];
397 auto & block_constraint_rows = elem_constraint_rows[n_elem_blocks_read];
399 block_global_nodes.resize(block_n_elem[n_elem_blocks_read]);
400 block_constraint_rows.resize(block_n_elem[n_elem_blocks_read]);
402 for (
auto e :
make_range(block_n_elem[n_elem_blocks_read]))
404 block_global_nodes[e].resize(block_n_nodes[n_elem_blocks_read]);
405 block_constraint_rows[e].resize(block_n_coef_vec[n_elem_blocks_read]);
409 n_elem_blocks_read++;
410 if (n_elem_blocks_read == n_elem_blocks)
412 n_elem_blocks_read = 0;
414 section = ELEM_NODES_LINES;
417 case ELEM_NODES_LINES:
419 const int end_node_to_read =
420 std::min(block_n_nodes[n_elem_blocks_read], n_elem_nodes_read +
max_ints_per_line);
421 for (; n_elem_nodes_read != end_node_to_read; ++n_elem_nodes_read)
426 elem_global_nodes[n_elem_blocks_read][n_elems_read][n_elem_nodes_read] = node_id;
431 libmesh_error_msg_if(stream.fail(),
"Failure to parse elem nodes\n");
434 if (n_elem_nodes_read == block_n_nodes[n_elem_blocks_read])
436 n_elem_nodes_read = 0;
437 section = ELEM_COEF_VEC_IDS;
441 case ELEM_COEF_VEC_IDS:
443 const int end_cvid_to_read =
444 std::min(block_n_coef_vec[n_elem_blocks_read], n_elem_cvids_read +
max_ints_per_line);
445 for (; n_elem_cvids_read != end_cvid_to_read; ++n_elem_cvids_read)
451 elem_constraint_rows[n_elem_blocks_read][n_elems_read][n_elem_cvids_read] = node_cvid;
456 libmesh_error_msg_if(stream.fail(),
"Failure to parse elem cvids\n");
458 if (n_elem_cvids_read == block_n_nodes[n_elem_blocks_read])
460 n_elem_cvids_read = 0;
462 section = ELEM_NODES_LINES;
463 if (n_elems_read == block_n_elem[n_elem_blocks_read])
466 n_elem_blocks_read++;
467 if (n_elem_blocks_read == n_elem_blocks)
468 section = N_COEF_BLOCKS;
475 stream >> n_dense_coef_vec_blocks;
477 stream >> n_sparse_coef_vec_blocks;
479 libmesh_error_msg_if(stream.fail(),
"Failure to parse n_*_coef_vec_blocks\n");
481 if (n_sparse_coef_vec_blocks != 0)
482 libmesh_not_implemented();
484 dense_constraint_vecs.resize(n_dense_coef_vec_blocks);
485 n_coef_vecs_in_subblock.resize(n_dense_coef_vec_blocks);
486 n_coef_comp.resize(n_dense_coef_vec_blocks);
488 section = N_VECS_PER_BLOCK;
491 case N_VECS_PER_BLOCK:
492 stream >> n_coef_vecs_in_subblock[n_coef_headers_read];
493 stream >> n_coef_comp[n_coef_headers_read];
495 libmesh_error_msg_if(stream.fail(),
"Failure to parse dense coef subblock header\n");
497 dense_constraint_vecs[n_coef_headers_read].resize
498 (n_coef_vecs_in_subblock[n_coef_headers_read]);
500 for (
auto & vec : dense_constraint_vecs[n_coef_headers_read])
501 vec.resize(n_coef_comp[n_coef_headers_read]);
503 n_coef_headers_read++;
504 if (n_coef_headers_read == n_dense_coef_vec_blocks)
506 n_coef_headers_read = 0;
507 section = COEF_VEC_COMPONENTS;
510 case COEF_VEC_COMPONENTS:
513 dense_constraint_vecs[n_coef_blocks_read][n_coef_vecs_read];
515 const int end_coef_to_read =
516 std::min(n_coef_comp[n_coef_blocks_read],
518 for (; n_coef_comp_read != end_coef_to_read; ++n_coef_comp_read)
523 current_vec[n_coef_comp_read] = coef_comp;
528 libmesh_error_msg_if(stream.fail(),
"Failure to parse coefficients\n");
530 if (n_coef_comp_read == n_coef_comp[n_coef_blocks_read])
532 n_coef_comp_read = 0;
534 if (n_coef_vecs_read == n_coef_vecs_in_subblock[n_coef_blocks_read])
536 n_coef_vecs_read = 0;
537 n_coef_blocks_read++;
538 if (n_coef_blocks_read == n_dense_coef_vec_blocks)
539 section = END_OF_FILE;
548 if (section == END_OF_FILE)
552 libmesh_error_msg(
"Premature end of file");
554 libmesh_error_msg(
"Input stream failure! Perhaps the file does not exist?");
558 if (n_dense_coef_vec_blocks)
559 for (
auto coef_vec_block :
562 auto & dcv0 = dense_constraint_vecs[0];
563 auto & dcvi = dense_constraint_vecs[coef_vec_block];
564 dcv0.insert(dcv0.end(),
565 std::make_move_iterator(dcvi.begin()),
566 std::make_move_iterator(dcvi.end()));
568 dense_constraint_vecs.resize(1);
572 std::vector<std::vector<std::vector<std::vector<Real>>>> elem_constraint_mat(n_elem_blocks);
578 std::map<std::vector<std::pair<dof_id_type, Real>>,
Node *> local_nodes;
582 for (
auto block_num :
make_range(n_elem_blocks))
584 elem_constraint_mat[block_num].resize(block_n_elem[block_num]);
592 block_dim[block_num],
596 libmesh_error_msg_if(elem->dim() != block_dim[block_num],
597 "Elem dim " << elem->dim() <<
598 " != block_dim " << block_dim[block_num]);
600 auto & my_constraint_rows = elem_constraint_rows[block_num][elem_num];
601 auto & my_global_nodes = elem_global_nodes[block_num][elem_num];
602 auto & my_constraint_mat = elem_constraint_mat[block_num][elem_num];
604 my_constraint_mat.resize(block_n_coef_vec[block_num]);
605 for (
auto spline_node_index :
607 my_constraint_mat[spline_node_index].resize(elem->n_nodes());
609 for (
auto spline_node_index :
614 my_constraint_rows[spline_node_index];
618 for (; elem_coef_vec_index >= first_block_coef_vec &&
619 coef_block_num != n_dense_coef_vec_blocks; ++coef_block_num)
621 first_block_coef_vec += n_coef_vecs_in_subblock[coef_block_num];
625 libmesh_error_msg_if(coef_block_num == n_dense_coef_vec_blocks &&
626 first_block_coef_vec <= elem_coef_vec_index,
627 "Can't find valid constraint coef vector");
632 (
dyna_int_type(elem->n_nodes()) != n_coef_comp[coef_block_num],
633 "Found " << n_coef_comp[coef_block_num] <<
634 " constraint coef vectors for " <<
635 elem->n_nodes() <<
" nodes");
637 for (
auto elem_node_index :
639 my_constraint_mat[spline_node_index][elem_node_index] =
640 dense_constraint_vecs[0][elem_coef_vec_index][elem_node_index];
643 for (
auto elem_node_index :
651 std::vector<std::pair<dof_id_type, Real>> key;
653 for (
auto spline_node_index :
657 my_constraint_rows[spline_node_index];
660 libmesh_vector_at(dense_constraint_vecs[0],
661 elem_coef_vec_index)[elem_node_index];
665 libmesh_error_msg_if(my_global_nodes[spline_node_index] <= global_node_idx,
666 "Found unsorted global node");
668 global_node_idx = my_global_nodes[spline_node_index];
671 key.emplace_back(global_node_idx, coef);
674 if (
const auto local_node_it = local_nodes.find(key);
675 local_node_it != local_nodes.end())
676 elem->set_node(elem_defn.
nodes[elem_node_index],
677 local_node_it->second);
682 std::vector<std::pair<std::pair<const Elem *, unsigned int>,
Real>> constraint_row;
684 for (
auto spline_node_index :
688 my_global_nodes[spline_node_index];
691 my_constraint_rows[spline_node_index];
698 libmesh_vector_at(dense_constraint_vecs[0],
699 elem_coef_vec_index)[elem_node_index];
701 w += coef * libmesh_vector_at(spline_weights,
710 constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
715 if (weight_control_flag)
717 local_nodes[key] = n;
718 elem->set_node(elem_defn.
nodes[elem_node_index], n);
720 constraint_rows[n] = constraint_row;
740 #ifdef LIBMESH_ENABLE_DEPRECATED 743 libmesh_deprecated();
747 #endif // LIBMESH_ENABLE_DEPRECATED 761 "Element of type " << dyna_elem <<
763 " degree " << p <<
" not found!");
765 return eletypes_it->second;
772 int libmesh_dbg_var(
dim),
773 int libmesh_dbg_var(p))
781 "Element of type " << libmesh_elem <<
785 libmesh_assert_equal_to(eletypes_it->second.dim,
dim);
786 libmesh_assert_equal_to(eletypes_it->second.p, p);
789 return eletypes_it->second;
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
ElemType
Defines an enum for geometric element types.
virtual void read(const std::string &name) override
Reads in a mesh in the Dyna format from the ASCII file given by name.
std::vector< Node * > spline_node_ptrs
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
virtual Node *& set_node(const unsigned int i)
A Node is like a Point, but with more information.
int32_t dyna_int_type
The integer type DYNA uses.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
double dyna_fp_type
The floating-point type DYNA uses.
bool ends_with(std::string_view superstring, std::string_view suffix)
Look for a substring at the very end of a string.
void read_mesh(std::istream &in)
Implementation of the read() function.
Defines mapping from libMesh element types to LS-DYNA element types or vice-versa.
This is the base class from which all geometric element types are derived.
void clear_spline_nodes()
Removes any spline nodes (both NodeElem and Node), leaving only the FE mesh generated from those spli...
void add_def(const ElementDefinition &eledef)
bool _keep_spline_nodes
Whether to keep or eventually discard spline nodes.
The libMesh namespace provides an interface to certain functionality in the library.
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
static const unsigned int type_to_n_nodes_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of nodes in the element...
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
This is the MeshBase class.
void add_spline_constraints(DofMap &dof_map, unsigned int sys_num, unsigned int var_num)
Constrains finite element degrees of freedom in terms of spline degrees of freedom by adding user-def...
This class handles the numbering of degrees of freedom on a mesh.
static const ElementDefinition & find_elem_definition(dyna_int_type dyna_elem, int dim, int p)
Finds the ElementDefinition corresponding to a particular element type.
const dof_id_type n_nodes
std::map< ElemType, ElementDefinition > out
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
void set_default_mapping_type(const ElemMappingType type)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
std::string unzip_file(std::string_view name)
Create an unzipped copy of a bz2 or xz file, returning the name of the now-unzipped file that can be ...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
void set_extra_datum(const unsigned int index, const T value)
Sets the value on this object of the extra datum associated with index, which should have been obtain...
static ElementMaps build_element_maps()
A static function used to construct the _element_maps struct, statically.
static const int max_fps_per_line
How many can we find on a line?
virtual void clear()
Deletes all the element and node data that is currently stored.
unsigned int add_node_datum(const std::string &name, bool allocate_data=true, const T *default_value=nullptr)
Register a datum (of type T) to be added to each node in the mesh.
DynaIO(MeshBase &mesh, bool keep_spline_nodes=true)
Constructor.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
std::unordered_map< Node *, Elem * > spline_nodeelem_ptrs
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
void set_default_mapping_data(const unsigned char data)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
static const int max_ints_per_line
How many can we find on a line?
struct which holds a map from LS-DYNA to libMesh element numberings and vice-versa.
std::vector< unsigned int > nodes
processor_id_type processor_id() const
std::map< Key, ElementDefinition > in
ElementDefinition(ElemType type_in, dyna_int_type dyna_type_in, dyna_int_type dim_in, dyna_int_type p_in)
A Point defines a location in LIBMESH_DIM dimensional Real space.