libMesh
dyna_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // Local includes
19 #include "libmesh/libmesh_config.h"
20 #include "libmesh/libmesh_logging.h"
21 #include "libmesh/dof_map.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/dyna_io.h"
24 #include "libmesh/mesh_base.h"
25 #include "libmesh/mesh_tools.h"
26 #include "libmesh/int_range.h"
27 #include "libmesh/utility.h"
28 
29 // TIMPI includes
31 
32 // gzstream for reading compressed files as a stream
33 #ifdef LIBMESH_HAVE_GZSTREAM
34 # include "gzstream.h"
35 #endif
36 
37 // C++ includes
38 #include <array>
39 #include <cstddef>
40 #include <fstream>
41 #include <iterator>
42 #include <numeric> // iota
43 
44 namespace libMesh
45 {
46 
47 // Initialize the static data member
49 
50 
51 
52 // Definition of the static function which constructs the ElementMaps object.
54 {
55  // Object to be filled up
56  ElementMaps em;
57 
58  // em.add_def(ElementDefinition(TRI3, 2, 2, 1)); // node mapping?
59  // em.add_def(ElementDefinition(TRI6, 2, 2, 2)); // node mapping?
60  // em.add_def(ElementDefinition(TET4, 2, 3, 1)); // node mapping?
61  // em.add_def(ElementDefinition(TET10, 2, 3, 2)); // node mapping?
62  // em.add_def(ElementDefinition(PRISM6, 3, 3, 1)); // node mapping?
63  // em.add_def(ElementDefinition(PRISM18, 3, 3, 2)); // node mapping?
64 
65  // Eventually we'll map both tensor-product and non-tensor-product
66  // BEXT elements to ours; for now we only support non-tensor-product
67  // Bezier elements.
68  // for (unsigned int i=0; i != 2; ++i)
69  for (unsigned int i=1; i != 2; ++i)
70  {
71  // We only have one element for whom node orders match...
72  em.add_def(ElementDefinition(EDGE2, i, 1, 1));
73 
74  em.add_def(ElementDefinition(EDGE3, i, 1, 2,
75  {0, 2, 1}));
76  em.add_def(ElementDefinition(EDGE4, i, 1, 2,
77  {0, 2, 3, 1}));
78 
79  em.add_def(ElementDefinition(QUAD4, i, 2, 1,
80  {0, 1, 3, 2}));
81  em.add_def(ElementDefinition(QUAD9, i, 2, 2,
82  {0, 4, 1, 7, 8, 5, 3, 6, 2}));
83 
84  em.add_def(ElementDefinition(HEX8, i, 3, 1,
85  {0, 1, 3, 2, 4, 5, 7, 6}));
86  em.add_def(ElementDefinition(HEX27, i, 3, 2,
87  { 0, 8, 1, 11, 20, 9, 3, 10, 2,
88  12, 21, 13, 24, 26, 22, 15, 23, 14,
89  4, 16, 5, 19, 25, 17, 7, 18, 6}));
90  }
91 
92  return em;
93 }
94 
95 
96 
98  (ElemType type_in,
99  dyna_int_type dyna_type_in,
100  dyna_int_type dim_in,
101  dyna_int_type p_in) :
102  type(type_in),
103  dyna_type(dyna_type_in),
104  dim(dim_in),
105  p(p_in)
106 {
107  const unsigned int n_nodes = Elem::type_to_n_nodes_map[type_in];
108  if (n_nodes == invalid_uint)
109  libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
110  nodes.resize(n_nodes);
111  std::iota(nodes.begin(), nodes.end(), 0);
112 }
113 
114 
116  (ElemType type_in,
117  dyna_int_type dyna_type_in,
118  dyna_int_type dim_in,
119  dyna_int_type p_in,
120  std::vector<unsigned int> && nodes_in) :
121  type(type_in),
122  dyna_type(dyna_type_in),
123  dim(dim_in),
124  p(p_in),
125  nodes(nodes_in)
126 {}
127 
128 
129 
131  bool keep_spline_nodes) :
133  _keep_spline_nodes (keep_spline_nodes)
134 {
135 }
136 
137 
138 
139 void DynaIO::read (const std::string & name)
140 {
141  const bool gzipped_file = Utility::ends_with(name, ".gz");
142 
143  std::unique_ptr<std::istream> in;
144 
145  if (gzipped_file)
146  {
147 #ifdef LIBMESH_HAVE_GZSTREAM
148  auto inf = std::make_unique<igzstream>();
149  libmesh_assert(inf);
150  inf->open(name.c_str(), std::ios::in);
151  in = std::move(inf);
152 #else
153  libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
154 #endif
155  }
156  else
157  {
158  auto inf = std::make_unique<std::ifstream>();
159  libmesh_assert(inf);
160 
161  std::string new_name = Utility::unzip_file(name);
162 
163  inf->open(new_name.c_str(), std::ios::in);
164  in = std::move(inf);
165  }
166 
167  libmesh_assert(in.get());
168 
169  this->read_mesh (*in);
170 }
171 
172 
173 
174 void DynaIO::read_mesh(std::istream & in)
175 {
177 
178  // This is a serial-only process for now;
179  // the Mesh should be read on processor 0 and
180  // broadcast later
181  libmesh_assert_equal_to (mesh.processor_id(), 0);
182 
183  libmesh_error_msg_if(!in.good(), "Can't read input stream");
184 
185  // clear any data in the mesh
186  mesh.clear();
187 
188  // clear any of our own data
189  spline_node_ptrs.clear();
190  spline_nodeelem_ptrs.clear();
191 
192  // Expect different sections, in this order, perhaps with blank
193  // lines and/or comments in between:
194 
195  enum FileSection {
196  FILE_HEADER,
197  PATCH_HEADER,
198  NODE_LINES,
199  // (repeat for each node)
200  N_ELEM_SUBBLOCKS,
201  ELEM_SUBBLOCK_HEADER,
202  // (repeat for each subblock)
203  ELEM_NODES_LINES,
204  ELEM_COEF_VEC_IDS,
205  // (repeat nodes lines + coef vec ids for each elem, subblock)
206  N_COEF_BLOCKS, // number of coef vec blocks of each type
207  N_VECS_PER_BLOCK, // number of coef vecs in each dense block
208  COEF_VEC_COMPONENTS,
209  // (repeat coef vec components as necessary)
210  // (repeat coef blocks as necessary)
211  //
212  // reserved for sparse block stuff we don't support yet
213  END_OF_FILE };
214 
215  FileSection section = FILE_HEADER;
216 
217  // Values to remember from section to section
218  dyna_int_type patch_id, n_spline_nodes, n_elem, n_coef_vec, weight_control_flag;
219  dyna_int_type n_elem_blocks, n_dense_coef_vec_blocks;
220  std::vector<dyna_int_type> // indexed from 0 to n_elem_blocks
221  block_elem_type,
222  block_n_elem,
223  block_n_nodes, // Number of *spline* nodes constraining elements
224  block_n_coef_vec, // Number of coefficient vectors for each elem
225  block_p,
226  block_dim;
227  std::vector<dyna_int_type> // indexed from 0 to n_dense_coef_vec_blocks
228  n_coef_vecs_in_subblock, n_coef_comp;
229  unsigned char weight_index = 0;
230  dyna_int_type n_nodes_read = 0,
231  n_elem_blocks_read = 0,
232  n_elems_read = 0,
233  n_elem_nodes_read = 0,
234  n_elem_cvids_read = 0,
235  n_coef_headers_read = 0,
236  n_coef_blocks_read = 0,
237  n_coef_comp_read = 0,
238  n_coef_vecs_read = 0;
239 
240  // For reading the file line by line
241  std::string s;
242 
243  // For storing global (spline) weights, until we have enough data to
244  // use them for calculating local (Bezier element) nodes
245  std::vector<Real> spline_weights;
246 
247  // For storing elements' constraint equations:
248  // Global node indices (1-based in file, 0-based in memory):
249  // elem_global_nodes[block_num][elem_num][local_index] = global_index
250  std::vector<std::vector<std::vector<dof_id_type>>> elem_global_nodes;
251 
252  // Dense constraint vectors in the Dyna file
253  // When first read:
254  // dense_constraint_vecs[block_num][vec_in_block][column_num] = coef
255  // When used:
256  // dense_constraint_vecs[0][vec_num][column_num] = coef
257  std::vector<std::vector<std::vector<Real>>> dense_constraint_vecs;
258 
259  // Constraint vector indices (1-based in file, 0-based in memory):
260  // elem_constraint_rows[block_num][elem_num][row_num] = cv_index
261  std::vector<std::vector<std::vector<dof_id_type>>> elem_constraint_rows;
262 
263  while (true)
264  {
265  // Try to read something. This may set EOF!
266  std::getline(in, s);
267 
268  if (in)
269  {
270  // Process s...
271 
272  if (s.find("B E X T 2.0") == static_cast<std::string::size_type>(0))
273  {
274  libmesh_error_msg_if(section != FILE_HEADER,
275  "Found 'B E X T 2.0' outside file header?");
276 
277  section = PATCH_HEADER;
278  continue;
279  }
280 
281  // Ignore comments
282  if (s.find("$#") == static_cast<std::string::size_type>(0))
283  continue;
284 
285  // Ignore blank lines
286  if (s.find_first_not_of(" \t") == std::string::npos)
287  continue;
288 
289  // Parse each important section
290  std::stringstream stream(s);
291  switch (section) {
292  case PATCH_HEADER:
293  stream >> patch_id;
294  stream >> n_spline_nodes;
295  stream >> n_elem;
296  stream >> n_coef_vec;
297  stream >> weight_control_flag;
298  libmesh_error_msg_if(stream.fail(), "Failure to parse patch header\n");
299 
300  spline_node_ptrs.resize(n_spline_nodes);
301  spline_weights.resize(n_spline_nodes);
302 
303  // Even if we have w=1.0, we still want RATIONAL_BERNSTEIN
304  // elements!
305  // if (weight_control_flag)
306  {
307  // If we ever add more nodes that aren't in this file,
308  // merge this mesh with a non-spline mesh, etc., 1.0
309  // is a good default weight
310  const Real default_weight = 1.0;
311  weight_index = cast_int<unsigned char>
312  (mesh.add_node_datum<Real>("rational_weight", true,
313  &default_weight));
315  mesh.set_default_mapping_data(weight_index);
316  }
317 
318  section = NODE_LINES;
319  break;
320  case NODE_LINES:
321  {
322  std::array<dyna_fp_type, 4> xyzw;
323  stream >> xyzw[0];
324  stream >> xyzw[1];
325  stream >> xyzw[2];
326  stream >> xyzw[3];
327 
328  if (weight_control_flag)
329  spline_weights[n_nodes_read] = xyzw[3];
330 
331  // We'll use the spline nodes via NodeElem as the "true"
332  // degrees of freedom, to which other Elem degrees of
333  // freedom will be tied via constraint equations.
334  Node *n = spline_node_ptrs[n_nodes_read] =
335  mesh.add_point(Point(xyzw[0], xyzw[1], xyzw[2]));
337  elem->set_node(0, n);
338  elem->subdomain_id() = 1; // Separate id to ease Exodus output
339  spline_nodeelem_ptrs[n] = elem;
340  }
341  ++n_nodes_read;
342 
343  libmesh_error_msg_if(stream.fail(), "Failure to parse node line\n");
344 
345  if (n_nodes_read >= n_spline_nodes)
346  section = N_ELEM_SUBBLOCKS;
347  break;
348  case N_ELEM_SUBBLOCKS:
349  stream >> n_elem_blocks;
350  libmesh_error_msg_if(stream.fail(), "Failure to parse n_elem_blocks\n");
351 
352  block_elem_type.resize(n_elem_blocks);
353  block_n_elem.resize(n_elem_blocks);
354  block_n_nodes.resize(n_elem_blocks);
355  block_n_coef_vec.resize(n_elem_blocks);
356  block_p.resize(n_elem_blocks);
357  block_dim.resize(n_elem_blocks);
358 
359  elem_global_nodes.resize(n_elem_blocks);
360  elem_constraint_rows.resize(n_elem_blocks);
361 
362  n_elem_blocks_read = 0;
363  section = ELEM_SUBBLOCK_HEADER;
364  break;
365  case ELEM_SUBBLOCK_HEADER:
366  stream >> block_elem_type[n_elem_blocks_read];
367  stream >> block_n_elem[n_elem_blocks_read];
368  stream >> block_n_nodes[n_elem_blocks_read];
369  stream >> block_n_coef_vec[n_elem_blocks_read];
370  stream >> block_p[n_elem_blocks_read];
371 
372  libmesh_error_msg_if(stream.fail(), "Failure to parse elem block\n");
373 
374  block_dim[n_elem_blocks_read] = 1; // All blocks here are at least 1D
375 
376  dyna_int_type block_other_p; // Check for isotropic p
377  stream >> block_other_p;
378  if (!stream.fail())
379  {
380  block_dim[n_elem_blocks_read] = 2; // Found a second dimension!
381 
382  if (block_other_p != block_p[n_elem_blocks_read])
383  libmesh_not_implemented(); // We don't support p anisotropy
384 
385  stream >> block_other_p;
386  if (!stream.fail())
387  {
388  block_dim[n_elem_blocks_read] = 3; // Found a third dimension!
389 
390  if (block_other_p != block_p[n_elem_blocks_read])
391  libmesh_not_implemented();
392  }
393  }
394 
395  {
396  auto & block_global_nodes = elem_global_nodes[n_elem_blocks_read];
397  auto & block_constraint_rows = elem_constraint_rows[n_elem_blocks_read];
398 
399  block_global_nodes.resize(block_n_elem[n_elem_blocks_read]);
400  block_constraint_rows.resize(block_n_elem[n_elem_blocks_read]);
401 
402  for (auto e : make_range(block_n_elem[n_elem_blocks_read]))
403  {
404  block_global_nodes[e].resize(block_n_nodes[n_elem_blocks_read]);
405  block_constraint_rows[e].resize(block_n_coef_vec[n_elem_blocks_read]);
406  }
407  }
408 
409  n_elem_blocks_read++;
410  if (n_elem_blocks_read == n_elem_blocks)
411  {
412  n_elem_blocks_read = 0;
413  n_elems_read = 0;
414  section = ELEM_NODES_LINES;
415  }
416  break;
417  case ELEM_NODES_LINES:
418  {
419  const int end_node_to_read =
420  std::min(block_n_nodes[n_elem_blocks_read], n_elem_nodes_read + max_ints_per_line);
421  for (; n_elem_nodes_read != end_node_to_read; ++n_elem_nodes_read)
422  {
423  dyna_int_type node_id;
424  stream >> node_id;
425  node_id--;
426  elem_global_nodes[n_elem_blocks_read][n_elems_read][n_elem_nodes_read] = node_id;
427 
428  // Let's assume that our *only* mid-line breaks are
429  // due to the max_ints_per_line limit. This should be
430  // less flexible but better for error detection.
431  libmesh_error_msg_if(stream.fail(), "Failure to parse elem nodes\n");
432  }
433 
434  if (n_elem_nodes_read == block_n_nodes[n_elem_blocks_read])
435  {
436  n_elem_nodes_read = 0;
437  section = ELEM_COEF_VEC_IDS;
438  }
439  }
440  break;
441  case ELEM_COEF_VEC_IDS:
442  {
443  const int end_cvid_to_read =
444  std::min(block_n_coef_vec[n_elem_blocks_read], n_elem_cvids_read + max_ints_per_line);
445  for (; n_elem_cvids_read != end_cvid_to_read; ++n_elem_cvids_read)
446  {
447  dyna_int_type node_cvid;
448  stream >> node_cvid;
449  node_cvid--;
450 
451  elem_constraint_rows[n_elem_blocks_read][n_elems_read][n_elem_cvids_read] = node_cvid;
452 
453  // Let's assume that our *only* mid-line breaks are
454  // due to the max_ints_per_line limit. This should be
455  // less flexible but better for error detection.
456  libmesh_error_msg_if(stream.fail(), "Failure to parse elem cvids\n");
457  }
458  if (n_elem_cvids_read == block_n_nodes[n_elem_blocks_read])
459  {
460  n_elem_cvids_read = 0;
461  n_elems_read++;
462  section = ELEM_NODES_LINES; // Read another elem, nodes first
463  if (n_elems_read == block_n_elem[n_elem_blocks_read])
464  {
465  n_elems_read = 0;
466  n_elem_blocks_read++;
467  if (n_elem_blocks_read == n_elem_blocks)
468  section = N_COEF_BLOCKS; // Move on to coefficient vectors
469  }
470  }
471  }
472  break;
473  case N_COEF_BLOCKS:
474  {
475  stream >> n_dense_coef_vec_blocks;
476  dyna_int_type n_sparse_coef_vec_blocks;
477  stream >> n_sparse_coef_vec_blocks;
478 
479  libmesh_error_msg_if(stream.fail(), "Failure to parse n_*_coef_vec_blocks\n");
480 
481  if (n_sparse_coef_vec_blocks != 0)
482  libmesh_not_implemented();
483 
484  dense_constraint_vecs.resize(n_dense_coef_vec_blocks);
485  n_coef_vecs_in_subblock.resize(n_dense_coef_vec_blocks);
486  n_coef_comp.resize(n_dense_coef_vec_blocks);
487 
488  section = N_VECS_PER_BLOCK;
489  }
490  break;
491  case N_VECS_PER_BLOCK:
492  stream >> n_coef_vecs_in_subblock[n_coef_headers_read];
493  stream >> n_coef_comp[n_coef_headers_read];
494 
495  libmesh_error_msg_if(stream.fail(), "Failure to parse dense coef subblock header\n");
496 
497  dense_constraint_vecs[n_coef_headers_read].resize
498  (n_coef_vecs_in_subblock[n_coef_headers_read]);
499 
500  for (auto & vec : dense_constraint_vecs[n_coef_headers_read])
501  vec.resize(n_coef_comp[n_coef_headers_read]);
502 
503  n_coef_headers_read++;
504  if (n_coef_headers_read == n_dense_coef_vec_blocks)
505  {
506  n_coef_headers_read = 0;
507  section = COEF_VEC_COMPONENTS;
508  }
509  break;
510  case COEF_VEC_COMPONENTS:
511  {
512  auto & current_vec =
513  dense_constraint_vecs[n_coef_blocks_read][n_coef_vecs_read];
514 
515  const int end_coef_to_read =
516  std::min(n_coef_comp[n_coef_blocks_read],
517  n_coef_comp_read + max_fps_per_line);
518  for (; n_coef_comp_read != end_coef_to_read; ++n_coef_comp_read)
519  {
520  dyna_fp_type coef_comp;
521  stream >> coef_comp;
522 
523  current_vec[n_coef_comp_read] = coef_comp;
524 
525  // Let's assume that our *only* mid-line breaks are
526  // due to the max_fps_per_line limit. This should be
527  // less flexible but better for error detection.
528  libmesh_error_msg_if(stream.fail(), "Failure to parse coefficients\n");
529  }
530  if (n_coef_comp_read == n_coef_comp[n_coef_blocks_read])
531  {
532  n_coef_comp_read = 0;
533  n_coef_vecs_read++;
534  if (n_coef_vecs_read == n_coef_vecs_in_subblock[n_coef_blocks_read])
535  {
536  n_coef_vecs_read = 0;
537  n_coef_blocks_read++;
538  if (n_coef_blocks_read == n_dense_coef_vec_blocks)
539  section = END_OF_FILE;
540  }
541  }
542  }
543  break;
544  default:
545  libmesh_error();
546  }
547 
548  if (section == END_OF_FILE)
549  break;
550  } // if (in)
551  else if (in.eof())
552  libmesh_error_msg("Premature end of file");
553  else
554  libmesh_error_msg("Input stream failure! Perhaps the file does not exist?");
555  }
556 
557  // Merge dense_constraint_vecs blocks
558  if (n_dense_coef_vec_blocks)
559  for (auto coef_vec_block :
560  IntRange<dyna_int_type>(1, n_dense_coef_vec_blocks))
561  {
562  auto & dcv0 = dense_constraint_vecs[0];
563  auto & dcvi = dense_constraint_vecs[coef_vec_block];
564  dcv0.insert(dcv0.end(),
565  std::make_move_iterator(dcvi.begin()),
566  std::make_move_iterator(dcvi.end()));
567  }
568  dense_constraint_vecs.resize(1);
569 
570  // Constraint matrices:
571  // elem_constraint_mat[block_num][elem_num][local_node_index][elem_global_nodes_index] = c
572  std::vector<std::vector<std::vector<std::vector<Real>>>> elem_constraint_mat(n_elem_blocks);
573 
574  // We need to calculate local nodes on the fly, and we'll be
575  // calculating them from constraint matrix columns, and we'll need
576  // to make sure that the same node is found each time it's
577  // calculated from multiple neighboring elements.
578  std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
579 
580  auto & constraint_rows = mesh.get_constraint_rows();
581 
582  for (auto block_num : make_range(n_elem_blocks))
583  {
584  elem_constraint_mat[block_num].resize(block_n_elem[block_num]);
585 
586  for (auto elem_num :
587  make_range(block_n_elem[block_num]))
588  {
589  // Consult the import element table to determine which element to build
590  const ElementDefinition & elem_defn =
591  find_elem_definition(block_elem_type[block_num],
592  block_dim[block_num],
593  block_p[block_num]);
594 
595  auto elem = Elem::build(elem_defn.type);
596  libmesh_error_msg_if(elem->dim() != block_dim[block_num],
597  "Elem dim " << elem->dim() <<
598  " != block_dim " << block_dim[block_num]);
599 
600  auto & my_constraint_rows = elem_constraint_rows[block_num][elem_num];
601  auto & my_global_nodes = elem_global_nodes[block_num][elem_num];
602  auto & my_constraint_mat = elem_constraint_mat[block_num][elem_num];
603 
604  my_constraint_mat.resize(block_n_coef_vec[block_num]);
605  for (auto spline_node_index :
606  make_range(block_n_coef_vec[block_num]))
607  my_constraint_mat[spline_node_index].resize(elem->n_nodes());
608 
609  for (auto spline_node_index :
610  make_range(block_n_coef_vec[block_num]))
611  {
612  // Find which coef block this elem's vectors are from
613  const dyna_int_type elem_coef_vec_index =
614  my_constraint_rows[spline_node_index];
615 
616  dyna_int_type coef_block_num = 0;
617  dyna_int_type first_block_coef_vec = 0;
618  for (; elem_coef_vec_index >= first_block_coef_vec &&
619  coef_block_num != n_dense_coef_vec_blocks; ++coef_block_num)
620  {
621  first_block_coef_vec += n_coef_vecs_in_subblock[coef_block_num];
622  }
623 
624  // Make sure we did find a valid coef block
625  libmesh_error_msg_if(coef_block_num == n_dense_coef_vec_blocks &&
626  first_block_coef_vec <= elem_coef_vec_index,
627  "Can't find valid constraint coef vector");
628 
629  coef_block_num--;
630 
631  libmesh_error_msg_if
632  (dyna_int_type(elem->n_nodes()) != n_coef_comp[coef_block_num],
633  "Found " << n_coef_comp[coef_block_num] <<
634  " constraint coef vectors for " <<
635  elem->n_nodes() << " nodes");
636 
637  for (auto elem_node_index :
638  make_range(elem->n_nodes()))
639  my_constraint_mat[spline_node_index][elem_node_index] =
640  dense_constraint_vecs[0][elem_coef_vec_index][elem_node_index];
641  }
642 
643  for (auto elem_node_index :
644  make_range(elem->n_nodes()))
645  {
646  dof_id_type global_node_idx = DofObject::invalid_id;
647 
648  // New finite element node data: dot product of
649  // constraint matrix columns with spline node data.
650  // Store that column as a key.
651  std::vector<std::pair<dof_id_type, Real>> key;
652 
653  for (auto spline_node_index :
654  make_range(block_n_coef_vec[block_num]))
655  {
656  const dyna_int_type elem_coef_vec_index =
657  my_constraint_rows[spline_node_index];
658 
659  const Real coef =
660  libmesh_vector_at(dense_constraint_vecs[0],
661  elem_coef_vec_index)[elem_node_index];
662 
663  // Global nodes are supposed to be in sorted order
664  if (global_node_idx != DofObject::invalid_id)
665  libmesh_error_msg_if(my_global_nodes[spline_node_index] <= global_node_idx,
666  "Found unsorted global node");
667 
668  global_node_idx = my_global_nodes[spline_node_index];
669 
670  if (coef != 0) // Ignore irrelevant spline nodes
671  key.emplace_back(global_node_idx, coef);
672  }
673 
674  if (const auto local_node_it = local_nodes.find(key);
675  local_node_it != local_nodes.end())
676  elem->set_node(elem_defn.nodes[elem_node_index],
677  local_node_it->second);
678  else
679  {
680  Point p(0);
681  Real w = 0;
682  std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
683 
684  for (auto spline_node_index :
685  make_range(block_n_coef_vec[block_num]))
686  {
687  const dof_id_type my_node_idx =
688  my_global_nodes[spline_node_index];
689 
690  const dyna_int_type elem_coef_vec_index =
691  my_constraint_rows[spline_node_index];
692 
693  Node * spline_node =
694  libmesh_vector_at(spline_node_ptrs,
695  my_node_idx);
696 
697  const Real coef =
698  libmesh_vector_at(dense_constraint_vecs[0],
699  elem_coef_vec_index)[elem_node_index];
700  p.add_scaled(*spline_node, coef);
701  w += coef * libmesh_vector_at(spline_weights,
702  my_node_idx);
703 
704  // We don't need to store 0 entries;
705  // constraint_rows is a sparse structure.
706  if (coef)
707  {
708  Elem * nodeelem =
709  libmesh_map_find(spline_nodeelem_ptrs, spline_node);
710  constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
711  }
712  }
713 
714  Node *n = mesh.add_point(p);
715  if (weight_control_flag)
716  n->set_extra_datum<Real>(weight_index, w);
717  local_nodes[key] = n;
718  elem->set_node(elem_defn.nodes[elem_node_index], n);
719 
720  constraint_rows[n] = constraint_row;
721  }
722  }
723 
724  mesh.add_elem(std::move(elem));
725  }
726  }
727 
728  if (!_keep_spline_nodes)
730 }
731 
732 
734  unsigned int,
735  unsigned int)
736 {
737 }
738 
739 
740 #ifdef LIBMESH_ENABLE_DEPRECATED
742 {
743  libmesh_deprecated();
744 
746 }
747 #endif // LIBMESH_ENABLE_DEPRECATED
748 
749 
750 
753  int dim, int p)
754 {
755  auto eletypes_it =
756  _element_maps.in.find(std::make_tuple(dyna_elem, dim, p));
757 
758  // Make sure we actually found something
759  libmesh_error_msg_if
760  (eletypes_it == _element_maps.in.end(),
761  "Element of type " << dyna_elem <<
762  " dim " << dim <<
763  " degree " << p << " not found!");
764 
765  return eletypes_it->second;
766 }
767 
768 
769 
772  int libmesh_dbg_var(dim),
773  int libmesh_dbg_var(p))
774 {
775  auto eletypes_it =
776  _element_maps.out.find(libmesh_elem);
777 
778  // Make sure we actually found something
779  libmesh_error_msg_if
780  (eletypes_it == _element_maps.out.end(),
781  "Element of type " << libmesh_elem <<
782  " not found!");
783 
784  // Make sure we found the right thing
785  libmesh_assert_equal_to(eletypes_it->second.dim, dim);
786  libmesh_assert_equal_to(eletypes_it->second.p, p);
787 
788 
789  return eletypes_it->second;
790 }
791 
792 
793 
794 
795 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
ElemType
Defines an enum for geometric element types.
virtual void read(const std::string &name) override
Reads in a mesh in the Dyna format from the ASCII file given by name.
Definition: dyna_io.C:139
std::vector< Node * > spline_node_ptrs
Definition: dyna_io.h:159
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
Definition: mesh_base.h:1709
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
A Node is like a Point, but with more information.
Definition: node.h:52
int32_t dyna_int_type
The integer type DYNA uses.
Definition: dyna_io.h:116
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
double dyna_fp_type
The floating-point type DYNA uses.
Definition: dyna_io.h:182
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:969
bool ends_with(std::string_view superstring, std::string_view suffix)
Look for a substring at the very end of a string.
Definition: utility.C:213
unsigned int dim
void read_mesh(std::istream &in)
Implementation of the read() function.
Definition: dyna_io.C:174
void clear_spline_nodes(MeshBase &)
Remove spline node (for IsoGeometric Analysis meshes) elements and nodes and constraints from the mes...
Definition: mesh_tools.C:1147
Defines mapping from libMesh element types to LS-DYNA element types or vice-versa.
Definition: dyna_io.h:125
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void clear_spline_nodes()
Removes any spline nodes (both NodeElem and Node), leaving only the FE mesh generated from those spli...
Definition: dyna_io.C:741
void add_def(const ElementDefinition &eledef)
Definition: dyna_io.h:196
bool _keep_spline_nodes
Whether to keep or eventually discard spline nodes.
Definition: dyna_io.h:165
The libMesh namespace provides an interface to certain functionality in the library.
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
static const unsigned int type_to_n_nodes_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of nodes in the element...
Definition: elem.h:650
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
This is the MeshBase class.
Definition: mesh_base.h:75
void add_spline_constraints(DofMap &dof_map, unsigned int sys_num, unsigned int var_num)
Constrains finite element degrees of freedom in terms of spline degrees of freedom by adding user-def...
Definition: dyna_io.C:733
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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:752
const dof_id_type n_nodes
Definition: tecplot_io.C:67
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:57
std::map< ElemType, ElementDefinition > out
Definition: dyna_io.h:202
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:444
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:821
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:164
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
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:1126
static ElementMaps build_element_maps()
A static function used to construct the _element_maps struct, statically.
Definition: dyna_io.C:53
static const int max_fps_per_line
How many can we find on a line?
Definition: dyna_io.h:187
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
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:2347
DynaIO(MeshBase &mesh, bool keep_spline_nodes=true)
Constructor.
Definition: dyna_io.C:130
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
std::unordered_map< Node *, Elem * > spline_nodeelem_ptrs
Definition: dyna_io.h:160
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
Definition: dyna_io.h:213
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:140
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:839
static const int max_ints_per_line
How many can we find on a line?
Definition: dyna_io.h:177
struct which holds a map from LS-DYNA to libMesh element numberings and vice-versa.
Definition: dyna_io.h:193
std::vector< unsigned int > nodes
Definition: dyna_io.h:142
processor_id_type processor_id() const
std::map< Key, ElementDefinition > in
Definition: dyna_io.h:206
ElementDefinition(ElemType type_in, dyna_int_type dyna_type_in, dyna_int_type dim_in, dyna_int_type p_in)
Definition: dyna_io.C:98
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
uint8_t dof_id_type
Definition: id_types.h:67