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:

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.

Public Member Functions
	DynaIO (MeshBase &mesh, bool keep_spline_nodes=true)
	Constructor.
virtual void	read (const std::string &name) override
	Reads in a mesh in the Dyna format from the ASCII file given by name.
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.
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.

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

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

Private Types
typedef double	dyna_fp_type
	The floating-point type DYNA uses.

Private Member Functions
void	read_mesh (std::istream &in)
	Implementation of the read() function.

Static Private Member Functions
static ElementMaps	build_element_maps ()
	A static function used to construct the _element_maps struct, statically.

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.
MeshBase *	_obj
	A pointer to a non-const object object.
const bool	_is_parallel_format
	Flag specifying whether this format is parallel-capable.

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

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

dyna_int_type

typedef int32_t libMesh::DynaIO::dyna_int_type

The integer type DYNA uses.

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

{
}

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.

{
  // 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;
}

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

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

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

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

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

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.

{ return this->_is_parallel_format; }

mesh()

MeshBase & libMesh::MeshInput< MeshBase >::mesh ( )

inlineprotectedinherited

Returns

The object as a writable reference.

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

{
  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);
}

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

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

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.

{
  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());
}

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

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.

{ this->mesh().set_n_partitions() = n_parts; }

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 113 of file mesh_input.h.

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

Referenced by find_elem_definition().

_is_parallel_format

const bool libMesh::MeshInput< MeshBase >::_is_parallel_format

privateinherited

Flag specifying whether this format is parallel-capable.

If this is false (default) I/O is only permitted when the mesh has been serialized.

Definition at line 130 of file mesh_input.h.

_keep_spline_nodes

bool libMesh::DynaIO::_keep_spline_nodes

private

Whether to keep or eventually discard spline nodes.

Definition at line 154 of file dyna_io.h.

Referenced by read_mesh().

_obj

MeshBase * libMesh::MeshInput< MeshBase >::_obj

privateinherited

A pointer to a non-const object object.

This allows us to read the object from file.

Definition at line 123 of file mesh_input.h.

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.

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

Referenced by read_mesh().

spline_node_ptrs

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

private

Definition at line 148 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 149 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