libMesh
Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
libMesh::DynaIO Class Reference

Reading and writing meshes in (a subset of) LS-DYNA format. More...

#include <dyna_io.h>

Inheritance diagram for libMesh::DynaIO:
[legend]

Classes

struct  ElementDefinition
 Defines mapping from libMesh element types to LS-DYNA element types or vice-versa. More...
 
struct  ElementMaps
 struct which holds a map from LS-DYNA to libMesh element numberings and vice-versa. More...
 

Public Types

typedef int32_t dyna_int_type
 The integer type DYNA uses. More...
 

Public Member Functions

 DynaIO (MeshBase &mesh, bool keep_spline_nodes=true)
 Constructor. More...
 
virtual void read (const std::string &name) override
 Reads in a mesh in the Dyna format from the ASCII file given by name. More...
 
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-defined constraint rows to sys. More...
 
void clear_spline_nodes ()
 Removes any spline nodes (both NodeElem and Node), leaving only the FE mesh generated from those splines. More...
 
bool is_parallel_format () const
 Returns true iff this mesh file format and input class are parallelized, so that all processors can read their share of the data at once. More...
 

Static Public Member Functions

static const ElementDefinitionfind_elem_definition (dyna_int_type dyna_elem, int dim, int p)
 Finds the ElementDefinition corresponding to a particular element type. More...
 
static const ElementDefinitionfind_elem_definition (ElemType libmesh_elem, int dim, int p)
 

Protected Member Functions

MeshBasemesh ()
 
void set_n_partitions (unsigned int n_parts)
 Sets the number of partitions in the mesh. More...
 
void skip_comment_lines (std::istream &in, const char comment_start)
 Reads input from in, skipping all the lines that start with the character comment_start. More...
 

Protected Attributes

std::vector< bool > elems_of_dimension
 A vector of bools describing what dimension elements have been encountered when reading a mesh. More...
 

Private Types

typedef double dyna_fp_type
 The floating-point type DYNA uses. More...
 

Private Member Functions

void read_mesh (std::istream &in)
 Implementation of the read() function. More...
 

Static Private Member Functions

static ElementMaps build_element_maps ()
 A static function used to construct the _element_maps struct, statically. More...
 

Private Attributes

std::vector< Node * > spline_node_ptrs
 
std::unordered_map< Node *, Elem * > spline_nodeelem_ptrs
 
bool _keep_spline_nodes
 Whether to keep or eventually discard spline nodes. More...
 

Static Private Attributes

static const int max_ints_per_line = 10
 How many can we find on a line? More...
 
static const int max_fps_per_line = 5
 How many can we find on a line? More...
 
static ElementMaps _element_maps = DynaIO::build_element_maps()
 A static ElementMaps object that is built statically and used by all instances of this class. More...
 

Detailed Description

Reading and writing meshes in (a subset of) LS-DYNA format.

The initial implementation only handles cards in the format described in "Geometry import to LS-DYNA", for isogeometric analysis.

Author
Roy H. Stogner
Date
2019

Definition at line 52 of file dyna_io.h.

Member Typedef Documentation

◆ dyna_fp_type

typedef double libMesh::DynaIO::dyna_fp_type
private

The floating-point type DYNA uses.

Definition at line 178 of file dyna_io.h.

◆ dyna_int_type

The integer type DYNA uses.

Definition at line 112 of file dyna_io.h.

Constructor & Destructor Documentation

◆ DynaIO()

libMesh::DynaIO::DynaIO ( MeshBase mesh,
bool  keep_spline_nodes = true 
)
explicit

Constructor.

Takes a non-const Mesh reference which it will fill up with elements via the read() command.

By default keep_spline_nodes is true, and when reading a BEXT file we retain the spline nodes and their constraints on FE nodes. If keep_spline_nodes is false, we only keep the FE elements, which are no longer constrained to have higher continuity.

Definition at line 130 of file dyna_io.C.

131  :
132  MeshInput<MeshBase> (mesh),
133  _keep_spline_nodes (keep_spline_nodes)
134 {
135 }
bool _keep_spline_nodes
Whether to keep or eventually discard spline nodes.
Definition: dyna_io.h:161

Member Function Documentation

◆ add_spline_constraints()

void libMesh::DynaIO::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-defined constraint rows to sys.

Definition at line 737 of file dyna_io.C.

740 {
741 }

◆ build_element_maps()

DynaIO::ElementMaps libMesh::DynaIO::build_element_maps ( )
staticprivate

A static function used to construct the _element_maps struct, statically.

Definition at line 55 of file dyna_io.C.

References libMesh::DynaIO::ElementMaps::add_def(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX27, libMesh::HEX8, libMesh::QUAD4, and libMesh::QUAD9.

56 {
57  // Object to be filled up
58  ElementMaps em;
59 
60  // em.add_def(ElementDefinition(TRI3, 2, 2, 1)); // node mapping?
61  // em.add_def(ElementDefinition(TRI6, 2, 2, 2)); // node mapping?
62  // em.add_def(ElementDefinition(TET4, 2, 3, 1)); // node mapping?
63  // em.add_def(ElementDefinition(TET10, 2, 3, 2)); // node mapping?
64  // em.add_def(ElementDefinition(PRISM6, 3, 3, 1)); // node mapping?
65  // em.add_def(ElementDefinition(PRISM18, 3, 3, 2)); // node mapping?
66 
67  // Eventually we'll map both tensor-product and non-tensor-product
68  // BEXT elements to ours; for now we only support non-tensor-product
69  // Bezier elements.
70  // for (unsigned int i=0; i != 2; ++i)
71  for (unsigned int i=1; i != 2; ++i)
72  {
73  // We only have one element for whom node orders match...
74  em.add_def(ElementDefinition(EDGE2, i, 1, 1));
75 
76  em.add_def(ElementDefinition(EDGE3, i, 1, 2,
77  {0, 2, 1}));
78  em.add_def(ElementDefinition(EDGE4, i, 1, 2,
79  {0, 2, 3, 1}));
80 
81  em.add_def(ElementDefinition(QUAD4, i, 2, 1,
82  {0, 1, 3, 2}));
83  em.add_def(ElementDefinition(QUAD9, i, 2, 2,
84  {0, 4, 1, 7, 8, 5, 3, 6, 2}));
85 
86  em.add_def(ElementDefinition(HEX8, i, 3, 1,
87  {0, 1, 3, 2, 4, 5, 7, 6}));
88  em.add_def(ElementDefinition(HEX27, i, 3, 2,
89  { 0, 8, 1, 11, 20, 9, 3, 10, 2,
90  12, 21, 13, 24, 26, 22, 15, 23, 14,
91  4, 16, 5, 19, 25, 17, 7, 18, 6}));
92  }
93 
94  return em;
95 }

◆ clear_spline_nodes()

void libMesh::DynaIO::clear_spline_nodes ( )

Removes any spline nodes (both NodeElem and Node), leaving only the FE mesh generated from those splines.

Also removes node constraints to the now-missing nodes.

Definition at line 744 of file dyna_io.C.

References libMesh::MeshTools::clear_spline_nodes().

745 {
746  libmesh_deprecated();
747 
749 }
void clear_spline_nodes(MeshBase &)
Remove spline node (for IsoGeometric Analysis meshes) elements and nodes and constraints from the mes...
Definition: mesh_tools.C:999

◆ find_elem_definition() [1/2]

const DynaIO::ElementDefinition & libMesh::DynaIO::find_elem_definition ( dyna_int_type  dyna_elem,
int  dim,
int  p 
)
static

Finds the ElementDefinition corresponding to a particular element type.

Definition at line 754 of file dyna_io.C.

References _element_maps, dim, and libMesh::DynaIO::ElementMaps::in.

Referenced by libMesh::ExodusII_IO::read(), and read_mesh().

756 {
757  auto eletypes_it =
758  _element_maps.in.find(std::make_tuple(dyna_elem, dim, p));
759 
760  // Make sure we actually found something
761  libmesh_error_msg_if
762  (eletypes_it == _element_maps.in.end(),
763  "Element of type " << dyna_elem <<
764  " dim " << dim <<
765  " degree " << p << " not found!");
766 
767  return eletypes_it->second;
768 }
unsigned int dim
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
Definition: dyna_io.h:209
std::map< Key, ElementDefinition > in
Definition: dyna_io.h:202

◆ find_elem_definition() [2/2]

static const ElementDefinition& libMesh::DynaIO::find_elem_definition ( ElemType  libmesh_elem,
int  dim,
int  p 
)
static

◆ is_parallel_format()

bool libMesh::MeshInput< MeshBase >::is_parallel_format ( ) const
inlineinherited

Returns true iff this mesh file format and input class are parallelized, so that all processors can read their share of the data at once.

Definition at line 87 of file mesh_input.h.

References libMesh::MeshInput< MT >::_is_parallel_format.

87 { return this->_is_parallel_format; }
const bool _is_parallel_format
Flag specifying whether this format is parallel-capable.
Definition: mesh_input.h:130

◆ mesh()

MeshBase & libMesh::MeshInput< MeshBase >::mesh ( )
inlineprotectedinherited
Returns
The object as a writable reference.

Definition at line 178 of file mesh_input.h.

Referenced by libMesh::GMVIO::_read_one_cell(), libMesh::VTKIO::cells_to_vtk(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::Nemesis_IO::copy_elemental_solution(), libMesh::ExodusII_IO::copy_nodal_solution(), libMesh::TetGenIO::element_in(), libMesh::UNVIO::elements_in(), libMesh::UNVIO::elements_out(), libMesh::VTKIO::get_local_node_values(), libMesh::ExodusII_IO::get_sideset_data_indices(), libMesh::UNVIO::groups_in(), libMesh::TetGenIO::node_in(), libMesh::UNVIO::nodes_in(), libMesh::UNVIO::nodes_out(), libMesh::VTKIO::nodes_to_vtk(), libMesh::Nemesis_IO::prepare_to_write_nodal_data(), libMesh::GMVIO::read(), libMesh::Nemesis_IO::read(), libMesh::ExodusII_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read(), libMesh::VTKIO::read(), libMesh::CheckpointIO::read_bcs(), libMesh::CheckpointIO::read_connectivity(), libMesh::ExodusII_IO::read_header(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::UCDIO::read_implementation(), libMesh::UNVIO::read_implementation(), libMesh::GmshIO::read_mesh(), read_mesh(), libMesh::CheckpointIO::read_nodes(), libMesh::CheckpointIO::read_nodesets(), libMesh::CheckpointIO::read_remote_elem(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::ExodusII_IO::read_sideset_data(), libMesh::OFFIO::read_stream(), libMesh::MatlabIO::read_stream(), libMesh::CheckpointIO::read_subdomain_names(), libMesh::TetGenIO::write(), libMesh::Nemesis_IO::write(), libMesh::XdrIO::write(), libMesh::CheckpointIO::write(), libMesh::ExodusII_IO::write(), libMesh::GMVIO::write_ascii_new_impl(), libMesh::GMVIO::write_ascii_old_impl(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::Nemesis_IO::write_element_data(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO::write_elemsets(), libMesh::UCDIO::write_header(), libMesh::UCDIO::write_implementation(), libMesh::UCDIO::write_interior_elems(), libMesh::GmshIO::write_mesh(), libMesh::UCDIO::write_nodal_data(), libMesh::VTKIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_common(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::UCDIO::write_nodes(), libMesh::CheckpointIO::write_nodesets(), libMesh::XdrIO::write_parallel(), libMesh::GmshIO::write_post(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::ExodusII_IO::write_sideset_data(), libMesh::UCDIO::write_soln(), and libMesh::CheckpointIO::write_subdomain_names().

179 {
180  libmesh_error_msg_if(_obj == nullptr, "ERROR: _obj should not be nullptr!");
181  return *_obj;
182 }
MeshBase * _obj
A pointer to a non-const object object.
Definition: mesh_input.h:123

◆ read()

void libMesh::DynaIO::read ( const std::string &  name)
overridevirtual

Reads in a mesh in the Dyna format from the ASCII file given by name.

Note
The user is responsible for calling Mesh::prepare_for_use() after reading the mesh and before using it.
To safely use DynaIO::add_spline_constraints with a DistributedMesh, currently the user must allow_remote_element_removal(false) and allow_renumbering(false) before the mesh is read.

The patch ids defined in the Dyna file are stored as subdomain ids.

The spline nodes defined in the Dyna file are added to the mesh with type NodeElem. The only connection between spline nodes and finite element nodes will be user constraint equations, so using a space-filling-curve partitioner for these meshes might be a good idea.

Implements libMesh::MeshInput< MeshBase >.

Definition at line 139 of file dyna_io.C.

References libMesh::libmesh_assert(), libMesh::Quality::name(), read_mesh(), and libMesh::Utility::unzip_file().

Referenced by main(), libMesh::NameBasedIO::read(), MeshInputTest::testDynaFileMappings(), MeshInputTest::testDynaNoSplines(), MeshInputTest::testDynaReadElem(), and MeshInputTest::testDynaReadPatch().

140 {
141  const bool gzipped_file = (name.rfind(".gz") == name.size() - 3);
142  // These will be handled in unzip_file:
143  // const bool bzipped_file = (name.size() - name.rfind(".bz2") == 4);
144  // const bool xzipped_file = (name.size() - name.rfind(".xz") == 3);
145 
146  std::unique_ptr<std::istream> in;
147 
148  if (gzipped_file)
149  {
150 #ifdef LIBMESH_HAVE_GZSTREAM
151  igzstream * inf = new igzstream;
152  libmesh_assert(inf);
153  in.reset(inf);
154  inf->open(name.c_str(), std::ios::in);
155 #else
156  libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
157 #endif
158  }
159  else
160  {
161  std::ifstream * inf = new std::ifstream;
162  libmesh_assert(inf);
163  in.reset(inf);
164 
165  std::string new_name = Utility::unzip_file(name);
166 
167  inf->open(new_name.c_str(), std::ios::in);
168  }
169 
170  libmesh_assert(in.get());
171 
172  this->read_mesh (*in);
173 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
void read_mesh(std::istream &in)
Implementation of the read() function.
Definition: dyna_io.C:177
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 ...
Definition: utility.C:156
libmesh_assert(ctx)

◆ read_mesh()

void libMesh::DynaIO::read_mesh ( std::istream &  in)
private

Implementation of the read() function.

This function is called by the public interface function and implements reading the file.

Definition at line 177 of file dyna_io.C.

References _keep_spline_nodes, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_node_datum(), libMesh::MeshBase::add_point(), libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::build(), libMesh::MeshBase::clear(), libMesh::MeshTools::clear_spline_nodes(), find_elem_definition(), libMesh::MeshBase::get_constraint_rows(), libMesh::DofObject::invalid_id, libMesh::make_range(), max_fps_per_line, max_ints_per_line, libMesh::MeshInput< MeshBase >::mesh(), libMesh::MeshInput< MT >::mesh(), libMesh::MeshTools::n_elem(), libMesh::NODEELEM, libMesh::DynaIO::ElementDefinition::nodes, libMesh::ParallelObject::processor_id(), libMesh::RATIONAL_BERNSTEIN_MAP, libMesh::Real, libMesh::MeshBase::set_default_mapping_data(), libMesh::MeshBase::set_default_mapping_type(), libMesh::DofObject::set_extra_datum(), libMesh::Elem::set_node(), spline_node_ptrs, spline_nodeelem_ptrs, libMesh::Elem::subdomain_id(), and libMesh::DynaIO::ElementDefinition::type.

Referenced by read().

178 {
179  MeshBase & mesh = MeshInput<MeshBase>::mesh();
180 
181  // This is a serial-only process for now;
182  // the Mesh should be read on processor 0 and
183  // broadcast later
184  libmesh_assert_equal_to (mesh.processor_id(), 0);
185 
186  libmesh_error_msg_if(!in.good(), "Can't read input stream");
187 
188  // clear any data in the mesh
189  mesh.clear();
190 
191  // clear any of our own data
192  spline_node_ptrs.clear();
193  spline_nodeelem_ptrs.clear();
194 
195  // Expect different sections, in this order, perhaps with blank
196  // lines and/or comments in between:
197 
198  enum FileSection {
199  FILE_HEADER,
200  PATCH_HEADER,
201  NODE_LINES,
202  // (repeat for each node)
203  N_ELEM_SUBBLOCKS,
204  ELEM_SUBBLOCK_HEADER,
205  // (repeat for each subblock)
206  ELEM_NODES_LINES,
207  ELEM_COEF_VEC_IDS,
208  // (repeat nodes lines + coef vec ids for each elem, subblock)
209  N_COEF_BLOCKS, // number of coef vec blocks of each type
210  N_VECS_PER_BLOCK, // number of coef vecs in each dense block
211  COEF_VEC_COMPONENTS,
212  // (repeat coef vec components as necessary)
213  // (repeat coef blocks as necessary)
214  //
215  // reserved for sparse block stuff we don't support yet
216  END_OF_FILE };
217 
218  FileSection section = FILE_HEADER;
219 
220  // Values to remember from section to section
221  dyna_int_type patch_id, n_spline_nodes, n_elem, n_coef_vec, weight_control_flag;
222  dyna_int_type n_elem_blocks, n_dense_coef_vec_blocks;
223  std::vector<dyna_int_type> // indexed from 0 to n_elem_blocks
224  block_elem_type,
225  block_n_elem,
226  block_n_nodes, // Number of *spline* nodes constraining elements
227  block_n_coef_vec, // Number of coefficient vectors for each elem
228  block_p,
229  block_dim;
230  std::vector<dyna_int_type> // indexed from 0 to n_dense_coef_vec_blocks
231  n_coef_vecs_in_subblock, n_coef_comp;
232  unsigned char weight_index = 0;
233  dyna_int_type n_nodes_read = 0,
234  n_elem_blocks_read = 0,
235  n_elems_read = 0,
236  n_elem_nodes_read = 0,
237  n_elem_cvids_read = 0,
238  n_coef_headers_read = 0,
239  n_coef_blocks_read = 0,
240  n_coef_comp_read = 0,
241  n_coef_vecs_read = 0;
242 
243  // For reading the file line by line
244  std::string s;
245 
246  // For storing global (spline) weights, until we have enough data to
247  // use them for calculating local (Bezier element) nodes
248  std::vector<Real> spline_weights;
249 
250  // For storing elements' constraint equations:
251  // Global node indices (1-based in file, 0-based in memory):
252  // elem_global_nodes[block_num][elem_num][local_index] = global_index
253  std::vector<std::vector<std::vector<dof_id_type>>> elem_global_nodes;
254 
255  // Dense constraint vectors in the Dyna file
256  // When first read:
257  // dense_constraint_vecs[block_num][vec_in_block][column_num] = coef
258  // When used:
259  // dense_constraint_vecs[0][vec_num][column_num] = coef
260  std::vector<std::vector<std::vector<Real>>> dense_constraint_vecs;
261 
262  // Constraint vector indices (1-based in file, 0-based in memory):
263  // elem_constraint_rows[block_num][elem_num][row_num] = cv_index
264  std::vector<std::vector<std::vector<dof_id_type>>> elem_constraint_rows;
265 
266  while (true)
267  {
268  // Try to read something. This may set EOF!
269  std::getline(in, s);
270 
271  if (in)
272  {
273  // Process s...
274 
275  if (s.find("B E X T 2.0") == static_cast<std::string::size_type>(0))
276  {
277  libmesh_error_msg_if(section != FILE_HEADER,
278  "Found 'B E X T 2.0' outside file header?");
279 
280  section = PATCH_HEADER;
281  continue;
282  }
283 
284  // Ignore comments
285  if (s.find("$#") == static_cast<std::string::size_type>(0))
286  continue;
287 
288  // Ignore blank lines
289  if (s.find_first_not_of(" \t") == std::string::npos)
290  continue;
291 
292  // Parse each important section
293  std::stringstream stream(s);
294  switch (section) {
295  case PATCH_HEADER:
296  stream >> patch_id;
297  stream >> n_spline_nodes;
298  stream >> n_elem;
299  stream >> n_coef_vec;
300  stream >> weight_control_flag;
301  libmesh_error_msg_if(stream.fail(), "Failure to parse patch header\n");
302 
303  spline_node_ptrs.resize(n_spline_nodes);
304  spline_weights.resize(n_spline_nodes);
305 
306  // Even if we have w=1.0, we still want RATIONAL_BERNSTEIN
307  // elements!
308  // if (weight_control_flag)
309  {
310  // If we ever add more nodes that aren't in this file,
311  // merge this mesh with a non-spline mesh, etc., 1.0
312  // is a good default weight
313  const Real default_weight = 1.0;
314  weight_index = cast_int<unsigned char>
315  (mesh.add_node_datum<Real>("rational_weight", true,
316  &default_weight));
318  mesh.set_default_mapping_data(weight_index);
319  }
320 
321  section = NODE_LINES;
322  break;
323  case NODE_LINES:
324  {
325  std::array<dyna_fp_type, 4> xyzw;
326  stream >> xyzw[0];
327  stream >> xyzw[1];
328  stream >> xyzw[2];
329  stream >> xyzw[3];
330 
331  if (weight_control_flag)
332  spline_weights[n_nodes_read] = xyzw[3];
333 
334  // We'll use the spline nodes via NodeElem as the "true"
335  // degrees of freedom, to which other Elem degrees of
336  // freedom will be tied via constraint equations.
337  Node *n = spline_node_ptrs[n_nodes_read] =
338  mesh.add_point(Point(xyzw[0], xyzw[1], xyzw[2]));
339  Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
340  elem->set_node(0) = n;
341  elem->subdomain_id() = 1; // Separate id to ease Exodus output
342  spline_nodeelem_ptrs[n] = elem;
343  }
344  ++n_nodes_read;
345 
346  libmesh_error_msg_if(stream.fail(), "Failure to parse node line\n");
347 
348  if (n_nodes_read >= n_spline_nodes)
349  section = N_ELEM_SUBBLOCKS;
350  break;
351  case N_ELEM_SUBBLOCKS:
352  stream >> n_elem_blocks;
353  libmesh_error_msg_if(stream.fail(), "Failure to parse n_elem_blocks\n");
354 
355  block_elem_type.resize(n_elem_blocks);
356  block_n_elem.resize(n_elem_blocks);
357  block_n_nodes.resize(n_elem_blocks);
358  block_n_coef_vec.resize(n_elem_blocks);
359  block_p.resize(n_elem_blocks);
360  block_dim.resize(n_elem_blocks);
361 
362  elem_global_nodes.resize(n_elem_blocks);
363  elem_constraint_rows.resize(n_elem_blocks);
364 
365  n_elem_blocks_read = 0;
366  section = ELEM_SUBBLOCK_HEADER;
367  break;
368  case ELEM_SUBBLOCK_HEADER:
369  stream >> block_elem_type[n_elem_blocks_read];
370  stream >> block_n_elem[n_elem_blocks_read];
371  stream >> block_n_nodes[n_elem_blocks_read];
372  stream >> block_n_coef_vec[n_elem_blocks_read];
373  stream >> block_p[n_elem_blocks_read];
374 
375  libmesh_error_msg_if(stream.fail(), "Failure to parse elem block\n");
376 
377  block_dim[n_elem_blocks_read] = 1; // All blocks here are at least 1D
378 
379  dyna_int_type block_other_p; // Check for isotropic p
380  stream >> block_other_p;
381  if (!stream.fail())
382  {
383  block_dim[n_elem_blocks_read] = 2; // Found a second dimension!
384 
385  if (block_other_p != block_p[n_elem_blocks_read])
386  libmesh_not_implemented(); // We don't support p anisotropy
387 
388  stream >> block_other_p;
389  if (!stream.fail())
390  {
391  block_dim[n_elem_blocks_read] = 3; // Found a third dimension!
392 
393  if (block_other_p != block_p[n_elem_blocks_read])
394  libmesh_not_implemented();
395  }
396  }
397 
398  {
399  auto & block_global_nodes = elem_global_nodes[n_elem_blocks_read];
400  auto & block_constraint_rows = elem_constraint_rows[n_elem_blocks_read];
401 
402  block_global_nodes.resize(block_n_elem[n_elem_blocks_read]);
403  block_constraint_rows.resize(block_n_elem[n_elem_blocks_read]);
404 
405  for (auto e : make_range(block_n_elem[n_elem_blocks_read]))
406  {
407  block_global_nodes[e].resize(block_n_nodes[n_elem_blocks_read]);
408  block_constraint_rows[e].resize(block_n_coef_vec[n_elem_blocks_read]);
409  }
410  }
411 
412  n_elem_blocks_read++;
413  if (n_elem_blocks_read == n_elem_blocks)
414  {
415  n_elem_blocks_read = 0;
416  n_elems_read = 0;
417  section = ELEM_NODES_LINES;
418  }
419  break;
420  case ELEM_NODES_LINES:
421  {
422  const int end_node_to_read =
423  std::min(block_n_nodes[n_elem_blocks_read], n_elem_nodes_read + max_ints_per_line);
424  for (; n_elem_nodes_read != end_node_to_read; ++n_elem_nodes_read)
425  {
426  dyna_int_type node_id;
427  stream >> node_id;
428  node_id--;
429  elem_global_nodes[n_elem_blocks_read][n_elems_read][n_elem_nodes_read] = node_id;
430 
431  // Let's assume that our *only* mid-line breaks are
432  // due to the max_ints_per_line limit. This should be
433  // less flexible but better for error detection.
434  libmesh_error_msg_if(stream.fail(), "Failure to parse elem nodes\n");
435  }
436 
437  if (n_elem_nodes_read == block_n_nodes[n_elem_blocks_read])
438  {
439  n_elem_nodes_read = 0;
440  section = ELEM_COEF_VEC_IDS;
441  }
442  }
443  break;
444  case ELEM_COEF_VEC_IDS:
445  {
446  const int end_cvid_to_read =
447  std::min(block_n_coef_vec[n_elem_blocks_read], n_elem_cvids_read + max_ints_per_line);
448  for (; n_elem_cvids_read != end_cvid_to_read; ++n_elem_cvids_read)
449  {
450  dyna_int_type node_cvid;
451  stream >> node_cvid;
452  node_cvid--;
453 
454  elem_constraint_rows[n_elem_blocks_read][n_elems_read][n_elem_cvids_read] = node_cvid;
455 
456  // Let's assume that our *only* mid-line breaks are
457  // due to the max_ints_per_line limit. This should be
458  // less flexible but better for error detection.
459  libmesh_error_msg_if(stream.fail(), "Failure to parse elem cvids\n");
460  }
461  if (n_elem_cvids_read == block_n_nodes[n_elem_blocks_read])
462  {
463  n_elem_cvids_read = 0;
464  n_elems_read++;
465  section = ELEM_NODES_LINES; // Read another elem, nodes first
466  if (n_elems_read == block_n_elem[n_elem_blocks_read])
467  {
468  n_elems_read = 0;
469  n_elem_blocks_read++;
470  if (n_elem_blocks_read == n_elem_blocks)
471  section = N_COEF_BLOCKS; // Move on to coefficient vectors
472  }
473  }
474  }
475  break;
476  case N_COEF_BLOCKS:
477  {
478  stream >> n_dense_coef_vec_blocks;
479  dyna_int_type n_sparse_coef_vec_blocks;
480  stream >> n_sparse_coef_vec_blocks;
481 
482  libmesh_error_msg_if(stream.fail(), "Failure to parse n_*_coef_vec_blocks\n");
483 
484  if (n_sparse_coef_vec_blocks != 0)
485  libmesh_not_implemented();
486 
487  dense_constraint_vecs.resize(n_dense_coef_vec_blocks);
488  n_coef_vecs_in_subblock.resize(n_dense_coef_vec_blocks);
489  n_coef_comp.resize(n_dense_coef_vec_blocks);
490 
491  section = N_VECS_PER_BLOCK;
492  }
493  break;
494  case N_VECS_PER_BLOCK:
495  stream >> n_coef_vecs_in_subblock[n_coef_headers_read];
496  stream >> n_coef_comp[n_coef_headers_read];
497 
498  libmesh_error_msg_if(stream.fail(), "Failure to parse dense coef subblock header\n");
499 
500  dense_constraint_vecs[n_coef_headers_read].resize
501  (n_coef_vecs_in_subblock[n_coef_headers_read]);
502 
503  for (auto & vec : dense_constraint_vecs[n_coef_headers_read])
504  vec.resize(n_coef_comp[n_coef_headers_read]);
505 
506  n_coef_headers_read++;
507  if (n_coef_headers_read == n_dense_coef_vec_blocks)
508  {
509  n_coef_headers_read = 0;
510  section = COEF_VEC_COMPONENTS;
511  }
512  break;
513  case COEF_VEC_COMPONENTS:
514  {
515  auto & current_vec =
516  dense_constraint_vecs[n_coef_blocks_read][n_coef_vecs_read];
517 
518  const int end_coef_to_read =
519  std::min(n_coef_comp[n_coef_blocks_read],
520  n_coef_comp_read + max_fps_per_line);
521  for (; n_coef_comp_read != end_coef_to_read; ++n_coef_comp_read)
522  {
523  dyna_fp_type coef_comp;
524  stream >> coef_comp;
525 
526  current_vec[n_coef_comp_read] = coef_comp;
527 
528  // Let's assume that our *only* mid-line breaks are
529  // due to the max_fps_per_line limit. This should be
530  // less flexible but better for error detection.
531  libmesh_error_msg_if(stream.fail(), "Failure to parse coefficients\n");
532  }
533  if (n_coef_comp_read == n_coef_comp[n_coef_blocks_read])
534  {
535  n_coef_comp_read = 0;
536  n_coef_vecs_read++;
537  if (n_coef_vecs_read == n_coef_vecs_in_subblock[n_coef_blocks_read])
538  {
539  n_coef_vecs_read = 0;
540  n_coef_blocks_read++;
541  if (n_coef_blocks_read == n_dense_coef_vec_blocks)
542  section = END_OF_FILE;
543  }
544  }
545  }
546  break;
547  default:
548  libmesh_error();
549  }
550 
551  if (section == END_OF_FILE)
552  break;
553  } // if (in)
554  else if (in.eof())
555  libmesh_error_msg("Premature end of file");
556  else
557  libmesh_error_msg("Input stream failure! Perhaps the file does not exist?");
558  }
559 
560  // Merge dense_constraint_vecs blocks
561  if (n_dense_coef_vec_blocks)
562  for (auto coef_vec_block :
563  IntRange<dyna_int_type>(1, n_dense_coef_vec_blocks))
564  {
565  auto & dcv0 = dense_constraint_vecs[0];
566  auto & dcvi = dense_constraint_vecs[coef_vec_block];
567  dcv0.insert(dcv0.end(),
568  std::make_move_iterator(dcvi.begin()),
569  std::make_move_iterator(dcvi.end()));
570  }
571  dense_constraint_vecs.resize(1);
572 
573  // Constraint matrices:
574  // elem_constraint_mat[block_num][elem_num][local_node_index][elem_global_nodes_index] = c
575  std::vector<std::vector<std::vector<std::vector<Real>>>> elem_constraint_mat(n_elem_blocks);
576 
577  // We need to calculate local nodes on the fly, and we'll be
578  // calculating them from constraint matrix columns, and we'll need
579  // to make sure that the same node is found each time it's
580  // calculated from multiple neighboring elements.
581  std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
582 
583  auto & constraint_rows = mesh.get_constraint_rows();
584 
585  for (auto block_num : make_range(n_elem_blocks))
586  {
587  elem_constraint_mat[block_num].resize(block_n_elem[block_num]);
588 
589  for (auto elem_num :
590  make_range(block_n_elem[block_num]))
591  {
592  // Consult the import element table to determine which element to build
593  const ElementDefinition & elem_defn =
594  find_elem_definition(block_elem_type[block_num],
595  block_dim[block_num],
596  block_p[block_num]);
597 
598  auto elem = Elem::build(elem_defn.type);
599  libmesh_error_msg_if(elem->dim() != block_dim[block_num],
600  "Elem dim " << elem->dim() <<
601  " != block_dim " << block_dim[block_num]);
602 
603  auto & my_constraint_rows = elem_constraint_rows[block_num][elem_num];
604  auto & my_global_nodes = elem_global_nodes[block_num][elem_num];
605  auto & my_constraint_mat = elem_constraint_mat[block_num][elem_num];
606 
607  my_constraint_mat.resize(block_n_coef_vec[block_num]);
608  for (auto spline_node_index :
609  make_range(block_n_coef_vec[block_num]))
610  my_constraint_mat[spline_node_index].resize(elem->n_nodes());
611 
612  for (auto spline_node_index :
613  make_range(block_n_coef_vec[block_num]))
614  {
615  // Find which coef block this elem's vectors are from
616  const dyna_int_type elem_coef_vec_index =
617  my_constraint_rows[spline_node_index];
618 
619  dyna_int_type coef_block_num = 0;
620  dyna_int_type first_block_coef_vec = 0;
621  for (; elem_coef_vec_index >= first_block_coef_vec &&
622  coef_block_num != n_dense_coef_vec_blocks; ++coef_block_num)
623  {
624  first_block_coef_vec += n_coef_vecs_in_subblock[coef_block_num];
625  }
626 
627  // Make sure we did find a valid coef block
628  libmesh_error_msg_if(coef_block_num == n_dense_coef_vec_blocks &&
629  first_block_coef_vec <= elem_coef_vec_index,
630  "Can't find valid constraint coef vector");
631 
632  coef_block_num--;
633 
634  libmesh_error_msg_if
635  (dyna_int_type(elem->n_nodes()) != n_coef_comp[coef_block_num],
636  "Found " << n_coef_comp[coef_block_num] <<
637  " constraint coef vectors for " <<
638  elem->n_nodes() << " nodes");
639 
640  for (auto elem_node_index :
641  make_range(elem->n_nodes()))
642  my_constraint_mat[spline_node_index][elem_node_index] =
643  dense_constraint_vecs[0][elem_coef_vec_index][elem_node_index];
644  }
645 
646  for (auto elem_node_index :
647  make_range(elem->n_nodes()))
648  {
649  dof_id_type global_node_idx = DofObject::invalid_id;
650 
651  // New finite element node data: dot product of
652  // constraint matrix columns with spline node data.
653  // Store that column as a key.
654  std::vector<std::pair<dof_id_type, Real>> key;
655 
656  for (auto spline_node_index :
657  make_range(block_n_coef_vec[block_num]))
658  {
659  const dyna_int_type elem_coef_vec_index =
660  my_constraint_rows[spline_node_index];
661 
662  const Real coef =
663  libmesh_vector_at(dense_constraint_vecs[0],
664  elem_coef_vec_index)[elem_node_index];
665 
666  // Global nodes are supposed to be in sorted order
667  if (global_node_idx != DofObject::invalid_id)
668  libmesh_error_msg_if(my_global_nodes[spline_node_index] <= global_node_idx,
669  "Found unsorted global node");
670 
671  global_node_idx = my_global_nodes[spline_node_index];
672 
673  if (coef != 0) // Ignore irrelevant spline nodes
674  key.emplace_back(global_node_idx, coef);
675  }
676 
677  auto local_node_it = local_nodes.find(key);
678 
679  if (local_node_it != local_nodes.end())
680  elem->set_node(elem_defn.nodes[elem_node_index]) =
681  local_node_it->second;
682  else
683  {
684  Point p(0);
685  Real w = 0;
686  std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
687 
688  for (auto spline_node_index :
689  make_range(block_n_coef_vec[block_num]))
690  {
691  const dof_id_type my_node_idx =
692  my_global_nodes[spline_node_index];
693 
694  const dyna_int_type elem_coef_vec_index =
695  my_constraint_rows[spline_node_index];
696 
697  Node * spline_node =
698  libmesh_vector_at(spline_node_ptrs,
699  my_node_idx);
700 
701  const Real coef =
702  libmesh_vector_at(dense_constraint_vecs[0],
703  elem_coef_vec_index)[elem_node_index];
704  p.add_scaled(*spline_node, coef);
705  w += coef * libmesh_vector_at(spline_weights,
706  my_node_idx);
707 
708  // We don't need to store 0 entries;
709  // constraint_rows is a sparse structure.
710  if (coef)
711  {
712  Elem * nodeelem =
713  libmesh_map_find(spline_nodeelem_ptrs, spline_node);
714  constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
715  }
716  }
717 
718  Node *n = mesh.add_point(p);
719  if (weight_control_flag)
720  n->set_extra_datum<Real>(weight_index, w);
721  local_nodes[key] = n;
722  elem->set_node(elem_defn.nodes[elem_node_index]) = n;
723 
724  constraint_rows[n] = constraint_row;
725  }
726  }
727 
728  mesh.add_elem(std::move(elem));
729  }
730  }
731 
732  if (!_keep_spline_nodes)
734 }
std::vector< Node * > spline_node_ptrs
Definition: dyna_io.h:155
constraint_rows_type & get_constraint_rows()
Definition: mesh_base.h:1683
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2381
int32_t dyna_int_type
The integer type DYNA uses.
Definition: dyna_io.h:112
double dyna_fp_type
The floating-point type DYNA uses.
Definition: dyna_io.h:178
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:850
void clear_spline_nodes(MeshBase &)
Remove spline node (for IsoGeometric Analysis meshes) elements and nodes and constraints from the mes...
Definition: mesh_tools.C:999
bool _keep_spline_nodes
Whether to keep or eventually discard spline nodes.
Definition: dyna_io.h:161
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 ElementDefinition & find_elem_definition(dyna_int_type dyna_elem, int dim, int p)
Finds the ElementDefinition corresponding to a particular element type.
Definition: dyna_io.C:754
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)
Definition: elem.C:273
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...
Definition: mesh_base.h:800
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:477
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...
Definition: dof_object.h:1119
static const int max_fps_per_line
How many can we find on a line?
Definition: dyna_io.h:183
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:862
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.
Definition: mesh_base.h:2246
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unordered_map< Node *, Elem * > spline_nodeelem_ptrs
Definition: dyna_io.h:156
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...
Definition: int_range.h:134
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...
Definition: mesh_base.h:818
static const int max_ints_per_line
How many can we find on a line?
Definition: dyna_io.h:173
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:67

◆ set_n_partitions()

void libMesh::MeshInput< MeshBase >::set_n_partitions ( unsigned int  n_parts)
inlineprotectedinherited

Sets the number of partitions in the mesh.

Typically this gets done by the partitioner, but some parallel file formats begin "pre-partitioned".

Definition at line 101 of file mesh_input.h.

References libMesh::MeshInput< MT >::mesh().

Referenced by libMesh::Nemesis_IO::read(), and libMesh::XdrIO::read_header().

101 { this->mesh().set_n_partitions() = n_parts; }
unsigned int & set_n_partitions()
Definition: mesh_base.h:1789

◆ skip_comment_lines()

void libMesh::MeshInput< MeshBase >::skip_comment_lines ( std::istream &  in,
const char  comment_start 
)
protectedinherited

Reads input from in, skipping all the lines that start with the character comment_start.

Definition at line 187 of file mesh_input.h.

Referenced by libMesh::TetGenIO::read(), and libMesh::UCDIO::read_implementation().

189 {
190  char c, line[256];
191 
192  while (in.get(c), c==comment_start)
193  in.getline (line, 255);
194 
195  // put back first character of
196  // first non-comment line
197  in.putback (c);
198 }

Member Data Documentation

◆ _element_maps

DynaIO::ElementMaps libMesh::DynaIO::_element_maps = DynaIO::build_element_maps()
staticprivate

A static ElementMaps object that is built statically and used by all instances of this class.

Definition at line 209 of file dyna_io.h.

Referenced by find_elem_definition().

◆ _keep_spline_nodes

bool libMesh::DynaIO::_keep_spline_nodes
private

Whether to keep or eventually discard spline nodes.

Definition at line 161 of file dyna_io.h.

Referenced by read_mesh().

◆ elems_of_dimension

std::vector<bool> libMesh::MeshInput< MeshBase >::elems_of_dimension
protectedinherited

◆ max_fps_per_line

const int libMesh::DynaIO::max_fps_per_line = 5
staticprivate

How many can we find on a line?

Definition at line 183 of file dyna_io.h.

Referenced by read_mesh().

◆ max_ints_per_line

const int libMesh::DynaIO::max_ints_per_line = 10
staticprivate

How many can we find on a line?

Definition at line 173 of file dyna_io.h.

Referenced by read_mesh().

◆ spline_node_ptrs

std::vector<Node *> libMesh::DynaIO::spline_node_ptrs
private

Definition at line 155 of file dyna_io.h.

Referenced by read_mesh().

◆ spline_nodeelem_ptrs

std::unordered_map<Node *, Elem *> libMesh::DynaIO::spline_nodeelem_ptrs
private

Definition at line 156 of file dyna_io.h.

Referenced by read_mesh().


The documentation for this class was generated from the following files: