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/abaqus_io.h"
20 : #include "libmesh/point.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/enum_to_string.h"
23 : #include "libmesh/boundary_info.h"
24 : #include "libmesh/utility.h"
25 :
26 : // gzstream for reading compressed files as a stream
27 : #ifdef LIBMESH_HAVE_GZSTREAM
28 : # include "gzstream.h"
29 : #endif
30 :
31 : // C++ includes
32 : #include <unordered_map>
33 : #include <string>
34 : #include <cstdlib> // std::strtol
35 : #include <sstream>
36 : #include <cctype> // isspace
37 :
38 : // Anonymous namespace to hold mapping Data for Abaqus/libMesh element types
39 : namespace
40 : {
41 : using namespace libMesh;
42 :
43 : /**
44 : * Data structure used for mapping Abaqus IDs to libMesh IDs, and
45 : * eventually (possibly) vice-versa.
46 : */
47 : struct ElementDefinition
48 : {
49 : // Maps (zero-based!) Abaqus local node numbers to libmesh local node numbers
50 : std::vector<unsigned> abaqus_zero_based_node_id_to_libmesh_node_id;
51 :
52 : // Maps (zero-based!) Abaqus side numbers to libmesh side numbers
53 : std::vector<unsigned short> abaqus_zero_based_side_id_to_libmesh_side_id;
54 : };
55 :
56 : /**
57 : * Locally-available map containing all element data.
58 : */
59 : std::map<ElemType, ElementDefinition> eletypes;
60 :
61 : /**
62 : * Helper function to fill up eletypes map
63 : */
64 168 : void add_eletype_entry(ElemType libmesh_elem_type,
65 : const unsigned * node_map,
66 : unsigned node_map_size,
67 : const unsigned short * side_map,
68 : unsigned side_map_size)
69 : {
70 : // If map entry does not exist, this will create it
71 168 : ElementDefinition & map_entry = eletypes[libmesh_elem_type];
72 :
73 :
74 : // Use the "swap trick" from Scott Meyer's "Effective STL" to swap
75 : // an unnamed temporary vector into the map_entry's vector. Note:
76 : // the vector(iter, iter) constructor is used.
77 336 : std::vector<unsigned>
78 168 : (node_map, node_map+node_map_size).swap
79 14 : (map_entry.abaqus_zero_based_node_id_to_libmesh_node_id);
80 :
81 336 : std::vector<unsigned short>
82 168 : (side_map, side_map+side_map_size).swap
83 14 : (map_entry.abaqus_zero_based_side_id_to_libmesh_side_id);
84 168 : }
85 :
86 :
87 : /**
88 : * Helper function to initialize the eletypes map.
89 : */
90 48 : void init_eletypes ()
91 : {
92 : // This should happen only once. The first time this method is
93 : // called the eletypes data structure will be empty, and we will
94 : // fill it. Any subsequent calls will find an initialized
95 : // eletypes map and will do nothing.
96 48 : if (eletypes.empty())
97 : {
98 : {
99 : // NODEELEM
100 12 : const unsigned int node_map[] = {0}; // identity
101 12 : const unsigned short side_map[] = {0}; // identity
102 12 : add_eletype_entry(NODEELEM, node_map, 1, side_map, 1);
103 : }
104 :
105 : {
106 : // EDGE2
107 12 : const unsigned int node_map[] = {0,1}; // identity
108 12 : const unsigned short side_map[] = {0,1}; // identity
109 12 : add_eletype_entry(EDGE2, node_map, 2, side_map, 2);
110 : }
111 :
112 : {
113 : // TRI3
114 12 : const unsigned int node_map[] = {0,1,2}; // identity
115 12 : const unsigned short side_map[] = {0,1,2}; // identity
116 12 : add_eletype_entry(TRI3, node_map, 3, side_map, 3);
117 : }
118 :
119 : {
120 : // TRI6
121 12 : const unsigned int node_map[] = {0,1,2,3,4,5}; // identity
122 12 : const unsigned short side_map[] = {0,1,2}; // identity
123 12 : add_eletype_entry(TRI6, node_map, 6, side_map, 3);
124 : }
125 :
126 : {
127 : // QUAD4
128 12 : const unsigned int node_map[] = {0,1,2,3}; // identity
129 12 : const unsigned short side_map[] = {0,1,2,3}; // identity
130 12 : add_eletype_entry(QUAD4, node_map, 4, side_map, 4);
131 : }
132 :
133 : {
134 : // QUAD8
135 12 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7}; // identity
136 12 : const unsigned short side_map[] = {0,1,2,3}; // identity
137 12 : add_eletype_entry(QUAD8, node_map, 8, side_map, 4);
138 : }
139 :
140 : {
141 : // TET4
142 12 : const unsigned int node_map[] = {0,1,2,3}; // identity
143 12 : const unsigned short side_map[] = {0,1,2,3}; // identity
144 12 : add_eletype_entry(TET4, node_map, 4, side_map, 4);
145 : }
146 :
147 : {
148 : // TET10
149 12 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,9}; // identity
150 12 : const unsigned short side_map[] = {0,1,2,3}; // identity
151 12 : add_eletype_entry(TET10, node_map, 10, side_map, 4);
152 : }
153 :
154 : {
155 : // HEX8
156 12 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7}; // identity
157 12 : const unsigned short side_map[] = {0,5,1,2,3,4}; // inverse = 0,2,3,4,5,1
158 12 : add_eletype_entry(HEX8, node_map, 8, side_map, 6);
159 : }
160 :
161 : {
162 : // HEX20
163 12 : const unsigned int node_map[] = // map is its own inverse
164 : {0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15};
165 12 : const unsigned short side_map[] = // inverse = 0,2,3,4,5,1
166 : {0,5,1,2,3,4};
167 12 : add_eletype_entry(HEX20, node_map, 20, side_map, 6);
168 : }
169 :
170 : {
171 : // HEX27
172 12 : const unsigned int node_map[] = // inverse = ...,21,23,24,25,26,22,20
173 : {0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15,26,20,25,21,22,23,24};
174 12 : const unsigned short side_map[] = // inverse = 0,2,3,4,5,1
175 : {0,5,1,2,3,4};
176 12 : add_eletype_entry(HEX27, node_map, 27, side_map, 6);
177 : }
178 :
179 : {
180 : // PRISM6
181 12 : const unsigned int node_map[] = {0,1,2,3,4,5}; // identity
182 12 : const unsigned short side_map[] = {0,4,1,2,3}; // inverse = 0,2,3,4,1
183 12 : add_eletype_entry(PRISM6, node_map, 6, side_map, 5);
184 : }
185 :
186 : {
187 : // PRISM15
188 12 : const unsigned int node_map[] = // map is its own inverse
189 : {0,1,2,3,4,5,6,7,8,12,13,14,9,10,11};
190 12 : const unsigned short side_map[] = // inverse = 0,2,3,4,1
191 : {0,4,1,2,3};
192 12 : add_eletype_entry(PRISM15, node_map, 15, side_map, 5);
193 : }
194 :
195 : {
196 : // PRISM18
197 12 : const unsigned int node_map[] = // map is its own inverse
198 : {0,1,2,3,4,5,6,7,8,12,13,14,9,10,11,15,16,17};
199 12 : const unsigned short side_map[] = // inverse = 0,2,3,4,1
200 : {0,4,1,2,3};
201 12 : add_eletype_entry(PRISM18, node_map, 18, side_map, 5);
202 : }
203 : } // if (eletypes.empty())
204 48 : }
205 : } // anonymous namespace
206 :
207 :
208 :
209 :
210 :
211 : namespace libMesh
212 : {
213 :
214 142 : AbaqusIO::AbaqusIO (MeshBase & mesh_in) :
215 : MeshInput<MeshBase> (mesh_in),
216 134 : build_sidesets_from_nodesets(false),
217 142 : _already_seen_part(false)
218 : {
219 142 : }
220 :
221 :
222 :
223 134 : AbaqusIO::~AbaqusIO () = default;
224 :
225 :
226 :
227 24 : void AbaqusIO::read (const std::string & fname)
228 : {
229 : // Get a reference to the mesh we are reading
230 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
231 :
232 : // Clear any existing mesh data
233 24 : the_mesh.clear();
234 :
235 : // Open stream for reading
236 24 : const bool gzipped_file = Utility::ends_with(fname, ".gz");
237 :
238 24 : if (gzipped_file)
239 : {
240 : #ifdef LIBMESH_HAVE_GZSTREAM
241 26 : auto inf = std::make_unique<igzstream>();
242 2 : libmesh_assert(inf);
243 2 : inf->open(fname.c_str(), std::ios::in);
244 22 : _in = std::move(inf); // class takes ownership as base class pointer
245 : #else
246 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
247 : #endif
248 20 : }
249 : else
250 : {
251 0 : auto inf = std::make_unique<std::ifstream>();
252 0 : libmesh_assert(inf);
253 :
254 0 : std::string new_name = Utility::unzip_file(fname);
255 0 : inf->open(new_name.c_str(), std::ios::in);
256 0 : libmesh_assert(inf->good());
257 0 : _in = std::move(inf); // class takes ownership as base class pointer
258 0 : }
259 :
260 : // Initialize the elems_of_dimension array. We will use this in a
261 : // "1-based" manner so that elems_of_dimension[d]==true means
262 : // elements of dimension d have been seen.
263 24 : elems_of_dimension.resize(4, false);
264 :
265 : // Read file line-by-line... this is based on a set of different
266 : // test input files. I have not looked at the full input file
267 : // specs for Abaqus.
268 4 : std::string s;
269 : while (true)
270 : {
271 : // Try to read something. This may set EOF!
272 1332 : std::getline(*_in, s);
273 :
274 1443 : if (*_in)
275 : {
276 : // Process s...
277 : //
278 : // There are many sections in Abaqus files, we read some
279 : // but others are just ignored... Some sections may occur
280 : // more than once. For example for a hybrid grid, you
281 : // will have multiple *Element sections...
282 :
283 : // Some Abaqus files use all upper-case for section names,
284 : // so we will just convert s to uppercase
285 1308 : std::string upper(s);
286 109 : std::transform(upper.begin(), upper.end(), upper.begin(), ::toupper);
287 :
288 : // 0.) Look for the "*Part" section
289 1308 : if (upper.find("*PART") == static_cast<std::string::size_type>(0))
290 : {
291 0 : libmesh_error_msg_if
292 : (_already_seen_part,
293 : "We currently don't support reading Abaqus files with multiple PART sections");
294 :
295 0 : _already_seen_part = true;
296 : }
297 :
298 : // 1.) Look for the "*Nodes" section
299 1308 : if (upper.find("*NODE") == static_cast<std::string::size_type>(0))
300 : {
301 : // Some sections that begin with *NODE are actually
302 : // unrelated sections which we want to skip. I
303 : // have only seen this with a single space, but it would
304 : // probably be more robust to remove whitespace before
305 : // making this check.
306 4 : bool skip_this_section = false;
307 238 : std::vector<std::string> keywords_to_ignore = {"OUTPUT", "PRINT", "FILE", "RESPONSE"};
308 240 : for (auto & affix : keywords_to_ignore) {
309 224 : std::stringstream keyword;
310 176 : keyword << "*NODE " << affix;
311 192 : if (upper.find(keyword.str()) == static_cast<std::string::size_type>(0))
312 2 : skip_this_section = true;
313 160 : }
314 48 : if (skip_this_section)
315 4 : continue;
316 :
317 : // Some *Node sections also specify an Nset name on the same line.
318 : // Look for one here.
319 50 : std::string nset_name = this->parse_label(s, "nset");
320 :
321 : // Process any lines of comments that may be present
322 24 : this->process_and_discard_comments();
323 :
324 : // Read a block of nodes
325 46 : this->read_nodes(nset_name);
326 40 : }
327 :
328 :
329 :
330 : // 2.) Look for the "*Element" section
331 1260 : else if (upper.find("*ELEMENT,") == static_cast<std::string::size_type>(0))
332 : {
333 : // Some sections that begin with *ELEMENT are actually
334 : // "*ELEMENT OUTPUT" sections which we want to skip. I
335 : // have only seen this with a single space, but it would
336 : // probably be more robust to remove whitespace before
337 : // making this check.
338 2 : bool skip_this_section = false;
339 120 : std::vector<std::string> keywords_to_ignore = {"OUTPUT", "MATRIX", "PROPERTIES", "RESPONSE"};
340 120 : for (auto & affix : keywords_to_ignore) {
341 112 : std::stringstream keyword;
342 88 : keyword << "*ELEMENT " << affix;
343 96 : if (upper.find(keyword.str()) == static_cast<std::string::size_type>(0))
344 0 : skip_this_section = true;
345 80 : }
346 24 : if (skip_this_section)
347 0 : continue;
348 :
349 : // Some *Element sections also specify an Elset name on the same line.
350 : // Look for one here.
351 50 : std::string elset_name = this->parse_label(s, "elset");
352 :
353 : // Process any lines of comments that may be present
354 24 : this->process_and_discard_comments();
355 :
356 : // Read a block of elements
357 46 : this->read_elements(upper, elset_name);
358 20 : }
359 :
360 :
361 :
362 : // 3.) Look for a Nodeset section
363 1236 : else if (upper.find("*NSET") == static_cast<std::string::size_type>(0))
364 : {
365 182 : std::string nset_name = this->parse_label(s, "nset");
366 :
367 : // I haven't seen an unnamed nset yet, but let's detect it
368 : // just in case...
369 84 : libmesh_error_msg_if(nset_name == "", "Unnamed nset encountered!");
370 :
371 : // Is this a "generated" nset, i.e. one which has three
372 : // entries corresponding to (first, last, stride)?
373 84 : bool is_generated = this->detect_generated_set(upper);
374 :
375 : // Process any lines of comments that may be present
376 84 : this->process_and_discard_comments();
377 :
378 : // Read the IDs, storing them in _nodeset_ids
379 84 : if (is_generated)
380 92 : this->generate_ids(nset_name, _nodeset_ids);
381 : else
382 69 : this->read_ids(nset_name, _nodeset_ids);
383 : } // *Nodeset
384 :
385 :
386 :
387 : // 4.) Look for an Elset section
388 1152 : else if (upper.find("*ELSET") == static_cast<std::string::size_type>(0))
389 : {
390 468 : std::string elset_name = this->parse_label(s, "elset");
391 :
392 : // I haven't seen an unnamed elset yet, but let's detect it
393 : // just in case...
394 216 : libmesh_error_msg_if(elset_name == "", "Unnamed elset encountered!");
395 :
396 : // Is this a "generated" elset, i.e. one which has three
397 : // entries corresponding to (first, last, stride)?
398 216 : bool is_generated = this->detect_generated_set(upper);
399 :
400 : // Process any lines of comments that may be present
401 216 : this->process_and_discard_comments();
402 :
403 : // Read the IDs, storing them in _elemset_ids
404 216 : if (is_generated)
405 230 : this->generate_ids(elset_name, _elemset_ids);
406 : else
407 184 : this->read_ids(elset_name, _elemset_ids);
408 : } // *Elset
409 :
410 :
411 :
412 : // 5.) Look for a Surface section. Need to be a little
413 : // careful, since there are also "surface interaction"
414 : // sections we don't want to read here.
415 936 : else if (upper.find("*SURFACE,") == static_cast<std::string::size_type>(0))
416 : {
417 : // Get the name from the Name=Foo label. This will be the map key.
418 0 : std::string sideset_name = this->parse_label(s, "name");
419 :
420 : // Some (all?) surfaces also declare a type, which can
421 : // help us figure out how they should be parsed, so we
422 : // pass this to the read_sideset() function.
423 0 : std::string sideset_type = this->parse_label(s, "type");
424 :
425 : // Process any lines of comments that may be present
426 0 : this->process_and_discard_comments();
427 :
428 : // Read the sideset IDs
429 0 : this->read_sideset(sideset_name, sideset_type, _sideset_ids);
430 : }
431 :
432 : // Derived classes could override this to add support for additional sections
433 1284 : extended_parsing(upper);
434 :
435 1177 : continue;
436 : } // if (*_in)
437 :
438 : // If !file, check to see if EOF was set. If so, break out
439 : // of while loop.
440 24 : if (_in->eof())
441 2 : break;
442 :
443 : // If !in and !in.eof(), stream is in a bad state!
444 0 : libmesh_error_msg("Stream is bad! Perhaps the file: " << fname << " does not exist?");
445 109 : } // while
446 :
447 : // Set the Mesh dimension based on the highest dimension element seen.
448 24 : the_mesh.set_mesh_dimension(this->max_elem_dimension_seen());
449 :
450 : // Set element IDs based on the element sets.
451 24 : this->assign_subdomain_ids();
452 :
453 : // Assign nodeset values to the BoundaryInfo object
454 24 : this->assign_boundary_node_ids();
455 :
456 : // Assign sideset values in the BoundaryInfo object
457 24 : this->assign_sideset_ids();
458 :
459 : // If the Abaqus file contains only nodesets, we can have libmesh
460 : // generate sidesets from them. This BoundaryInfo function currently
461 : // *overwrites* existing sidesets in surprising ways, so we don't
462 : // call it if there are already sidesets present in the original file.
463 24 : if (build_sidesets_from_nodesets && the_mesh.get_boundary_info().n_boundary_conds() == 0)
464 0 : the_mesh.get_boundary_info().build_side_list_from_node_list();
465 :
466 : // Delete lower-dimensional elements from the Mesh. We assume these
467 : // were only used for setting BCs, and aren't part of the actual
468 : // Mesh.
469 : {
470 24 : unsigned char max_dim = this->max_elem_dimension_seen();
471 :
472 44222 : for (auto & elem : the_mesh.element_ptr_range())
473 24096 : if (elem->dim() < max_dim)
474 20 : the_mesh.delete_elem(elem);
475 : }
476 144 : }
477 :
478 :
479 :
480 :
481 :
482 :
483 :
484 24 : void AbaqusIO::read_nodes(std::string nset_name)
485 : {
486 : // Get a reference to the mesh we are reading
487 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
488 :
489 : // In the input files I have, Abaqus neither tells what
490 : // the mesh dimension is nor how many nodes it has...
491 : //
492 : // The node line format is:
493 : // id, x, y, z
494 : // and you do have to parse out the commas.
495 : // The z-coordinate will only be present for 3D meshes
496 :
497 : // Temporary variables for parsing lines of text
498 : char c;
499 4 : std::string line;
500 :
501 : // We need to duplicate some of the read_ids code if this *NODE
502 : // section also defines an NSET. We'll set up the id_storage
503 : // pointer and push back IDs into this vector in the loop below...
504 2 : std::vector<dof_id_type> * id_storage = nullptr;
505 24 : if (nset_name != "")
506 0 : id_storage = &(_nodeset_ids[nset_name]);
507 :
508 : // We will read nodes until the next line begins with *, since that will be the
509 : // next section.
510 : // TODO: Is Abaqus guaranteed to start the line with '*' or can there be leading white space?
511 33516 : while (_in->peek() != '*' && _in->peek() != EOF)
512 : {
513 : // Read an entire line which corresponds to a single point's id
514 : // and (x,y,z) values.
515 33492 : std::getline(*_in, line);
516 :
517 : // Remove all whitespace characters from the line. This way we
518 : // can do the remaining parsing without worrying about tabs,
519 : // different numbers of spaces, etc.
520 33492 : strip_ws(line);
521 :
522 : // Make a stream out of the modified line so we can stream values
523 : // from it in the usual way.
524 33492 : std::stringstream ss(line);
525 :
526 : // Values to be read in from file
527 33492 : dof_id_type abaqus_node_id=0;
528 33492 : Real x=0, y=0, z=0;
529 :
530 : // Note: we assume *at least* 2D points here, should we worry about
531 : // trying to read 1D Abaqus meshes?
532 64193 : ss >> abaqus_node_id >> c >> x >> c >> y;
533 :
534 : // Peek at the next character. If it is a comma, then there is another
535 : // value to read!
536 33492 : if (ss.peek() == ',')
537 33492 : ss >> c >> z;
538 :
539 : // If this *NODE section defines an NSET, also store the abaqus ID in id_storage
540 33492 : if (id_storage)
541 0 : id_storage->push_back(abaqus_node_id);
542 :
543 : // Convert from Abaqus 1-based to libMesh 0-based numbering
544 33492 : libmesh_error_msg_if(abaqus_node_id < 1,
545 : "Invalid Abaqus node ID found");
546 33492 : const dof_id_type libmesh_node_id = abaqus_node_id-1;
547 :
548 33492 : libmesh_error_msg_if(the_mesh.query_node_ptr(libmesh_node_id),
549 : "Duplicate Abaqus node ID found");
550 :
551 : // Add the point to the mesh using libmesh's numbering,
552 : // and post-increment the libmesh node counter.
553 36283 : the_mesh.add_point(Point(x,y,z), libmesh_node_id);
554 27910 : } // while
555 24 : }
556 :
557 :
558 :
559 :
560 :
561 24 : void AbaqusIO::read_elements(std::string upper, std::string elset_name)
562 : {
563 : // Get a reference to the mesh we are reading
564 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
565 :
566 : // initialize the eletypes map (eletypes is a file-global variable)
567 24 : init_eletypes();
568 :
569 24 : ElemType elem_type = INVALID_ELEM;
570 2 : unsigned n_nodes_per_elem = 0;
571 :
572 : // Within s, we should have "type=XXXX"
573 48 : if (upper.find("T3D2") != std::string::npos ||
574 22 : upper.find("B31") != std::string::npos)
575 : {
576 0 : elem_type = EDGE2;
577 0 : n_nodes_per_elem = 2;
578 0 : elems_of_dimension[1] = true;
579 : }
580 24 : else if (upper.find("B32") != std::string::npos)
581 : {
582 0 : elem_type = EDGE3;
583 0 : n_nodes_per_elem = 3;
584 0 : elems_of_dimension[1] = true;
585 : }
586 46 : else if (upper.find("S3") != std::string::npos ||
587 48 : upper.find("CPE3") != std::string::npos ||
588 22 : upper.find("2D3") != std::string::npos)
589 : {
590 0 : elem_type = TRI3;
591 0 : n_nodes_per_elem = 3;
592 0 : elems_of_dimension[2] = true;
593 : }
594 46 : else if (upper.find("CPE4") != std::string::npos ||
595 46 : upper.find("S4") != std::string::npos ||
596 70 : upper.find("CPEG4") != std::string::npos ||
597 22 : upper.find("2D4") != std::string::npos)
598 : {
599 0 : elem_type = QUAD4;
600 0 : n_nodes_per_elem = 4;
601 0 : elems_of_dimension[2] = true;
602 : }
603 46 : else if (upper.find("CPE6") != std::string::npos ||
604 46 : upper.find("S6") != std::string::npos ||
605 70 : upper.find("CPEG6") != std::string::npos ||
606 22 : upper.find("2D6") != std::string::npos)
607 : {
608 0 : elem_type = TRI6;
609 0 : n_nodes_per_elem = 6;
610 0 : elems_of_dimension[2] = true;
611 : }
612 46 : else if (upper.find("CPE8") != std::string::npos ||
613 46 : upper.find("S8") != std::string::npos ||
614 70 : upper.find("CPEG8") != std::string::npos ||
615 22 : upper.find("2D8") != std::string::npos)
616 : {
617 0 : elem_type = QUAD8;
618 0 : n_nodes_per_elem = 8;
619 0 : elems_of_dimension[2] = true;
620 : }
621 24 : else if (upper.find("3D8") != std::string::npos)
622 : {
623 24 : elem_type = HEX8;
624 2 : n_nodes_per_elem = 8;
625 4 : elems_of_dimension[3] = true;
626 : }
627 0 : else if (upper.find("3D4") != std::string::npos)
628 : {
629 0 : elem_type = TET4;
630 0 : n_nodes_per_elem = 4;
631 0 : elems_of_dimension[3] = true;
632 : }
633 0 : else if (upper.find("3D20") != std::string::npos)
634 : {
635 0 : elem_type = HEX20;
636 0 : n_nodes_per_elem = 20;
637 0 : elems_of_dimension[3] = true;
638 : }
639 0 : else if (upper.find("3D27") != std::string::npos)
640 : {
641 0 : elem_type = HEX27;
642 0 : n_nodes_per_elem = 27;
643 0 : elems_of_dimension[3] = true;
644 : }
645 0 : else if (upper.find("3D6") != std::string::npos)
646 : {
647 0 : elem_type = PRISM6;
648 0 : n_nodes_per_elem = 6;
649 0 : elems_of_dimension[3] = true;
650 : }
651 0 : else if (upper.find("3D15") != std::string::npos)
652 : {
653 0 : elem_type = PRISM15;
654 0 : n_nodes_per_elem = 15;
655 0 : elems_of_dimension[3] = true;
656 : }
657 0 : else if (upper.find("3D10") != std::string::npos)
658 : {
659 0 : elem_type = TET10;
660 0 : n_nodes_per_elem = 10;
661 0 : elems_of_dimension[3] = true;
662 : }
663 0 : else if (upper.find("MASS") != std::string::npos)
664 : {
665 : // "MASS" elements are used to indicate point masses in Abaqus
666 : // models. These are mapped to NODEELEMs in libmesh.
667 0 : elem_type = NODEELEM;
668 0 : n_nodes_per_elem = 1;
669 0 : elems_of_dimension[0] = true;
670 : }
671 : else
672 0 : libmesh_error_msg("Unrecognized element type: " << upper);
673 :
674 : // Insert the elem type we detected into the set of all elem types for this mesh
675 22 : _elem_types.insert(elem_type);
676 :
677 : // Grab a reference to the element definition for this element type
678 24 : const ElementDefinition & eledef = eletypes[elem_type];
679 :
680 : // If the element definition was not found, the call above would have
681 : // created one with an uninitialized struct. Check for that here...
682 24 : libmesh_error_msg_if
683 : (eledef.abaqus_zero_based_node_id_to_libmesh_node_id.size() == 0,
684 : "No Abaqus->LibMesh mapping information for ElemType "
685 : << Utility::enum_to_string(elem_type) << "!");
686 :
687 : // We will read elements until the next line begins with *, since that will be the
688 : // next section.
689 24120 : while (_in->peek() != '*' && _in->peek() != EOF)
690 : {
691 : // Read the element ID, it is the first number on each line. It is
692 : // followed by a comma, so read that also. We will need this ID later
693 : // when we try to assign subdomain IDs
694 24096 : dof_id_type abaqus_elem_id = 0;
695 : char c;
696 24096 : *_in >> abaqus_elem_id >> c;
697 :
698 : // Add an element of the appropriate type to the Mesh, with the
699 : // abaqus element ID.
700 26104 : std::unique_ptr<Elem> new_elem = Elem::build(elem_type);
701 24096 : new_elem->set_id() = abaqus_elem_id-1;
702 :
703 28112 : Elem * elem = the_mesh.add_elem(std::move(new_elem));
704 :
705 : // The count of the total number of IDs read for the current element.
706 2008 : unsigned id_count=0;
707 :
708 : // Continue reading line-by-line until we have read enough nodes for this element
709 48192 : while (id_count < n_nodes_per_elem)
710 : {
711 : // Read entire line (up to carriage return) of comma-separated values
712 4016 : std::string csv_line;
713 24096 : std::getline(*_in, csv_line);
714 :
715 : // Create a stream object out of the current line
716 28112 : std::stringstream line_stream(csv_line);
717 :
718 : // Process the comma-separated values
719 4016 : std::string cell;
720 216864 : while (std::getline(line_stream, cell, ','))
721 : {
722 : dof_id_type abaqus_global_node_id;
723 192768 : bool success = string_to_num(cell, abaqus_global_node_id);
724 :
725 192768 : if (success)
726 : {
727 : // Map the id'th element ID (Abaqus 1-based numbering) to LibMesh numbering
728 192768 : libmesh_error_msg_if(abaqus_global_node_id < 1,
729 : "Invalid Abaqus node ID found");
730 192768 : const dof_id_type libmesh_global_node_id = abaqus_global_node_id-1;
731 :
732 : // Grab the node pointer from the mesh for this ID
733 192768 : Node * node = the_mesh.node_ptr(libmesh_global_node_id);
734 :
735 : // If node_ptr() returns nullptr, it may mean we have not yet read the
736 : // *Nodes section, though I assumed that always came before the *Elements section...
737 192768 : libmesh_error_msg_if
738 : (node == nullptr,
739 : "Error! Mesh::node_ptr() returned nullptr. Either no node exists with ID "
740 : << libmesh_global_node_id
741 : << " or perhaps this input file has *Elements defined before *Nodes?");
742 :
743 : // Note: id_count is the zero-based abaqus (elem local) node index. We therefore map
744 : // it to a libmesh elem local node index using the element definition map
745 : unsigned libmesh_elem_local_node_id =
746 192768 : eledef.abaqus_zero_based_node_id_to_libmesh_node_id[id_count];
747 :
748 : // Set this node pointer within the element.
749 192768 : elem->set_node(libmesh_elem_local_node_id, node);
750 :
751 : // Increment the count of IDs read for this element
752 192768 : id_count++;
753 : } // end if (success)
754 : } // end while getline(',')
755 20080 : } // end while (id_count)
756 :
757 : // Ensure that we read *exactly* as many nodes as we were expecting to, no more.
758 24096 : libmesh_error_msg_if
759 : (id_count != n_nodes_per_elem,
760 : "Error: Needed to read "
761 : << n_nodes_per_elem
762 : << " nodes, but read "
763 : << id_count
764 : << " instead!");
765 :
766 : // If we are recording Elset IDs, add this element to the correct set for later processing.
767 : // Make sure to add it with the Abaqus ID, not the libmesh one!
768 24096 : if (elset_name != "")
769 0 : _elemset_ids[elset_name].push_back(abaqus_elem_id);
770 20080 : } // end while (peek)
771 24 : }
772 :
773 :
774 :
775 :
776 348 : std::string AbaqusIO::parse_label(std::string line, std::string label_name) const
777 : {
778 : // Handle files which have weird line endings from e.g. windows.
779 : // You can check what kind of line endings you have with 'cat -vet'.
780 : // For example, some files may have two kinds of line endings like:
781 : //
782 : // 4997,^I496,^I532,^I487,^I948^M$
783 : //
784 : // and we don't want to deal with this when extracting a label, so
785 : // just remove all the space characters, which should include all
786 : // kinds of remaining newlines. (I don't think Abaqus allows
787 : // whitespace in label names.)
788 348 : strip_ws(line);
789 :
790 : // Do all string comparisons in upper-case
791 : std::string
792 348 : upper_line(line),
793 348 : upper_label_name(label_name);
794 29 : std::transform(upper_line.begin(), upper_line.end(), upper_line.begin(), ::toupper);
795 29 : std::transform(upper_label_name.begin(), upper_label_name.end(), upper_label_name.begin(), ::toupper);
796 :
797 : // Get index of start of "label="
798 348 : size_t label_index = upper_line.find(upper_label_name + "=");
799 :
800 348 : if (label_index != std::string::npos)
801 : {
802 : // Location of the first comma following "label="
803 275 : size_t comma_index = upper_line.find(",", label_index);
804 :
805 : // Construct iterators from which to build the sub-string.
806 : // Note: The +1 while initializing beg is to skip past the "=" which follows the label name
807 : std::string::iterator
808 50 : beg = line.begin() + label_name.size() + 1 + label_index,
809 300 : end = (comma_index == std::string::npos) ? line.end() : line.begin() + comma_index;
810 :
811 275 : return std::string(beg, end);
812 : }
813 :
814 : // The label index was not found, return the empty string
815 48 : return std::string("");
816 : }
817 :
818 :
819 :
820 :
821 300 : bool AbaqusIO::detect_generated_set(std::string upper) const
822 : {
823 : // Avoid issues with weird line endings, spaces before commas, etc.
824 300 : strip_ws(upper);
825 :
826 : // Check each comma-separated value in "upper" to see if it is the generate flag.
827 50 : std::string cell;
828 350 : std::stringstream line_stream(upper);
829 1008 : while (std::getline(line_stream, cell, ','))
830 876 : if (cell == "GENERATE")
831 14 : return true;
832 :
833 11 : return false;
834 250 : }
835 :
836 :
837 :
838 132 : void AbaqusIO::read_ids(std::string set_name, container_t & container)
839 : {
840 : // Grab a reference to a vector that will hold all the IDs
841 132 : std::vector<dof_id_type> & id_storage = container[set_name];
842 :
843 : // Read until the start of another section is detected, or EOF is encountered
844 1176 : while (_in->peek() != '*' && _in->peek() != EOF)
845 : {
846 : // Read entire comma-separated line into a string
847 174 : std::string csv_line;
848 1044 : std::getline(*_in, csv_line);
849 :
850 : // On that line, use std::getline again to parse each
851 : // comma-separated entry.
852 174 : std::string cell;
853 1218 : std::stringstream line_stream(csv_line);
854 17208 : while (std::getline(line_stream, cell, ','))
855 : {
856 : dof_id_type id;
857 16164 : bool success = string_to_num(cell, id);
858 :
859 : // Note that lists of comma-separated values in abaqus also
860 : // *end* with a comma, so the last call to getline on a given
861 : // line will get an empty string, which we must detect.
862 16164 : if (success)
863 16164 : id_storage.push_back(id);
864 : }
865 870 : }
866 132 : }
867 :
868 :
869 :
870 :
871 168 : void AbaqusIO::generate_ids(std::string set_name, container_t & container)
872 : {
873 : // Grab a reference to a vector that will hold all the IDs
874 168 : std::vector<dof_id_type> & id_storage = container[set_name];
875 :
876 : // Read until the start of another section is detected, or EOF is
877 : // encountered. "generate" sections seem to only have one line,
878 : // although I suppose it's possible they could have more.
879 336 : while (_in->peek() != '*' && _in->peek() != EOF)
880 : {
881 : // Read entire comma-separated line into a string
882 28 : std::string csv_line;
883 168 : std::getline(*_in, csv_line);
884 :
885 : // Remove all whitespaces from csv_line.
886 168 : strip_ws(csv_line);
887 :
888 : // Create a new stringstream object from the string, and stream
889 : // in the comma-separated values.
890 : char c;
891 : dof_id_type start, end, stride;
892 196 : std::stringstream line_stream(csv_line);
893 322 : line_stream >> start >> c >> end >> c >> stride;
894 :
895 : // Generate entries in the id_storage. Note: each element can
896 : // only belong to a single Elset (since this corresponds to the
897 : // subdomain_id) so if an element appears in multiple Elsets,
898 : // the "last" one (alphabetically, based on set name) in the
899 : // _elemset_ids map will "win".
900 29160 : for (dof_id_type current = start; current <= end; current += stride)
901 28992 : id_storage.push_back(current);
902 140 : }
903 168 : }
904 :
905 :
906 :
907 :
908 0 : void AbaqusIO::read_sideset(const std::string & sideset_name,
909 : const std::string & sideset_type,
910 : sideset_container_t & container)
911 : {
912 : // Grab a reference to a vector that will hold all the IDs
913 0 : std::vector<std::pair<dof_id_type, unsigned>> & id_storage = container[sideset_name];
914 :
915 : // Variables for storing values read in from file
916 0 : unsigned side_id=0;
917 : char c;
918 0 : std::string elem_id_or_set, dummy;
919 :
920 : // Read until the start of another section is detected, or EOF is encountered
921 0 : while (_in->peek() != '*' && _in->peek() != EOF)
922 : {
923 : // Read first string up to and including the comma, which is discarded.
924 0 : std::getline(*_in, elem_id_or_set, ',');
925 :
926 : // Strip any leading or trailing trailing whitespace from this
927 : // string, since some Abaqus files may have this.
928 0 : strip_ws(elem_id_or_set);
929 :
930 : // Handle sidesets of type "NODE"
931 0 : if (sideset_type == "NODE")
932 : {
933 : // Create a sideset from a nodeset. This is not currently
934 : // supported, so print a warning and keep going.
935 0 : if (_nodeset_ids.count(elem_id_or_set))
936 0 : libMesh::out << "Warning: skipped creating a sideset from a "
937 0 : << "nodeset, this is not yet supported."
938 0 : << std::endl;
939 : }
940 :
941 : // Otherwise, we assume it's a sideset of the form: "391, S2" or "Elset_1, S3".
942 : // If it's not one of these forms, we'll throw an error instead
943 : // of letting the stream get into a bad state.
944 : else // sideset_type != "NODE"
945 : {
946 : // Read the character "S", followed by the side id. Note: the >> operator
947 : // eats whitespace until it reaches a valid character, so this should work
948 : // whether or not there is a space after the previous comma.
949 0 : *_in >> c >> side_id;
950 :
951 : // Try to convert first string to an integer.
952 : dof_id_type elem_id;
953 0 : bool success = string_to_num(elem_id_or_set, elem_id);
954 :
955 0 : if (success)
956 : {
957 : // if the side set is of the form of "391, S2"
958 0 : id_storage.emplace_back(elem_id, side_id);
959 : }
960 :
961 0 : if (!success)
962 : {
963 : // if the side set is of the form of "Elset_1, S3"
964 0 : const auto & vec = libmesh_map_find(_elemset_ids, elem_id_or_set);
965 0 : for (const auto & elem_id_in_elset : vec)
966 0 : id_storage.emplace_back(elem_id_in_elset, side_id);
967 : }
968 : }
969 :
970 : // Successful or not, we extract the remaining characters on the
971 : // line, including the newline, to (hopefully) go to the next section.
972 0 : std::getline(*_in, dummy);
973 : } // while
974 0 : }
975 :
976 :
977 :
978 :
979 24 : void AbaqusIO::assign_subdomain_ids()
980 : {
981 : // Get a reference to the mesh we are reading
982 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
983 :
984 : // The number of elemsets we've found while reading
985 2 : std::size_t n_elemsets = _elemset_ids.size();
986 :
987 : // Fill in a temporary map with (ElemType, index) pairs based on the _elem_types set. This
988 : // will allow us to easily look up this index in the loop below.
989 4 : std::map<ElemType, unsigned> elem_types_map;
990 : {
991 2 : unsigned ctr=0;
992 48 : for (const auto & type : _elem_types)
993 24 : elem_types_map[type] = ctr++;
994 : }
995 :
996 : // Loop over each Elemset and assign subdomain IDs to Mesh elements
997 : {
998 : // The maximum element dimension seen while reading the Mesh
999 24 : unsigned char max_dim = this->max_elem_dimension_seen();
1000 :
1001 : // The elemset_id counter assigns a logical numbering to the _elemset_ids keys
1002 2 : container_t::iterator it = _elemset_ids.begin(), end = _elemset_ids.end();
1003 240 : for (unsigned elemset_id=0; it != end; ++it, ++elemset_id)
1004 : {
1005 : // Grab a reference to the vector of IDs
1006 18 : std::vector<dof_id_type> & id_vector = it->second;
1007 :
1008 : // Loop over this vector
1009 18168 : for (const auto & id : id_vector)
1010 : {
1011 : // Map the id'th element ID (Abaqus 1-based numbering) to LibMesh numbering
1012 17952 : libmesh_error_msg_if(id < 1,
1013 : "Invalid Abaqus element ID found");
1014 17952 : const dof_id_type libmesh_elem_id = id-1;
1015 :
1016 : // Get reference to that element
1017 17952 : Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
1018 :
1019 : // We won't assign subdomain ids to lower-dimensional
1020 : // elements, as they are assumed to represent boundary
1021 : // conditions. Since lower-dimensional elements can
1022 : // appear in multiple sidesets, it doesn't make sense to
1023 : // assign them a single subdomain id... only the last one
1024 : // assigned would actually "stick".
1025 17952 : if (elem.dim() < max_dim)
1026 0 : break;
1027 :
1028 : // Compute the proper subdomain ID, based on the formula in the
1029 : // documentation for this function.
1030 : subdomain_id_type computed_id = cast_int<subdomain_id_type>
1031 17952 : (elemset_id + (elem_types_map[elem.type()] * n_elemsets));
1032 :
1033 : // Assign this ID to the element in question
1034 17952 : elem.subdomain_id() = computed_id;
1035 :
1036 : // We will also assign a unique name to the computed_id,
1037 : // which is created by appending the geometric element
1038 : // name to the elset name provided by the user in the
1039 : // Abaqus file.
1040 37400 : std::string computed_name = it->first + "_" + Utility::enum_to_string(elem.type());
1041 17952 : the_mesh.subdomain_name(computed_id) = computed_name;
1042 : }
1043 : }
1044 : }
1045 24 : }
1046 :
1047 :
1048 :
1049 :
1050 24 : void AbaqusIO::assign_boundary_node_ids()
1051 : {
1052 : // Get a reference to the mesh we are reading
1053 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
1054 :
1055 : // Iterate over the container of nodesets
1056 2 : container_t::iterator it = _nodeset_ids.begin(), end = _nodeset_ids.end();
1057 108 : for (boundary_id_type current_id=0; it != end; ++it, ++current_id)
1058 : {
1059 : // Associate current_id with the name we determined earlier
1060 91 : the_mesh.get_boundary_info().nodeset_name(current_id) = it->first;
1061 :
1062 : // Get a reference to the current vector of nodeset ID values
1063 7 : std::vector<dof_id_type> & nodeset_ids = it->second;
1064 :
1065 27288 : for (const auto & id : nodeset_ids)
1066 : {
1067 : // Map the id'th element ID (Abaqus 1-based numbering) to LibMesh numbering
1068 27204 : libmesh_error_msg_if(id < 1,
1069 : "Invalid Abaqus node ID found");
1070 27204 : const dof_id_type libmesh_global_node_id = id-1;
1071 :
1072 : // Get node pointer from the mesh
1073 27204 : Node * node = the_mesh.node_ptr(libmesh_global_node_id);
1074 :
1075 27204 : libmesh_error_msg_if(node == nullptr,
1076 : "Error! Mesh::node_ptr() returned nullptr!");
1077 :
1078 : // Add this node with the current_id (which is determined by the
1079 : // alphabetical ordering of the map) to the BoundaryInfo object
1080 27204 : the_mesh.get_boundary_info().add_node(node, current_id);
1081 : }
1082 : }
1083 24 : }
1084 :
1085 :
1086 :
1087 :
1088 24 : void AbaqusIO::assign_sideset_ids()
1089 : {
1090 : // Get a reference to the mesh we are reading
1091 24 : MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
1092 :
1093 : // initialize the eletypes map (eletypes is a file-global variable)
1094 24 : init_eletypes();
1095 :
1096 : // Iterate over the container of sidesets
1097 : {
1098 2 : sideset_container_t::iterator it = _sideset_ids.begin(), end = _sideset_ids.end();
1099 24 : for (boundary_id_type current_id=0; it != end; ++it, ++current_id)
1100 : {
1101 : // Associate current_id with the name we determined earlier
1102 0 : the_mesh.get_boundary_info().sideset_name(current_id) = it->first;
1103 :
1104 : // Get a reference to the current vector of nodeset ID values
1105 0 : std::vector<std::pair<dof_id_type,unsigned>> & sideset_ids = it->second;
1106 :
1107 0 : for (const auto & [abaqus_elem_id, abaqus_side_number] : sideset_ids)
1108 : {
1109 : // Map the id'th element ID (Abaqus 1-based numbering) to LibMesh numbering
1110 0 : libmesh_error_msg_if(abaqus_elem_id < 1,
1111 : "Invalid Abaqus element ID found");
1112 0 : const dof_id_type libmesh_elem_id = abaqus_elem_id-1;
1113 :
1114 : // Get a reference to that element
1115 0 : Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
1116 :
1117 : // Grab a reference to the element definition for this element type
1118 0 : const ElementDefinition & eledef = eletypes[elem.type()];
1119 :
1120 : // If the element definition was not found, the call above would have
1121 : // created one with an uninitialized struct. Check for that here...
1122 0 : libmesh_error_msg_if(eledef.abaqus_zero_based_side_id_to_libmesh_side_id.size() == 0,
1123 : "No Abaqus->LibMesh mapping information for ElemType "
1124 : << Utility::enum_to_string(elem.type()) << "!");
1125 :
1126 : // Add this node with the current_id (which is determined by the
1127 : // alphabetical ordering of the map). Side numbers in Abaqus are 1-based,
1128 : // so we subtract 1 here before passing the abaqus side number to the
1129 : // mapping array
1130 0 : the_mesh.get_boundary_info().add_side
1131 0 : (&elem,
1132 0 : eledef.abaqus_zero_based_side_id_to_libmesh_side_id[abaqus_side_number-1],
1133 : current_id);
1134 : }
1135 : }
1136 : }
1137 :
1138 :
1139 : // Some elsets (if they contain lower-dimensional elements) also
1140 : // define sidesets. So loop over them and build a searchable data
1141 : // structure we can use to assign sidesets.
1142 : {
1143 24 : unsigned char max_dim = this->max_elem_dimension_seen();
1144 :
1145 : // multimap from lower-dimensional-element-hash-key to
1146 : // pair(lower-dimensional-element, boundary_id). The
1147 : // lower-dimensional element is used to verify the results of the
1148 : // hash table search. The boundary_id will be used to set a
1149 : // boundary ID on a higher-dimensional element. We use a multimap
1150 : // because the lower-dimensional elements can belong to more than
1151 : // 1 sideset, and multiple lower-dimensional elements can hash to
1152 : // the same value, but this is very rare.
1153 4 : std::unordered_multimap<dof_id_type, std::pair<Elem *, boundary_id_type>> provide_bcs;
1154 :
1155 : // The elemset_id counter assigns a logical numbering to the
1156 : // _elemset_ids keys. We are going to use these ids as boundary
1157 : // ids, so elemset_id is of type boundary_id_type.
1158 2 : container_t::iterator it = _elemset_ids.begin(), end = _elemset_ids.end();
1159 240 : for (boundary_id_type elemset_id=0; it != end; ++it, ++elemset_id)
1160 : {
1161 : // Grab a reference to the vector of IDs
1162 18 : std::vector<dof_id_type> & id_vector = it->second;
1163 :
1164 : // Loop over this vector
1165 216 : for (const auto & id : id_vector)
1166 : {
1167 : // Map the id'th element ID (Abaqus 1-based numbering) to LibMesh numbering
1168 216 : libmesh_error_msg_if(id < 1,
1169 : "Invalid Abaqus element ID found");
1170 216 : const dof_id_type libmesh_elem_id = id-1;
1171 :
1172 : // Get a reference to that element
1173 216 : Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
1174 :
1175 : // If the element dimension is equal to the maximum
1176 : // dimension seen, we can break out of this for loop --
1177 : // this elset does not contain sideset information.
1178 216 : if (elem.dim() == max_dim)
1179 18 : break;
1180 :
1181 : // If the element dimension is zero, this elset contains
1182 : // NodeElems which AFAIK are not used for sidesets, they
1183 : // are only used to define point masses.
1184 0 : if (elem.dim() == 0)
1185 0 : break;
1186 :
1187 : // We can only handle elements that are *exactly*
1188 : // one dimension lower than the max element
1189 : // dimension. Not sure if "edge" BCs in 3D
1190 : // actually make sense/are required...
1191 0 : libmesh_error_msg_if(elem.dim()+1 != max_dim,
1192 : "ERROR: Expected boundary element of dimension " << max_dim-1 << " but got " << elem.dim());
1193 :
1194 : // Insert the current (key, pair(elem,id)) into the multimap.
1195 0 : provide_bcs.emplace(elem.key(), std::make_pair(&elem, elemset_id));
1196 :
1197 : // Associate the name of this sideset with the ID we've
1198 : // chosen. It's not necessary to do this for every
1199 : // element in the set, but it's convenient to do it here
1200 : // since we have all the necessary information...
1201 0 : the_mesh.get_boundary_info().sideset_name(elemset_id) = it->first;
1202 : }
1203 : }
1204 :
1205 : // Loop over elements and try to assign boundary information
1206 44222 : for (auto & elem : the_mesh.active_element_ptr_range())
1207 24096 : if (elem->dim() == max_dim)
1208 170680 : for (auto sn : elem->side_index_range())
1209 : {
1210 : // This is a max-dimension element that may require BCs.
1211 : // For each of its sides, including internal sides, we'll
1212 : // see if a lower-dimensional element provides boundary
1213 : // information for it. Note that we have not yet called
1214 : // find_neighbors(), so we can't use elem->neighbor(sn) in
1215 : // this algorithm...
1216 144576 : auto bounds = provide_bcs.equal_range (elem->key(sn));
1217 :
1218 : // Add boundary information for each side in the range.
1219 144576 : for (const auto & pr : as_range(bounds))
1220 : {
1221 : // We'll need to compare the lower dimensional element against the current side.
1222 0 : std::unique_ptr<Elem> side (elem->build_side_ptr(sn));
1223 :
1224 : // Extract the relevant data. We don't need the key for anything.
1225 0 : Elem * lower_dim_elem = pr.second.first;
1226 0 : boundary_id_type bid = pr.second.second;
1227 :
1228 : // This was a hash, so it might not be perfect. Let's verify...
1229 0 : if (*lower_dim_elem == *side)
1230 0 : the_mesh.get_boundary_info().add_side(elem, sn, bid);
1231 0 : }
1232 20 : }
1233 : }
1234 24 : }
1235 :
1236 :
1237 :
1238 348 : void AbaqusIO::process_and_discard_comments()
1239 : {
1240 58 : std::string dummy;
1241 : while (true)
1242 : {
1243 : // We assume we are at the beginning of a line that may be
1244 : // comments or may be data. We need to only discard the line if
1245 : // it begins with **, but we must avoid calling std::getline()
1246 : // since there's no way to put that back.
1247 348 : if (_in->peek() == '*')
1248 : {
1249 : // The first character was a star, so actually read it from the stream.
1250 0 : _in->get();
1251 :
1252 : // Peek at the next character...
1253 0 : if (_in->peek() == '*')
1254 : {
1255 : // OK, second character was star also, by definition this
1256 : // line must be a comment! Read the rest of the line and discard!
1257 0 : std::getline(*_in, dummy);
1258 : }
1259 : else
1260 : {
1261 : // The second character was _not_ a star, so put back the first star
1262 : // we pulled out so that the line can be parsed correctly by somebody
1263 : // else!
1264 0 : _in->unget();
1265 :
1266 : // Finally, break out of the while loop, we are done parsing comments
1267 0 : break;
1268 : }
1269 : }
1270 : else
1271 : {
1272 : // First character was not *, so this line must be data! Break out of the
1273 : // while loop!
1274 29 : break;
1275 : }
1276 : }
1277 348 : }
1278 :
1279 :
1280 :
1281 96 : unsigned char AbaqusIO::max_elem_dimension_seen ()
1282 : {
1283 8 : unsigned char max_dim = 0;
1284 :
1285 : unsigned char elem_dimensions_size = cast_int<unsigned char>
1286 8 : (elems_of_dimension.size());
1287 : // The elems_of_dimension array is 1-based in the UNV reader
1288 384 : for (unsigned char i=1; i<elem_dimensions_size; ++i)
1289 312 : if (elems_of_dimension[i])
1290 8 : max_dim = i;
1291 :
1292 96 : return max_dim;
1293 : }
1294 :
1295 208932 : bool AbaqusIO::string_to_num(const std::string &input, dof_id_type &output) const
1296 : {
1297 : char *endptr;
1298 208932 : output = cast_int<dof_id_type>(std::strtol(input.c_str(), &endptr, /*base=*/10));
1299 :
1300 208932 : return (output != 0 || endptr != input.c_str());
1301 : }
1302 :
1303 34308 : void AbaqusIO::strip_ws(std::string &line) const
1304 : {
1305 28590 : line.erase(std::remove_if(line.begin(), line.end(),
1306 138761 : [](unsigned char const c)
1307 1529029 : { return std::isspace(c); }),
1308 8577 : line.end());
1309 34308 : }
1310 :
1311 : } // namespace
|
