Line data Source code
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
30 : #include "timpi/parallel_implementation.h"
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
48 : DynaIO::ElementMaps DynaIO::_element_maps = DynaIO::build_element_maps();
49 :
50 :
51 :
52 : // Definition of the static function which constructs the ElementMaps object.
53 16389 : DynaIO::ElementMaps DynaIO::build_element_maps()
54 : {
55 : // Object to be filled up
56 462 : 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 32778 : for (unsigned int i=1; i != 2; ++i)
70 : {
71 : // We only have one element for whom node orders match...
72 16389 : em.add_def(ElementDefinition(EDGE2, i, 1, 1));
73 :
74 16851 : em.add_def(ElementDefinition(EDGE3, i, 1, 2,
75 : {0, 2, 1}));
76 16851 : em.add_def(ElementDefinition(EDGE4, i, 1, 2,
77 : {0, 2, 3, 1}));
78 :
79 16851 : em.add_def(ElementDefinition(QUAD4, i, 2, 1,
80 : {0, 1, 3, 2}));
81 16851 : em.add_def(ElementDefinition(QUAD9, i, 2, 2,
82 : {0, 4, 1, 7, 8, 5, 3, 6, 2}));
83 :
84 16851 : em.add_def(ElementDefinition(HEX8, i, 3, 1,
85 : {0, 1, 3, 2, 4, 5, 7, 6}));
86 32778 : 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 16389 : return em;
93 0 : }
94 :
95 :
96 :
97 16389 : DynaIO::ElementDefinition::ElementDefinition
98 : (ElemType type_in,
99 : dyna_int_type dyna_type_in,
100 : dyna_int_type dim_in,
101 16389 : dyna_int_type p_in) :
102 15465 : type(type_in),
103 15465 : dyna_type(dyna_type_in),
104 15465 : dim(dim_in),
105 16389 : p(p_in)
106 : {
107 16389 : const unsigned int n_nodes = Elem::type_to_n_nodes_map[type_in];
108 16389 : if (n_nodes == invalid_uint)
109 0 : libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
110 16389 : nodes.resize(n_nodes);
111 462 : std::iota(nodes.begin(), nodes.end(), 0);
112 16389 : }
113 :
114 :
115 98334 : DynaIO::ElementDefinition::ElementDefinition
116 : (ElemType type_in,
117 : dyna_int_type dyna_type_in,
118 : dyna_int_type dim_in,
119 : dyna_int_type p_in,
120 98334 : std::vector<unsigned int> && nodes_in) :
121 92790 : type(type_in),
122 92790 : dyna_type(dyna_type_in),
123 92790 : dim(dim_in),
124 92790 : p(p_in),
125 98334 : nodes(nodes_in)
126 98334 : {}
127 :
128 :
129 :
130 805 : DynaIO::DynaIO (MeshBase & mesh,
131 805 : bool keep_spline_nodes) :
132 : MeshInput<MeshBase> (mesh),
133 805 : _keep_spline_nodes (keep_spline_nodes)
134 : {
135 805 : }
136 :
137 :
138 :
139 156 : void DynaIO::read (const std::string & name)
140 : {
141 156 : const bool gzipped_file = Utility::ends_with(name, ".gz");
142 :
143 143 : std::unique_ptr<std::istream> in;
144 :
145 156 : if (gzipped_file)
146 : {
147 : #ifdef LIBMESH_HAVE_GZSTREAM
148 117 : auto inf = std::make_unique<igzstream>();
149 9 : libmesh_assert(inf);
150 9 : inf->open(name.c_str(), std::ios::in);
151 99 : in = std::move(inf);
152 : #else
153 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
154 : #endif
155 90 : }
156 : else
157 : {
158 52 : auto inf = std::make_unique<std::ifstream>();
159 4 : libmesh_assert(inf);
160 :
161 52 : std::string new_name = Utility::unzip_file(name);
162 :
163 48 : inf->open(new_name.c_str(), std::ios::in);
164 44 : in = std::move(inf);
165 40 : }
166 :
167 13 : libmesh_assert(in.get());
168 :
169 156 : this->read_mesh (*in);
170 156 : }
171 :
172 :
173 :
174 156 : void DynaIO::read_mesh(std::istream & in)
175 : {
176 156 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
177 :
178 : // This is a serial-only process for now;
179 : // the Mesh should be read on processor 0 and
180 : // broadcast later
181 13 : libmesh_assert_equal_to (mesh.processor_id(), 0);
182 :
183 169 : libmesh_error_msg_if(!in.good(), "Can't read input stream");
184 :
185 : // clear any data in the mesh
186 156 : mesh.clear();
187 :
188 : // clear any of our own data
189 156 : spline_node_ptrs.clear();
190 13 : 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 13 : 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 26 : block_elem_type,
222 26 : block_n_elem,
223 26 : block_n_nodes, // Number of *spline* nodes constraining elements
224 26 : block_n_coef_vec, // Number of coefficient vectors for each elem
225 26 : block_p,
226 26 : block_dim;
227 : std::vector<dyna_int_type> // indexed from 0 to n_dense_coef_vec_blocks
228 26 : n_coef_vecs_in_subblock, n_coef_comp;
229 13 : unsigned char weight_index = 0;
230 13 : dyna_int_type n_nodes_read = 0,
231 13 : n_elem_blocks_read = 0,
232 13 : n_elems_read = 0,
233 13 : n_elem_nodes_read = 0,
234 13 : n_elem_cvids_read = 0,
235 13 : n_coef_headers_read = 0,
236 13 : n_coef_blocks_read = 0,
237 13 : n_coef_comp_read = 0,
238 13 : n_coef_vecs_read = 0;
239 :
240 : // For reading the file line by line
241 26 : 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 26 : 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 39 : 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 39 : 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 39 : 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 111756 : std::getline(in, s);
267 :
268 121069 : if (in)
269 : {
270 : // Process s...
271 :
272 111756 : if (s.find("B E X T 2.0") == static_cast<std::string::size_type>(0))
273 : {
274 156 : libmesh_error_msg_if(section != FILE_HEADER,
275 : "Found 'B E X T 2.0' outside file header?");
276 :
277 13 : section = PATCH_HEADER;
278 1872 : continue;
279 130 : }
280 :
281 : // Ignore comments
282 111600 : if (s.find("$#") == static_cast<std::string::size_type>(0))
283 1573 : continue;
284 :
285 : // Ignore blank lines
286 109884 : if (s.find_first_not_of(" \t") == std::string::npos)
287 0 : continue;
288 :
289 : // Parse each important section
290 119028 : std::stringstream stream(s);
291 109884 : switch (section) {
292 156 : case PATCH_HEADER:
293 156 : stream >> patch_id;
294 156 : stream >> n_spline_nodes;
295 156 : stream >> n_elem;
296 156 : stream >> n_coef_vec;
297 156 : stream >> weight_control_flag;
298 156 : libmesh_error_msg_if(stream.fail(), "Failure to parse patch header\n");
299 :
300 156 : spline_node_ptrs.resize(n_spline_nodes);
301 156 : 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 156 : const Real default_weight = 1.0;
311 : weight_index = cast_int<unsigned char>
312 299 : (mesh.add_node_datum<Real>("rational_weight", true,
313 : &default_weight));
314 13 : mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
315 26 : mesh.set_default_mapping_data(weight_index);
316 : }
317 :
318 13 : section = NODE_LINES;
319 156 : break;
320 2172 : case NODE_LINES:
321 : {
322 : std::array<dyna_fp_type, 4> xyzw;
323 2172 : stream >> xyzw[0];
324 2172 : stream >> xyzw[1];
325 2172 : stream >> xyzw[2];
326 2172 : stream >> xyzw[3];
327 :
328 26064 : if (weight_control_flag)
329 23868 : 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 28236 : Node *n = spline_node_ptrs[n_nodes_read] =
335 28236 : mesh.add_point(Point(xyzw[0], xyzw[1], xyzw[2]));
336 26064 : Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
337 26064 : elem->set_node(0, n);
338 26064 : elem->subdomain_id() = 1; // Separate id to ease Exodus output
339 26064 : spline_nodeelem_ptrs[n] = elem;
340 : }
341 26064 : ++n_nodes_read;
342 :
343 26064 : libmesh_error_msg_if(stream.fail(), "Failure to parse node line\n");
344 :
345 26064 : if (n_nodes_read >= n_spline_nodes)
346 13 : section = N_ELEM_SUBBLOCKS;
347 2172 : break;
348 156 : case N_ELEM_SUBBLOCKS:
349 156 : stream >> n_elem_blocks;
350 156 : libmesh_error_msg_if(stream.fail(), "Failure to parse n_elem_blocks\n");
351 :
352 156 : block_elem_type.resize(n_elem_blocks);
353 156 : block_n_elem.resize(n_elem_blocks);
354 156 : block_n_nodes.resize(n_elem_blocks);
355 156 : block_n_coef_vec.resize(n_elem_blocks);
356 156 : block_p.resize(n_elem_blocks);
357 156 : block_dim.resize(n_elem_blocks);
358 :
359 156 : elem_global_nodes.resize(n_elem_blocks);
360 156 : elem_constraint_rows.resize(n_elem_blocks);
361 :
362 13 : n_elem_blocks_read = 0;
363 13 : section = ELEM_SUBBLOCK_HEADER;
364 13 : break;
365 156 : case ELEM_SUBBLOCK_HEADER:
366 169 : stream >> block_elem_type[n_elem_blocks_read];
367 169 : stream >> block_n_elem[n_elem_blocks_read];
368 169 : stream >> block_n_nodes[n_elem_blocks_read];
369 169 : stream >> block_n_coef_vec[n_elem_blocks_read];
370 169 : stream >> block_p[n_elem_blocks_read];
371 :
372 156 : libmesh_error_msg_if(stream.fail(), "Failure to parse elem block\n");
373 :
374 156 : 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 156 : stream >> block_other_p;
378 156 : if (!stream.fail())
379 : {
380 156 : block_dim[n_elem_blocks_read] = 2; // Found a second dimension!
381 :
382 169 : if (block_other_p != block_p[n_elem_blocks_read])
383 0 : libmesh_not_implemented(); // We don't support p anisotropy
384 :
385 156 : stream >> block_other_p;
386 156 : if (!stream.fail())
387 : {
388 48 : block_dim[n_elem_blocks_read] = 3; // Found a third dimension!
389 :
390 52 : if (block_other_p != block_p[n_elem_blocks_read])
391 0 : libmesh_not_implemented();
392 : }
393 : }
394 :
395 : {
396 26 : auto & block_global_nodes = elem_global_nodes[n_elem_blocks_read];
397 26 : auto & block_constraint_rows = elem_constraint_rows[n_elem_blocks_read];
398 :
399 169 : block_global_nodes.resize(block_n_elem[n_elem_blocks_read]);
400 169 : block_constraint_rows.resize(block_n_elem[n_elem_blocks_read]);
401 :
402 16669 : for (auto e : make_range(block_n_elem[n_elem_blocks_read]))
403 : {
404 19250 : block_global_nodes[e].resize(block_n_nodes[n_elem_blocks_read]);
405 19250 : block_constraint_rows[e].resize(block_n_coef_vec[n_elem_blocks_read]);
406 0 : }
407 : }
408 :
409 156 : n_elem_blocks_read++;
410 156 : if (n_elem_blocks_read == n_elem_blocks)
411 : {
412 13 : n_elem_blocks_read = 0;
413 13 : n_elems_read = 0;
414 13 : section = ELEM_NODES_LINES;
415 : }
416 13 : break;
417 18036 : case ELEM_NODES_LINES:
418 : {
419 : const int end_node_to_read =
420 18164 : std::min(block_n_nodes[n_elem_blocks_read], n_elem_nodes_read + max_ints_per_line);
421 180360 : for (; n_elem_nodes_read != end_node_to_read; ++n_elem_nodes_read)
422 : {
423 : dyna_int_type node_id;
424 162324 : stream >> node_id;
425 162324 : node_id--;
426 202905 : 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 162324 : libmesh_error_msg_if(stream.fail(), "Failure to parse elem nodes\n");
432 : }
433 :
434 19539 : if (n_elem_nodes_read == block_n_nodes[n_elem_blocks_read])
435 : {
436 1375 : n_elem_nodes_read = 0;
437 1375 : section = ELEM_COEF_VEC_IDS;
438 : }
439 : }
440 1503 : break;
441 18036 : case ELEM_COEF_VEC_IDS:
442 : {
443 : const int end_cvid_to_read =
444 18164 : std::min(block_n_coef_vec[n_elem_blocks_read], n_elem_cvids_read + max_ints_per_line);
445 180360 : for (; n_elem_cvids_read != end_cvid_to_read; ++n_elem_cvids_read)
446 : {
447 : dyna_int_type node_cvid;
448 162324 : stream >> node_cvid;
449 162324 : node_cvid--;
450 :
451 202905 : 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 162324 : libmesh_error_msg_if(stream.fail(), "Failure to parse elem cvids\n");
457 : }
458 19539 : if (n_elem_cvids_read == block_n_nodes[n_elem_blocks_read])
459 : {
460 1375 : n_elem_cvids_read = 0;
461 16500 : n_elems_read++;
462 1375 : section = ELEM_NODES_LINES; // Read another elem, nodes first
463 17875 : if (n_elems_read == block_n_elem[n_elem_blocks_read])
464 : {
465 13 : n_elems_read = 0;
466 156 : n_elem_blocks_read++;
467 156 : if (n_elem_blocks_read == n_elem_blocks)
468 13 : section = N_COEF_BLOCKS; // Move on to coefficient vectors
469 : }
470 : }
471 : }
472 1503 : break;
473 156 : case N_COEF_BLOCKS:
474 : {
475 156 : stream >> n_dense_coef_vec_blocks;
476 : dyna_int_type n_sparse_coef_vec_blocks;
477 156 : stream >> n_sparse_coef_vec_blocks;
478 :
479 156 : libmesh_error_msg_if(stream.fail(), "Failure to parse n_*_coef_vec_blocks\n");
480 :
481 156 : if (n_sparse_coef_vec_blocks != 0)
482 0 : libmesh_not_implemented();
483 :
484 156 : dense_constraint_vecs.resize(n_dense_coef_vec_blocks);
485 156 : n_coef_vecs_in_subblock.resize(n_dense_coef_vec_blocks);
486 156 : n_coef_comp.resize(n_dense_coef_vec_blocks);
487 :
488 0 : section = N_VECS_PER_BLOCK;
489 : }
490 156 : break;
491 156 : case N_VECS_PER_BLOCK:
492 169 : stream >> n_coef_vecs_in_subblock[n_coef_headers_read];
493 169 : stream >> n_coef_comp[n_coef_headers_read];
494 :
495 156 : libmesh_error_msg_if(stream.fail(), "Failure to parse dense coef subblock header\n");
496 :
497 26 : dense_constraint_vecs[n_coef_headers_read].resize
498 169 : (n_coef_vecs_in_subblock[n_coef_headers_read]);
499 :
500 10536 : for (auto & vec : dense_constraint_vecs[n_coef_headers_read])
501 11245 : vec.resize(n_coef_comp[n_coef_headers_read]);
502 :
503 156 : n_coef_headers_read++;
504 156 : if (n_coef_headers_read == n_dense_coef_vec_blocks)
505 : {
506 13 : n_coef_headers_read = 0;
507 13 : section = COEF_VEC_COMPONENTS;
508 : }
509 13 : break;
510 46968 : case COEF_VEC_COMPONENTS:
511 : {
512 : auto & current_vec =
513 50882 : dense_constraint_vecs[n_coef_blocks_read][n_coef_vecs_read];
514 :
515 : const int end_coef_to_read =
516 7828 : std::min(n_coef_comp[n_coef_blocks_read],
517 50017 : n_coef_comp_read + max_fps_per_line);
518 258324 : for (; n_coef_comp_read != end_coef_to_read; ++n_coef_comp_read)
519 : {
520 : dyna_fp_type coef_comp;
521 17613 : stream >> coef_comp;
522 :
523 228969 : 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 211356 : libmesh_error_msg_if(stream.fail(), "Failure to parse coefficients\n");
529 : }
530 50882 : if (n_coef_comp_read == n_coef_comp[n_coef_blocks_read])
531 : {
532 865 : n_coef_comp_read = 0;
533 10380 : n_coef_vecs_read++;
534 11245 : if (n_coef_vecs_read == n_coef_vecs_in_subblock[n_coef_blocks_read])
535 : {
536 13 : n_coef_vecs_read = 0;
537 156 : n_coef_blocks_read++;
538 156 : if (n_coef_blocks_read == n_dense_coef_vec_blocks)
539 13 : section = END_OF_FILE;
540 : }
541 : }
542 : }
543 3914 : break;
544 0 : default:
545 0 : libmesh_error();
546 : }
547 :
548 109884 : if (section == END_OF_FILE)
549 13 : break;
550 91570 : } // if (in)
551 0 : else if (in.eof())
552 0 : libmesh_error_msg("Premature end of file");
553 : else
554 0 : libmesh_error_msg("Input stream failure! Perhaps the file does not exist?");
555 9300 : }
556 :
557 : // Merge dense_constraint_vecs blocks
558 156 : if (n_dense_coef_vec_blocks)
559 0 : for (auto coef_vec_block :
560 156 : IntRange<dyna_int_type>(1, n_dense_coef_vec_blocks))
561 : {
562 0 : auto & dcv0 = dense_constraint_vecs[0];
563 0 : auto & dcvi = dense_constraint_vecs[coef_vec_block];
564 0 : dcv0.insert(dcv0.end(),
565 : std::make_move_iterator(dcvi.begin()),
566 0 : std::make_move_iterator(dcvi.end()));
567 : }
568 156 : 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 312 : 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 26 : std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
579 :
580 13 : auto & constraint_rows = mesh.get_constraint_rows();
581 :
582 312 : for (auto block_num : make_range(n_elem_blocks))
583 : {
584 182 : elem_constraint_mat[block_num].resize(block_n_elem[block_num]);
585 :
586 1505 : for (auto elem_num :
587 18044 : 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 16500 : find_elem_definition(block_elem_type[block_num],
592 2750 : block_dim[block_num],
593 6875 : block_p[block_num]);
594 :
595 16500 : auto elem = Elem::build(elem_defn.type);
596 16500 : libmesh_error_msg_if(elem->dim() != block_dim[block_num],
597 : "Elem dim " << elem->dim() <<
598 : " != block_dim " << block_dim[block_num]);
599 :
600 17875 : auto & my_constraint_rows = elem_constraint_rows[block_num][elem_num];
601 4125 : auto & my_global_nodes = elem_global_nodes[block_num][elem_num];
602 4125 : auto & my_constraint_mat = elem_constraint_mat[block_num][elem_num];
603 :
604 17875 : my_constraint_mat.resize(block_n_coef_vec[block_num]);
605 27277 : for (auto spline_node_index :
606 193726 : make_range(block_n_coef_vec[block_num]))
607 175851 : my_constraint_mat[spline_node_index].resize(elem->n_nodes());
608 :
609 27277 : for (auto spline_node_index :
610 193726 : 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 175851 : my_constraint_rows[spline_node_index];
615 :
616 13527 : dyna_int_type coef_block_num = 0;
617 13527 : dyna_int_type first_block_coef_vec = 0;
618 324648 : for (; elem_coef_vec_index >= first_block_coef_vec &&
619 162324 : coef_block_num != n_dense_coef_vec_blocks; ++coef_block_num)
620 : {
621 175851 : 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 162324 : 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 162324 : coef_block_num--;
630 :
631 162324 : 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 152847 : for (auto elem_node_index :
638 1996488 : make_range(elem->n_nodes()))
639 2139858 : my_constraint_mat[spline_node_index][elem_node_index] =
640 1987011 : dense_constraint_vecs[0][elem_coef_vec_index][elem_node_index];
641 : }
642 :
643 13527 : for (auto elem_node_index :
644 192351 : make_range(elem->n_nodes()))
645 : {
646 162324 : 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 27054 : std::vector<std::pair<dof_id_type, Real>> key;
652 :
653 288117 : for (auto spline_node_index :
654 2162862 : make_range(block_n_coef_vec[block_num]))
655 : {
656 : const dyna_int_type elem_coef_vec_index =
657 1834164 : my_constraint_rows[spline_node_index];
658 :
659 : const Real coef =
660 1834164 : libmesh_vector_at(dense_constraint_vecs[0],
661 1834164 : elem_coef_vec_index)[elem_node_index];
662 :
663 : // Global nodes are supposed to be in sorted order
664 1834164 : if (global_node_idx != DofObject::invalid_id)
665 1811160 : libmesh_error_msg_if(my_global_nodes[spline_node_index] <= global_node_idx,
666 : "Found unsorted global node");
667 :
668 1834164 : global_node_idx = my_global_nodes[spline_node_index];
669 :
670 1834164 : if (coef != 0) // Ignore irrelevant spline nodes
671 417852 : key.emplace_back(global_node_idx, coef);
672 : }
673 :
674 162324 : if (const auto local_node_it = local_nodes.find(key);
675 13527 : local_node_it != local_nodes.end())
676 89934 : elem->set_node(elem_defn.nodes[elem_node_index],
677 13836 : local_node_it->second);
678 : else
679 : {
680 6609 : Point p(0);
681 6609 : Real w = 0;
682 6609 : std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
683 :
684 143571 : for (auto spline_node_index :
685 1093170 : make_range(block_n_coef_vec[block_num]))
686 : {
687 : const dof_id_type my_node_idx =
688 929772 : my_global_nodes[spline_node_index];
689 :
690 : const dyna_int_type elem_coef_vec_index =
691 929772 : my_constraint_rows[spline_node_index];
692 :
693 : Node * spline_node =
694 929772 : libmesh_vector_at(spline_node_ptrs,
695 : my_node_idx);
696 :
697 : const Real coef =
698 929772 : libmesh_vector_at(dense_constraint_vecs[0],
699 929772 : elem_coef_vec_index)[elem_node_index];
700 929772 : p.add_scaled(*spline_node, coef);
701 929772 : 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 929772 : if (coef)
707 : {
708 : Elem * nodeelem =
709 162240 : libmesh_map_find(spline_nodeelem_ptrs, spline_node);
710 175760 : constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
711 : }
712 : }
713 :
714 79308 : Node *n = mesh.add_point(p);
715 79308 : if (weight_control_flag)
716 70308 : n->set_extra_datum<Real>(weight_index, w);
717 79308 : local_nodes[key] = n;
718 85917 : elem->set_node(elem_defn.nodes[elem_node_index], n);
719 :
720 79308 : constraint_rows[n] = constraint_row;
721 : }
722 : }
723 :
724 19250 : mesh.add_elem(std::move(elem));
725 13750 : }
726 : }
727 :
728 156 : if (!_keep_spline_nodes)
729 36 : MeshTools::clear_spline_nodes(MeshInput<MeshBase>::mesh());
730 286 : }
731 :
732 :
733 0 : void DynaIO::add_spline_constraints(DofMap &,
734 : unsigned int,
735 : unsigned int)
736 : {
737 0 : }
738 :
739 :
740 : #ifdef LIBMESH_ENABLE_DEPRECATED
741 0 : void DynaIO::clear_spline_nodes()
742 : {
743 : libmesh_deprecated();
744 :
745 0 : MeshTools::clear_spline_nodes(MeshInput<MeshBase>::mesh());
746 0 : }
747 : #endif // LIBMESH_ENABLE_DEPRECATED
748 :
749 :
750 :
751 : const DynaIO::ElementDefinition &
752 16500 : DynaIO::find_elem_definition(dyna_int_type dyna_elem,
753 : int dim, int p)
754 : {
755 : auto eletypes_it =
756 16500 : _element_maps.in.find(std::make_tuple(dyna_elem, dim, p));
757 :
758 : // Make sure we actually found something
759 16500 : 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 16500 : return eletypes_it->second;
766 : }
767 :
768 :
769 :
770 : const DynaIO::ElementDefinition &
771 65 : DynaIO::find_elem_definition(ElemType libmesh_elem,
772 : int libmesh_dbg_var(dim),
773 : int libmesh_dbg_var(p))
774 : {
775 : auto eletypes_it =
776 0 : _element_maps.out.find(libmesh_elem);
777 :
778 : // Make sure we actually found something
779 65 : 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 0 : libmesh_assert_equal_to(eletypes_it->second.dim, dim);
786 0 : libmesh_assert_equal_to(eletypes_it->second.p, p);
787 :
788 :
789 65 : return eletypes_it->second;
790 : }
791 :
792 :
793 :
794 :
795 : } // namespace libMesh
|
