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

#include <dyna_io.h>

Inheritance diagram for libMesh::DynaIO:

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 ElementDefinition &	find_elem_definition (dyna_int_type dyna_elem, int dim, int p)
	Finds the ElementDefinition corresponding to a particular element type. More...

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

Protected Member Functions
MeshBase &	mesh ()

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 182 of file dyna_io.h.

◆ dyna_int_type

typedef int32_t libMesh::DynaIO::dyna_int_type

The integer type DYNA uses.

Definition at line 116 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.

                                         :
   MeshInput<MeshBase> (mesh),
   _keep_spline_nodes  (keep_spline_nodes)
 {
 }

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 733 of file dyna_io.C.

736 {

737 }

◆ 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 53 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.

 {
   // Object to be filled up
   ElementMaps em;
 
   // em.add_def(ElementDefinition(TRI3, 2, 2, 1)); // node mapping?
   // em.add_def(ElementDefinition(TRI6, 2, 2, 2)); // node mapping?
   // em.add_def(ElementDefinition(TET4,  2, 3, 1));   // node mapping?
   // em.add_def(ElementDefinition(TET10,  2, 3, 2));  // node mapping?
   // em.add_def(ElementDefinition(PRISM6,  3, 3, 1)); // node mapping?
   // em.add_def(ElementDefinition(PRISM18, 3, 3, 2)); // node mapping?
 
   // Eventually we'll map both tensor-product and non-tensor-product
   // BEXT elements to ours; for now we only support non-tensor-product
   // Bezier elements.
   // for (unsigned int i=0; i != 2; ++i)
   for (unsigned int i=1; i != 2; ++i)
     {
       // We only have one element for whom node orders match...
       em.add_def(ElementDefinition(EDGE2, i, 1, 1));
 
       em.add_def(ElementDefinition(EDGE3, i, 1, 2,
                                    {0, 2, 1}));
       em.add_def(ElementDefinition(EDGE4, i, 1, 2,
                                    {0, 2, 3, 1}));
 
       em.add_def(ElementDefinition(QUAD4, i, 2, 1,
                                    {0, 1, 3, 2}));
       em.add_def(ElementDefinition(QUAD9, i, 2, 2,
                                    {0, 4, 1, 7, 8, 5, 3, 6, 2}));
 
       em.add_def(ElementDefinition(HEX8,  i, 3, 1,
                                    {0, 1, 3, 2, 4, 5, 7, 6}));
       em.add_def(ElementDefinition(HEX27, i, 3, 2,
                                    { 0,  8,  1, 11, 20,  9,  3, 10,  2,
                                     12, 21, 13, 24, 26, 22, 15, 23, 14,
                                      4, 16,  5, 19, 25, 17,  7, 18,  6}));
     }
 
   return em;
 }

◆ 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.

deprecated - use MeshTools::clear_spline_nodes(mesh) instead.

Definition at line 741 of file dyna_io.C.

References libMesh::MeshTools::clear_spline_nodes().

 {
   libmesh_deprecated();
 
   MeshTools::clear_spline_nodes(MeshInput<MeshBase>::mesh());
 }

◆ 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 752 of file dyna_io.C.

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

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

 {
   auto eletypes_it =
     _element_maps.in.find(std::make_tuple(dyna_elem, dim, p));
 
   // Make sure we actually found something
   libmesh_error_msg_if
     (eletypes_it == _element_maps.in.end(),
      "Element of type " << dyna_elem <<
      " dim " << dim <<
      " degree " << p << " not found!");
 
   return eletypes_it->second;
 }

◆ 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; }

libMesh::MeshInput< MeshBase >::_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.

 {
   libmesh_error_msg_if(_obj == nullptr, "ERROR: _obj should not be nullptr!");
   return *_obj;
 }

◆ 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::Utility::ends_with(), 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().

 {
   const bool gzipped_file = Utility::ends_with(name, ".gz");
 
   std::unique_ptr<std::istream> in;
 
   if (gzipped_file)
     {
 #ifdef LIBMESH_HAVE_GZSTREAM
       auto inf = std::make_unique<igzstream>();
       libmesh_assert(inf);
       inf->open(name.c_str(), std::ios::in);
       in = std::move(inf);
 #else
       libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
 #endif
     }
   else
     {
       auto inf = std::make_unique<std::ifstream>();
       libmesh_assert(inf);
 
       std::string new_name = Utility::unzip_file(name);
 
       inf->open(new_name.c_str(), std::ios::in);
       in = std::move(inf);
     }
 
   libmesh_assert(in.get());
 
   this->read_mesh (*in);
 }

◆ 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 174 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().

 {
   MeshBase & mesh = MeshInput<MeshBase>::mesh();
 
   // This is a serial-only process for now;
   // the Mesh should be read on processor 0 and
   // broadcast later
   libmesh_assert_equal_to (mesh.processor_id(), 0);
 
   libmesh_error_msg_if(!in.good(), "Can't read input stream");
 
   // clear any data in the mesh
   mesh.clear();
 
   // clear any of our own data
   spline_node_ptrs.clear();
   spline_nodeelem_ptrs.clear();
 
   // Expect different sections, in this order, perhaps with blank
   // lines and/or comments in between:
 
   enum FileSection {
     FILE_HEADER,
     PATCH_HEADER,
     NODE_LINES,
     // (repeat for each node)
     N_ELEM_SUBBLOCKS,
     ELEM_SUBBLOCK_HEADER,
     // (repeat for each subblock)
     ELEM_NODES_LINES,
     ELEM_COEF_VEC_IDS,
     // (repeat nodes lines + coef vec ids for each elem, subblock)
     N_COEF_BLOCKS, // number of coef vec blocks of each type
     N_VECS_PER_BLOCK, // number of coef vecs in each dense block
     COEF_VEC_COMPONENTS,
     //  (repeat coef vec components as necessary)
     //  (repeat coef blocks as necessary)
     //
     //  reserved for sparse block stuff we don't support yet
     END_OF_FILE };
 
   FileSection section = FILE_HEADER;
 
   // Values to remember from section to section
   dyna_int_type patch_id, n_spline_nodes, n_elem, n_coef_vec, weight_control_flag;
   dyna_int_type n_elem_blocks, n_dense_coef_vec_blocks;
   std::vector<dyna_int_type> // indexed from 0 to n_elem_blocks
     block_elem_type,
     block_n_elem,
     block_n_nodes,    // Number of *spline* nodes constraining elements
     block_n_coef_vec, // Number of coefficient vectors for each elem
     block_p,
     block_dim;
   std::vector<dyna_int_type> // indexed from 0 to n_dense_coef_vec_blocks
     n_coef_vecs_in_subblock, n_coef_comp;
   unsigned char weight_index = 0;
   dyna_int_type n_nodes_read = 0,
                 n_elem_blocks_read = 0,
                 n_elems_read = 0,
                 n_elem_nodes_read = 0,
                 n_elem_cvids_read = 0,
                 n_coef_headers_read = 0,
                 n_coef_blocks_read = 0,
                 n_coef_comp_read = 0,
                 n_coef_vecs_read = 0;
 
   // For reading the file line by line
   std::string s;
 
   // For storing global (spline) weights, until we have enough data to
   // use them for calculating local (Bezier element) nodes
   std::vector<Real> spline_weights;
 
   // For storing elements' constraint equations:
   // Global node indices (1-based in file, 0-based in memory):
   // elem_global_nodes[block_num][elem_num][local_index] = global_index
   std::vector<std::vector<std::vector<dof_id_type>>> elem_global_nodes;
 
   // Dense constraint vectors in the Dyna file
   // When first read:
   // dense_constraint_vecs[block_num][vec_in_block][column_num] = coef
   // When used:
   // dense_constraint_vecs[0][vec_num][column_num] = coef
   std::vector<std::vector<std::vector<Real>>> dense_constraint_vecs;
 
   // Constraint vector indices (1-based in file, 0-based in memory):
   // elem_constraint_rows[block_num][elem_num][row_num] = cv_index
   std::vector<std::vector<std::vector<dof_id_type>>> elem_constraint_rows;
 
   while (true)
     {
       // Try to read something.  This may set EOF!
       std::getline(in, s);
 
       if (in)
         {
           // Process s...
 
           if (s.find("B E X T 2.0") == static_cast<std::string::size_type>(0))
           {
             libmesh_error_msg_if(section != FILE_HEADER,
                                  "Found 'B E X T 2.0' outside file header?");
 
             section = PATCH_HEADER;
             continue;
           }
 
           // Ignore comments
           if (s.find("$#") == static_cast<std::string::size_type>(0))
             continue;
 
           // Ignore blank lines
           if (s.find_first_not_of(" \t") == std::string::npos)
             continue;
 
           // Parse each important section
           std::stringstream stream(s);
           switch (section) {
           case PATCH_HEADER:
             stream >> patch_id;
             stream >> n_spline_nodes;
             stream >> n_elem;
             stream >> n_coef_vec;
             stream >> weight_control_flag;
             libmesh_error_msg_if(stream.fail(), "Failure to parse patch header\n");
 
             spline_node_ptrs.resize(n_spline_nodes);
             spline_weights.resize(n_spline_nodes);
 
             // Even if we have w=1.0, we still want RATIONAL_BERNSTEIN
             // elements!
             // if (weight_control_flag)
               {
                 // If we ever add more nodes that aren't in this file,
                 // merge this mesh with a non-spline mesh, etc., 1.0
                 // is a good default weight
                 const Real default_weight = 1.0;
                 weight_index = cast_int<unsigned char>
                   (mesh.add_node_datum<Real>("rational_weight", true,
                                              &default_weight));
                 mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
                 mesh.set_default_mapping_data(weight_index);
               }
 
             section = NODE_LINES;
             break;
           case NODE_LINES:
             {
               std::array<dyna_fp_type, 4> xyzw;
               stream >> xyzw[0];
               stream >> xyzw[1];
               stream >> xyzw[2];
               stream >> xyzw[3];
 
               if (weight_control_flag)
                 spline_weights[n_nodes_read] = xyzw[3];
 
               // We'll use the spline nodes via NodeElem as the "true"
               // degrees of freedom, to which other Elem degrees of
               // freedom will be tied via constraint equations.
               Node *n = spline_node_ptrs[n_nodes_read] =
                 mesh.add_point(Point(xyzw[0], xyzw[1], xyzw[2]));
               Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
               elem->set_node(0, n);
               elem->subdomain_id() = 1; // Separate id to ease Exodus output
               spline_nodeelem_ptrs[n] = elem;
             }
             ++n_nodes_read;
 
             libmesh_error_msg_if(stream.fail(), "Failure to parse node line\n");
 
             if (n_nodes_read >= n_spline_nodes)
               section = N_ELEM_SUBBLOCKS;
             break;
           case N_ELEM_SUBBLOCKS:
             stream >> n_elem_blocks;
             libmesh_error_msg_if(stream.fail(), "Failure to parse n_elem_blocks\n");
 
             block_elem_type.resize(n_elem_blocks);
             block_n_elem.resize(n_elem_blocks);
             block_n_nodes.resize(n_elem_blocks);
             block_n_coef_vec.resize(n_elem_blocks);
             block_p.resize(n_elem_blocks);
             block_dim.resize(n_elem_blocks);
 
             elem_global_nodes.resize(n_elem_blocks);
             elem_constraint_rows.resize(n_elem_blocks);
 
             n_elem_blocks_read = 0;
             section = ELEM_SUBBLOCK_HEADER;
             break;
           case ELEM_SUBBLOCK_HEADER:
             stream >> block_elem_type[n_elem_blocks_read];
             stream >> block_n_elem[n_elem_blocks_read];
             stream >> block_n_nodes[n_elem_blocks_read];
             stream >> block_n_coef_vec[n_elem_blocks_read];
             stream >> block_p[n_elem_blocks_read];
 
             libmesh_error_msg_if(stream.fail(), "Failure to parse elem block\n");
 
             block_dim[n_elem_blocks_read] = 1; // All blocks here are at least 1D
 
             dyna_int_type block_other_p; // Check for isotropic p
             stream >> block_other_p;
             if (!stream.fail())
               {
                 block_dim[n_elem_blocks_read] = 2; // Found a second dimension!
 
                 if (block_other_p != block_p[n_elem_blocks_read])
                   libmesh_not_implemented(); // We don't support p anisotropy
 
                 stream >> block_other_p;
                 if (!stream.fail())
                   {
                     block_dim[n_elem_blocks_read] = 3; // Found a third dimension!
 
                     if (block_other_p != block_p[n_elem_blocks_read])
                       libmesh_not_implemented();
                   }
               }
 
             {
               auto & block_global_nodes = elem_global_nodes[n_elem_blocks_read];
               auto & block_constraint_rows = elem_constraint_rows[n_elem_blocks_read];
 
               block_global_nodes.resize(block_n_elem[n_elem_blocks_read]);
               block_constraint_rows.resize(block_n_elem[n_elem_blocks_read]);
 
               for (auto e : make_range(block_n_elem[n_elem_blocks_read]))
                 {
                   block_global_nodes[e].resize(block_n_nodes[n_elem_blocks_read]);
                   block_constraint_rows[e].resize(block_n_coef_vec[n_elem_blocks_read]);
                 }
             }
 
             n_elem_blocks_read++;
             if (n_elem_blocks_read == n_elem_blocks)
               {
                 n_elem_blocks_read = 0;
                 n_elems_read = 0;
                 section = ELEM_NODES_LINES;
               }
             break;
           case ELEM_NODES_LINES:
             {
               const int end_node_to_read =
                 std::min(block_n_nodes[n_elem_blocks_read], n_elem_nodes_read + max_ints_per_line);
               for (; n_elem_nodes_read != end_node_to_read; ++n_elem_nodes_read)
                 {
                   dyna_int_type node_id;
                   stream >> node_id;
                   node_id--;
                   elem_global_nodes[n_elem_blocks_read][n_elems_read][n_elem_nodes_read] = node_id;
 
                   // Let's assume that our *only* mid-line breaks are
                   // due to the max_ints_per_line limit.  This should be
                   // less flexible but better for error detection.
                   libmesh_error_msg_if(stream.fail(), "Failure to parse elem nodes\n");
                 }
 
               if (n_elem_nodes_read == block_n_nodes[n_elem_blocks_read])
                 {
                   n_elem_nodes_read = 0;
                   section = ELEM_COEF_VEC_IDS;
                 }
             }
             break;
           case ELEM_COEF_VEC_IDS:
             {
               const int end_cvid_to_read =
                 std::min(block_n_coef_vec[n_elem_blocks_read], n_elem_cvids_read + max_ints_per_line);
               for (; n_elem_cvids_read != end_cvid_to_read; ++n_elem_cvids_read)
                 {
                   dyna_int_type node_cvid;
                   stream >> node_cvid;
                   node_cvid--;
 
                   elem_constraint_rows[n_elem_blocks_read][n_elems_read][n_elem_cvids_read] = node_cvid;
 
                   // Let's assume that our *only* mid-line breaks are
                   // due to the max_ints_per_line limit.  This should be
                   // less flexible but better for error detection.
                   libmesh_error_msg_if(stream.fail(), "Failure to parse elem cvids\n");
                 }
               if (n_elem_cvids_read == block_n_nodes[n_elem_blocks_read])
                 {
                   n_elem_cvids_read = 0;
                   n_elems_read++;
                   section = ELEM_NODES_LINES; // Read another elem, nodes first
                   if (n_elems_read == block_n_elem[n_elem_blocks_read])
                     {
                       n_elems_read = 0;
                       n_elem_blocks_read++;
                       if (n_elem_blocks_read == n_elem_blocks)
                         section = N_COEF_BLOCKS; // Move on to coefficient vectors
                     }
                 }
             }
             break;
           case N_COEF_BLOCKS:
             {
               stream >> n_dense_coef_vec_blocks;
               dyna_int_type n_sparse_coef_vec_blocks;
               stream >> n_sparse_coef_vec_blocks;
 
               libmesh_error_msg_if(stream.fail(), "Failure to parse n_*_coef_vec_blocks\n");
 
               if (n_sparse_coef_vec_blocks != 0)
                 libmesh_not_implemented();
 
               dense_constraint_vecs.resize(n_dense_coef_vec_blocks);
               n_coef_vecs_in_subblock.resize(n_dense_coef_vec_blocks);
               n_coef_comp.resize(n_dense_coef_vec_blocks);
 
               section = N_VECS_PER_BLOCK;
             }
             break;
           case N_VECS_PER_BLOCK:
             stream >> n_coef_vecs_in_subblock[n_coef_headers_read];
             stream >> n_coef_comp[n_coef_headers_read];
 
             libmesh_error_msg_if(stream.fail(), "Failure to parse dense coef subblock header\n");
 
             dense_constraint_vecs[n_coef_headers_read].resize
               (n_coef_vecs_in_subblock[n_coef_headers_read]);
 
             for (auto & vec : dense_constraint_vecs[n_coef_headers_read])
               vec.resize(n_coef_comp[n_coef_headers_read]);
 
             n_coef_headers_read++;
             if (n_coef_headers_read == n_dense_coef_vec_blocks)
               {
                 n_coef_headers_read = 0;
                 section = COEF_VEC_COMPONENTS;
               }
             break;
           case COEF_VEC_COMPONENTS:
             {
               auto & current_vec =
                 dense_constraint_vecs[n_coef_blocks_read][n_coef_vecs_read];
 
               const int end_coef_to_read =
                 std::min(n_coef_comp[n_coef_blocks_read],
                          n_coef_comp_read + max_fps_per_line);
               for (; n_coef_comp_read != end_coef_to_read; ++n_coef_comp_read)
                 {
                   dyna_fp_type coef_comp;
                   stream >> coef_comp;
 
                   current_vec[n_coef_comp_read] = coef_comp;
 
                   // Let's assume that our *only* mid-line breaks are
                   // due to the max_fps_per_line limit.  This should be
                   // less flexible but better for error detection.
                   libmesh_error_msg_if(stream.fail(), "Failure to parse coefficients\n");
                 }
               if (n_coef_comp_read == n_coef_comp[n_coef_blocks_read])
                 {
                   n_coef_comp_read = 0;
                   n_coef_vecs_read++;
                   if (n_coef_vecs_read == n_coef_vecs_in_subblock[n_coef_blocks_read])
                     {
                       n_coef_vecs_read = 0;
                       n_coef_blocks_read++;
                       if (n_coef_blocks_read == n_dense_coef_vec_blocks)
                         section = END_OF_FILE;
                     }
                 }
             }
             break;
           default:
             libmesh_error();
           }
 
           if (section == END_OF_FILE)
             break;
         } // if (in)
       else if (in.eof())
         libmesh_error_msg("Premature end of file");
       else
         libmesh_error_msg("Input stream failure! Perhaps the file does not exist?");
     }
 
   // Merge dense_constraint_vecs blocks
   if (n_dense_coef_vec_blocks)
     for (auto coef_vec_block :
          IntRange<dyna_int_type>(1, n_dense_coef_vec_blocks))
       {
         auto & dcv0 = dense_constraint_vecs[0];
         auto & dcvi = dense_constraint_vecs[coef_vec_block];
         dcv0.insert(dcv0.end(),
                     std::make_move_iterator(dcvi.begin()),
                     std::make_move_iterator(dcvi.end()));
       }
   dense_constraint_vecs.resize(1);
 
   // Constraint matrices:
   // elem_constraint_mat[block_num][elem_num][local_node_index][elem_global_nodes_index] = c
   std::vector<std::vector<std::vector<std::vector<Real>>>> elem_constraint_mat(n_elem_blocks);
 
   // We need to calculate local nodes on the fly, and we'll be
   // calculating them from constraint matrix columns, and we'll need
   // to make sure that the same node is found each time it's
   // calculated from multiple neighboring elements.
   std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
 
   auto & constraint_rows = mesh.get_constraint_rows();
 
   for (auto block_num : make_range(n_elem_blocks))
     {
       elem_constraint_mat[block_num].resize(block_n_elem[block_num]);
 
       for (auto elem_num :
            make_range(block_n_elem[block_num]))
         {
           // Consult the import element table to determine which element to build
           const ElementDefinition & elem_defn =
             find_elem_definition(block_elem_type[block_num],
                                  block_dim[block_num],
                                  block_p[block_num]);
 
           auto elem = Elem::build(elem_defn.type);
           libmesh_error_msg_if(elem->dim() != block_dim[block_num],
                                "Elem dim " << elem->dim() <<
                                " != block_dim " << block_dim[block_num]);
 
           auto & my_constraint_rows = elem_constraint_rows[block_num][elem_num];
           auto & my_global_nodes    = elem_global_nodes[block_num][elem_num];
           auto & my_constraint_mat  = elem_constraint_mat[block_num][elem_num];
 
           my_constraint_mat.resize(block_n_coef_vec[block_num]);
           for (auto spline_node_index :
                make_range(block_n_coef_vec[block_num]))
             my_constraint_mat[spline_node_index].resize(elem->n_nodes());
 
           for (auto spline_node_index :
                make_range(block_n_coef_vec[block_num]))
             {
               // Find which coef block this elem's vectors are from
               const dyna_int_type elem_coef_vec_index =
                 my_constraint_rows[spline_node_index];
 
               dyna_int_type coef_block_num = 0;
               dyna_int_type first_block_coef_vec = 0;
               for (; elem_coef_vec_index >= first_block_coef_vec &&
                    coef_block_num != n_dense_coef_vec_blocks; ++coef_block_num)
                 {
                   first_block_coef_vec += n_coef_vecs_in_subblock[coef_block_num];
                 }
 
               // Make sure we did find a valid coef block
               libmesh_error_msg_if(coef_block_num == n_dense_coef_vec_blocks &&
                                    first_block_coef_vec <= elem_coef_vec_index,
                                    "Can't find valid constraint coef vector");
 
               coef_block_num--;
 
               libmesh_error_msg_if
                 (dyna_int_type(elem->n_nodes()) != n_coef_comp[coef_block_num],
                  "Found " << n_coef_comp[coef_block_num] <<
                  " constraint coef vectors for " <<
                  elem->n_nodes() << " nodes");
 
               for (auto elem_node_index :
                    make_range(elem->n_nodes()))
                 my_constraint_mat[spline_node_index][elem_node_index] =
                   dense_constraint_vecs[0][elem_coef_vec_index][elem_node_index];
             }
 
           for (auto elem_node_index :
                make_range(elem->n_nodes()))
             {
               dof_id_type global_node_idx = DofObject::invalid_id;
 
               // New finite element node data: dot product of
               // constraint matrix columns with spline node data.
               // Store that column as a key.
               std::vector<std::pair<dof_id_type, Real>> key;
 
               for (auto spline_node_index :
                    make_range(block_n_coef_vec[block_num]))
                 {
                   const dyna_int_type elem_coef_vec_index =
                     my_constraint_rows[spline_node_index];
 
                   const Real coef =
                     libmesh_vector_at(dense_constraint_vecs[0],
                                       elem_coef_vec_index)[elem_node_index];
 
                   // Global nodes are supposed to be in sorted order
                   if (global_node_idx != DofObject::invalid_id)
                     libmesh_error_msg_if(my_global_nodes[spline_node_index] <= global_node_idx,
                                          "Found unsorted global node");
 
                   global_node_idx = my_global_nodes[spline_node_index];
 
                   if (coef != 0) // Ignore irrelevant spline nodes
                     key.emplace_back(global_node_idx, coef);
                 }
 
               if (const auto local_node_it = local_nodes.find(key);
                   local_node_it != local_nodes.end())
                 elem->set_node(elem_defn.nodes[elem_node_index],
                   local_node_it->second);
               else
                 {
                   Point p(0);
                   Real w = 0;
                   std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
 
                   for (auto spline_node_index :
                        make_range(block_n_coef_vec[block_num]))
                     {
                       const dof_id_type my_node_idx =
                         my_global_nodes[spline_node_index];
 
                       const dyna_int_type elem_coef_vec_index =
                         my_constraint_rows[spline_node_index];
 
                       Node * spline_node =
                         libmesh_vector_at(spline_node_ptrs,
                                           my_node_idx);
 
                       const Real coef =
                         libmesh_vector_at(dense_constraint_vecs[0],
                                           elem_coef_vec_index)[elem_node_index];
                       p.add_scaled(*spline_node, coef);
                       w += coef * libmesh_vector_at(spline_weights,
                                                     my_node_idx);
 
                       // We don't need to store 0 entries;
                       // constraint_rows is a sparse structure.
                       if (coef)
                         {
                           Elem * nodeelem =
                             libmesh_map_find(spline_nodeelem_ptrs, spline_node);
                           constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
                         }
                     }
 
                   Node *n = mesh.add_point(p);
                   if (weight_control_flag)
                     n->set_extra_datum<Real>(weight_index, w);
                   local_nodes[key] = n;
                   elem->set_node(elem_defn.nodes[elem_node_index], n);
 
                   constraint_rows[n] = constraint_row;
                 }
             }
 
           mesh.add_elem(std::move(elem));
         }
     }
 
   if (!_keep_spline_nodes)
     MeshTools::clear_spline_nodes(MeshInput<MeshBase>::mesh());
 }

◆ 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; }

libMesh::MeshInput< MeshBase >::mesh

MeshBase & mesh()

Definition: mesh_input.h:178

libMesh::MeshBase::set_n_partitions

unsigned int & set_n_partitions()

Definition: mesh_base.h:1871

◆ 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().

 {
   char c, line[256];
 
   while (in.get(c), c==comment_start)
     in.getline (line, 255);
 
   // put back first character of
   // first non-comment line
   in.putback (c);
 }

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 213 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 165 of file dyna_io.h.

Referenced by read_mesh().

◆ elems_of_dimension

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

protectedinherited

A vector of bools describing what dimension elements have been encountered when reading a mesh.

Definition at line 107 of file mesh_input.h.

Referenced by libMesh::GMVIO::_read_one_cell(), libMesh::UNVIO::elements_in(), libMesh::UNVIO::max_elem_dimension_seen(), libMesh::AbaqusIO::max_elem_dimension_seen(), libMesh::AbaqusIO::read(), libMesh::GMVIO::read(), libMesh::Nemesis_IO::read(), libMesh::ExodusII_IO::read(), libMesh::VTKIO::read(), libMesh::AbaqusIO::read_elements(), libMesh::UCDIO::read_implementation(), libMesh::UNVIO::read_implementation(), and libMesh::XdrIO::read_serialized_connectivity().

◆ 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 187 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 177 of file dyna_io.h.

Referenced by read_mesh().

◆ spline_node_ptrs

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

private

Definition at line 159 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 160 of file dyna_io.h.

Referenced by read_mesh().

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

include/mesh/dyna_io.h
src/mesh/dyna_io.C

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

Detailed Description

Member Typedef Documentation

◆ dyna_fp_type

◆ dyna_int_type

Constructor & Destructor Documentation

◆ DynaIO()

Member Function Documentation

◆ add_spline_constraints()

◆ build_element_maps()

◆ clear_spline_nodes()

◆ find_elem_definition() [1/2]

◆ find_elem_definition() [2/2]

◆ is_parallel_format()

◆ mesh()

◆ read()

◆ read_mesh()

◆ set_n_partitions()

◆ skip_comment_lines()

Member Data Documentation

◆ _element_maps

◆ _keep_spline_nodes

◆ elems_of_dimension

◆ max_fps_per_line

◆ max_ints_per_line

◆ spline_node_ptrs

◆ spline_nodeelem_ptrs