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