abaqus_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 = 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("B32") != std::string::npos)
517  {
518  elem_type = EDGE3;
519  n_nodes_per_elem = 3;
520  elems_of_dimension[1] = true;
521  }
522  else if (upper.find("S3") != std::string::npos ||
523  upper.find("CPE3") != std::string::npos ||
524  upper.find("2D3") != 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("CPE4") != std::string::npos ||
531  upper.find("S4") != std::string::npos ||
532  upper.find("CPEG4") != std::string::npos ||
533  upper.find("2D4") != std::string::npos)
534  {
535  elem_type = QUAD4;
536  n_nodes_per_elem = 4;
537  elems_of_dimension[2] = true;
538  }
539  else if (upper.find("CPE6") != std::string::npos ||
540  upper.find("S6") != std::string::npos ||
541  upper.find("CPEG6") != std::string::npos ||
542  upper.find("2D6") != std::string::npos)
543  {
544  elem_type = TRI6;
545  n_nodes_per_elem = 6;
546  elems_of_dimension[2] = true;
547  }
548  else if (upper.find("CPE8") != std::string::npos ||
549  upper.find("S8") != std::string::npos ||
550  upper.find("CPEG8") != std::string::npos ||
551  upper.find("2D8") != std::string::npos)
552  {
553  elem_type = QUAD8;
554  n_nodes_per_elem = 8;
555  elems_of_dimension[2] = true;
556  }
557  else if (upper.find("3D8") != std::string::npos)
558  {
559  elem_type = HEX8;
560  n_nodes_per_elem = 8;
561  elems_of_dimension[3] = true;
562  }
563  else if (upper.find("3D4") != std::string::npos)
564  {
565  elem_type = TET4;
566  n_nodes_per_elem = 4;
567  elems_of_dimension[3] = true;
568  }
569  else if (upper.find("3D20") != std::string::npos)
570  {
571  elem_type = HEX20;
572  n_nodes_per_elem = 20;
573  elems_of_dimension[3] = true;
574  }
575  else if (upper.find("3D27") != std::string::npos)
576  {
577  elem_type = HEX27;
578  n_nodes_per_elem = 27;
579  elems_of_dimension[3] = true;
580  }
581  else if (upper.find("3D6") != std::string::npos)
582  {
583  elem_type = PRISM6;
584  n_nodes_per_elem = 6;
585  elems_of_dimension[3] = true;
586  }
587  else if (upper.find("3D15") != std::string::npos)
588  {
589  elem_type = PRISM15;
590  n_nodes_per_elem = 15;
591  elems_of_dimension[3] = true;
592  }
593  else if (upper.find("3D10") != std::string::npos)
594  {
595  elem_type = TET10;
596  n_nodes_per_elem = 10;
597  elems_of_dimension[3] = true;
598  }
599  else
600  libmesh_error_msg("Unrecognized element type: " << upper);
601 
602  // Insert the elem type we detected into the set of all elem types for this mesh
603  _elem_types.insert(elem_type);
604 
605  // Grab a reference to the element definition for this element type
606  const ElementDefinition & eledef = eletypes[elem_type];
607 
608  // If the element definition was not found, the call above would have
609  // created one with an uninitialized struct. Check for that here...
610  if (eledef.abaqus_zero_based_node_id_to_libmesh_node_id.size() == 0)
611  libmesh_error_msg("No Abaqus->LibMesh mapping information for ElemType " \
612  << Utility::enum_to_string(elem_type) \
613  << "!");
614 
615  // We will read elements until the next line begins with *, since that will be the
616  // next section.
617  while (_in.peek() != '*' && _in.peek() != EOF)
618  {
619  // Read the element ID, it is the first number on each line. It is
620  // followed by a comma, so read that also. We will need this ID later
621  // when we try to assign subdomain IDs
622  dof_id_type abaqus_elem_id = 0;
623  char c;
624  _in >> abaqus_elem_id >> c;
625 
626  // Add an element of the appropriate type to the Mesh.
627  Elem * elem = the_mesh.add_elem(Elem::build(elem_type).release());
628 
629  // Associate the ID returned from libmesh with the abaqus element ID
630  //_libmesh_to_abaqus_elem_mapping[elem->id()] = abaqus_elem_id;
631  _abaqus_to_libmesh_elem_mapping[abaqus_elem_id] = elem->id();
632 
633  // The count of the total number of IDs read for the current element.
634  unsigned id_count=0;
635 
636  // Continue reading line-by-line until we have read enough nodes for this element
637  while (id_count < n_nodes_per_elem)
638  {
639  // Read entire line (up to carriage return) of comma-separated values
640  std::string csv_line;
641  std::getline(_in, csv_line);
642 
643  // Create a stream object out of the current line
644  std::stringstream line_stream(csv_line);
645 
646  // Process the comma-separated values
647  std::string cell;
648  while (std::getline(line_stream, cell, ','))
649  {
650  // FIXME: factor out this strtol stuff into a utility function.
651  char * endptr;
652  dof_id_type abaqus_global_node_id = cast_int<dof_id_type>
653  (std::strtol(cell.c_str(), &endptr, /*base=*/10));
654 
655  if (abaqus_global_node_id!=0 || cell.c_str() != endptr)
656  {
657  // Use the global node number mapping to determine the corresponding libmesh global node id
658  dof_id_type libmesh_global_node_id = _abaqus_to_libmesh_node_mapping[abaqus_global_node_id];
659 
660  // Grab the node pointer from the mesh for this ID
661  Node * node = the_mesh.node_ptr(libmesh_global_node_id);
662 
663  // If node_ptr() returns nullptr, it may mean we have not yet read the
664  // *Nodes section, though I assumed that always came before the *Elements section...
665  if (node == nullptr)
666  libmesh_error_msg("Error! Mesh::node_ptr() returned nullptr. Either no node exists with ID " \
667  << libmesh_global_node_id \
668  << " or perhaps this input file has *Elements defined before *Nodes?");
669 
670  // Note: id_count is the zero-based abaqus (elem local) node index. We therefore map
671  // it to a libmesh elem local node index using the element definition map
672  unsigned libmesh_elem_local_node_id =
673  eledef.abaqus_zero_based_node_id_to_libmesh_node_id[id_count];
674 
675  // Set this node pointer within the element.
676  elem->set_node(libmesh_elem_local_node_id) = node;
677 
678  // Increment the count of IDs read for this element
679  id_count++;
680  } // end if strtol success
681  } // end while getline(',')
682  } // end while (id_count)
683 
684  // Ensure that we read *exactly* as many nodes as we were expecting to, no more.
685  if (id_count != n_nodes_per_elem)
686  libmesh_error_msg("Error: Needed to read " \
687  << n_nodes_per_elem \
688  << " nodes, but read " \
689  << id_count \
690  << " instead!");
691 
692  // If we are recording Elset IDs, add this element to the correct set for later processing.
693  // Make sure to add it with the Abaqus ID, not the libmesh one!
694  if (elset_name != "")
695  _elemset_ids[elset_name].push_back(abaqus_elem_id);
696  } // end while (peek)
697 }
698 
699 
700 
701 
702 std::string AbaqusIO::parse_label(std::string line, std::string label_name) const
703 {
704  // Handle files which have weird line endings from e.g. windows.
705  // You can check what kind of line endings you have with 'cat -vet'.
706  // For example, some files may have two kinds of line endings like:
707  //
708  // 4997,^I496,^I532,^I487,^I948^M$
709  //
710  // and we don't want to deal with this when extracting a label, so
711  // just remove all the space characters, which should include all
712  // kinds of remaining newlines. (I don't think Abaqus allows
713  // whitespace in label names.)
714  line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end());
715 
716  // Do all string comparisons in upper-case
717  std::string
718  upper_line(line),
719  upper_label_name(label_name);
720  std::transform(upper_line.begin(), upper_line.end(), upper_line.begin(), ::toupper);
721  std::transform(upper_label_name.begin(), upper_label_name.end(), upper_label_name.begin(), ::toupper);
722 
723  // Get index of start of "label="
724  size_t label_index = upper_line.find(upper_label_name + "=");
725 
726  if (label_index != std::string::npos)
727  {
728  // Location of the first comma following "label="
729  size_t comma_index = upper_line.find(",", label_index);
730 
731  // Construct iterators from which to build the sub-string.
732  // Note: The +1 while initializing beg is to skip past the "=" which follows the label name
733  std::string::iterator
734  beg = line.begin() + label_name.size() + 1 + label_index,
735  end = (comma_index == std::string::npos) ? line.end() : line.begin() + comma_index;
736 
737  return std::string(beg, end);
738  }
739 
740  // The label index was not found, return the empty string
741  return std::string("");
742 }
743 
744 
745 
746 
747 bool AbaqusIO::detect_generated_set(std::string upper) const
748 {
749  // Avoid issues with weird line endings, spaces before commas, etc.
750  upper.erase(std::remove_if(upper.begin(), upper.end(), isspace), upper.end());
751 
752  // Check each comma-separated value in "upper" to see if it is the generate flag.
753  std::string cell;
754  std::stringstream line_stream(upper);
755  while (std::getline(line_stream, cell, ','))
756  if (cell == "GENERATE")
757  return true;
758 
759  return false;
760 }
761 
762 
763 
764 void AbaqusIO::read_ids(std::string set_name, container_t & container)
765 {
766  // Grab a reference to a vector that will hold all the IDs
767  std::vector<dof_id_type> & id_storage = container[set_name];
768 
769  // Read until the start of another section is detected, or EOF is encountered
770  while (_in.peek() != '*' && _in.peek() != EOF)
771  {
772  // Read entire comma-separated line into a string
773  std::string csv_line;
774  std::getline(_in, csv_line);
775 
776  // On that line, use std::getline again to parse each
777  // comma-separated entry.
778  std::string cell;
779  std::stringstream line_stream(csv_line);
780  while (std::getline(line_stream, cell, ','))
781  {
782  // If no conversion can be performed by strtol, 0 is returned.
783  //
784  // If endptr is not nullptr, strtol() stores the address of the
785  // first invalid character in *endptr. If there were no
786  // digits at all, however, strtol() stores the original
787  // value of str in *endptr.
788  char * endptr;
789 
790  // FIXME - this needs to be updated for 64-bit inputs
791  dof_id_type id = cast_int<dof_id_type>
792  (std::strtol(cell.c_str(), &endptr, /*base=*/10));
793 
794  // Note that lists of comma-separated values in abaqus also
795  // *end* with a comma, so the last call to getline on a given
796  // line will get an empty string, which we must detect.
797  if (id != 0 || cell.c_str() != endptr)
798  {
799  // 'cell' is now a string with an integer id in it
800  id_storage.push_back( id );
801  }
802  }
803  }
804 }
805 
806 
807 
808 
809 void AbaqusIO::generate_ids(std::string set_name, container_t & container)
810 {
811  // Grab a reference to a vector that will hold all the IDs
812  std::vector<dof_id_type> & id_storage = container[set_name];
813 
814  // Read until the start of another section is detected, or EOF is
815  // encountered. "generate" sections seem to only have one line,
816  // although I suppose it's possible they could have more.
817  while (_in.peek() != '*' && _in.peek() != EOF)
818  {
819  // Read entire comma-separated line into a string
820  std::string csv_line;
821  std::getline(_in, csv_line);
822 
823  // Remove all whitespaces from csv_line.
824  csv_line.erase(std::remove_if(csv_line.begin(), csv_line.end(), isspace), csv_line.end());
825 
826  // Create a new stringstream object from the string, and stream
827  // in the comma-separated values.
828  char c;
829  dof_id_type start, end, stride;
830  std::stringstream line_stream(csv_line);
831  line_stream >> start >> c >> end >> c >> stride;
832 
833  // Generate entries in the id_storage. Note: each element can
834  // only belong to a single Elset (since this corresponds to the
835  // subdomain_id) so if an element appears in multiple Elsets,
836  // the "last" one (alphabetically, based on set name) in the
837  // _elemset_ids map will "win".
838  for (dof_id_type current = start; current <= end; current += stride)
839  id_storage.push_back(current);
840  }
841 }
842 
843 
844 
845 
846 void AbaqusIO::read_sideset(std::string sideset_name, sideset_container_t & container)
847 {
848  // Grab a reference to a vector that will hold all the IDs
849  std::vector<std::pair<dof_id_type, unsigned>> & id_storage = container[sideset_name];
850 
851  // Variables for storing values read in from file
852  dof_id_type elem_id=0;
853  unsigned side_id=0;
854  char c;
855  std::string dummy;
856 
857  // Read until the start of another section is detected, or EOF is encountered
858  while (_in.peek() != '*' && _in.peek() != EOF)
859  {
860  // The strings are of the form: "391, S2"
861 
862  // Read the element ID and the leading comma
863  _in >> elem_id >> c;
864 
865  // Read another character (the 'S') and finally the side ID
866  _in >> c >> side_id;
867 
868  // Store this pair of data in the vector
869  id_storage.push_back( std::make_pair(elem_id, side_id) );
870 
871  // Extract remaining characters on line including newline
872  std::getline(_in, dummy);
873  } // while
874 }
875 
876 
877 
878 
880 {
881  // Get a reference to the mesh we are reading
882  MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
883 
884  // The number of elemsets we've found while reading
885  std::size_t n_elemsets = _elemset_ids.size();
886 
887  // Fill in a temporary map with (ElemType, index) pairs based on the _elem_types set. This
888  // will allow us to easily look up this index in the loop below.
889  std::map<ElemType, unsigned> elem_types_map;
890  {
891  unsigned ctr=0;
892  for (const auto & type : _elem_types)
893  elem_types_map[type] = ctr++;
894  }
895 
896  // Loop over each Elemset and assign subdomain IDs to Mesh elements
897  {
898  // The maximum element dimension seen while reading the Mesh
899  unsigned char max_dim = this->max_elem_dimension_seen();
900 
901  // The elemset_id counter assigns a logical numbering to the _elemset_ids keys
902  container_t::iterator it = _elemset_ids.begin();
903  for (unsigned elemset_id=0; it != _elemset_ids.end(); ++it, ++elemset_id)
904  {
905  // Grab a reference to the vector of IDs
906  std::vector<dof_id_type> & id_vector = it->second;
907 
908  // Loop over this vector
909  for (const auto & id : id_vector)
910  {
911  // Map the id'th element ID (Abaqus numbering) to LibMesh numbering
912  dof_id_type libmesh_elem_id = _abaqus_to_libmesh_elem_mapping[id];
913 
914  // Get reference to that element
915  Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
916 
917  // We won't assign subdomain ids to lower-dimensional
918  // elements, as they are assumed to represent boundary
919  // conditions. Since lower-dimensional elements can
920  // appear in multiple sidesets, it doesn't make sense to
921  // assign them a single subdomain id... only the last one
922  // assigned would actually "stick".
923  if (elem.dim() < max_dim)
924  break;
925 
926  // Compute the proper subdomain ID, based on the formula in the
927  // documentation for this function.
928  subdomain_id_type computed_id = cast_int<subdomain_id_type>
929  (elemset_id + (elem_types_map[elem.type()] * n_elemsets));
930 
931  // Assign this ID to the element in question
932  elem.subdomain_id() = computed_id;
933 
934  // We will also assign a unique name to the computed_id,
935  // which is created by appending the geometric element
936  // name to the elset name provided by the user in the
937  // Abaqus file.
938  std::string computed_name = it->first + "_" + Utility::enum_to_string(elem.type());
939  the_mesh.subdomain_name(computed_id) = computed_name;
940  }
941  }
942  }
943 }
944 
945 
946 
947 
949 {
950  // Get a reference to the mesh we are reading
951  MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
952 
953  // Iterate over the container of nodesets
954  container_t::iterator it = _nodeset_ids.begin();
955  for (unsigned short current_id=0; it != _nodeset_ids.end(); ++it, ++current_id)
956  {
957  // Associate current_id with the name we determined earlier
958  the_mesh.get_boundary_info().nodeset_name(current_id) = it->first;
959 
960  // Get a reference to the current vector of nodeset ID values
961  std::vector<dof_id_type> & nodeset_ids = it->second;
962 
963  for (const auto & id : nodeset_ids)
964  {
965  // Map the Abaqus global node ID to the libmesh node ID
966  dof_id_type libmesh_global_node_id = _abaqus_to_libmesh_node_mapping[id];
967 
968  // Get node pointer from the mesh
969  Node * node = the_mesh.node_ptr(libmesh_global_node_id);
970 
971  if (node == nullptr)
972  libmesh_error_msg("Error! Mesh::node_ptr() returned nullptr!");
973 
974  // Add this node with the current_id (which is determined by the
975  // alphabetical ordering of the map) to the BoundaryInfo object
976  the_mesh.get_boundary_info().add_node(node, current_id);
977  }
978  }
979 }
980 
981 
982 
983 
985 {
986  // Get a reference to the mesh we are reading
987  MeshBase & the_mesh = MeshInput<MeshBase>::mesh();
988 
989  // initialize the eletypes map (eletypes is a file-global variable)
990  init_eletypes();
991 
992  // Iterate over the container of sidesets
993  {
994  sideset_container_t::iterator it = _sideset_ids.begin();
995  for (unsigned short current_id=0; it != _sideset_ids.end(); ++it, ++current_id)
996  {
997  // Associate current_id with the name we determined earlier
998  the_mesh.get_boundary_info().sideset_name(current_id) = it->first;
999 
1000  // Get a reference to the current vector of nodeset ID values
1001  std::vector<std::pair<dof_id_type,unsigned>> & sideset_ids = it->second;
1002 
1003  for (const auto & id : sideset_ids)
1004  {
1005  // sideset_ids is a vector of pairs (elem id, side id). Pull them out
1006  // now to make the code below more readable.
1007  dof_id_type abaqus_elem_id = id.first;
1008  unsigned abaqus_side_number = id.second;
1009 
1010  // Map the Abaqus element ID to LibMesh numbering
1011  dof_id_type libmesh_elem_id = _abaqus_to_libmesh_elem_mapping[ abaqus_elem_id ];
1012 
1013  // Get a reference to that element
1014  Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
1015 
1016  // Grab a reference to the element definition for this element type
1017  const ElementDefinition & eledef = eletypes[elem.type()];
1018 
1019  // If the element definition was not found, the call above would have
1020  // created one with an uninitialized struct. Check for that here...
1021  if (eledef.abaqus_zero_based_side_id_to_libmesh_side_id.size() == 0)
1022  libmesh_error_msg("No Abaqus->LibMesh mapping information for ElemType " \
1023  << Utility::enum_to_string(elem.type()) \
1024  << "!");
1025 
1026  // Add this node with the current_id (which is determined by the
1027  // alphabetical ordering of the map). Side numbers in Abaqus are 1-based,
1028  // so we subtract 1 here before passing the abaqus side number to the
1029  // mapping array
1030  the_mesh.get_boundary_info().add_side
1031  (&elem,
1032  eledef.abaqus_zero_based_side_id_to_libmesh_side_id[abaqus_side_number-1],
1033  current_id);
1034  }
1035  }
1036  }
1037 
1038 
1039  // Some elsets (if they contain lower-dimensional elements) also
1040  // define sidesets. So loop over them and build a searchable data
1041  // structure we can use to assign sidesets.
1042  {
1043  unsigned char max_dim = this->max_elem_dimension_seen();
1044 
1045  // multimap from lower-dimensional-element-hash-key to
1046  // pair(lower-dimensional-element, boundary_id). The
1047  // lower-dimensional element is used to verify the results of the
1048  // hash table search. The boundary_id will be used to set a
1049  // boundary ID on a higher-dimensional element. We use a multimap
1050  // because the lower-dimensional elements can belong to more than
1051  // 1 sideset, and multiple lower-dimensional elements can hash to
1052  // the same value, but this is very rare.
1053  std::unordered_multimap<dof_id_type, std::pair<Elem *, boundary_id_type>> provide_bcs;
1054 
1055  // The elemset_id counter assigns a logical numbering to the
1056  // _elemset_ids keys. We are going to use these ids as boundary
1057  // ids, so elemset_id is of type boundary_id_type.
1058  container_t::iterator it = _elemset_ids.begin();
1059  for (boundary_id_type elemset_id=0; it != _elemset_ids.end(); ++it, ++elemset_id)
1060  {
1061  // Grab a reference to the vector of IDs
1062  std::vector<dof_id_type> & id_vector = it->second;
1063 
1064  // Loop over this vector
1065  for (const auto & id : id_vector)
1066  {
1067  // Map the id_vector[i]'th element ID (Abaqus numbering) to LibMesh numbering
1068  dof_id_type libmesh_elem_id = _abaqus_to_libmesh_elem_mapping[id];
1069 
1070  // Get a reference to that element
1071  Elem & elem = the_mesh.elem_ref(libmesh_elem_id);
1072 
1073  // If the element dimension is equal to the maximum
1074  // dimension seen, we can break out of this for loop --
1075  // this elset does not contain sideset information.
1076  if (elem.dim() == max_dim)
1077  break;
1078 
1079  // We can only handle elements that are *exactly*
1080  // one dimension lower than the max element
1081  // dimension. Not sure if "edge" BCs in 3D
1082  // actually make sense/are required...
1083  if (elem.dim()+1 != max_dim)
1084  libmesh_error_msg("ERROR: Expected boundary element of dimension " << max_dim-1 << " but got " << elem.dim());
1085 
1086  // Insert the current (key, pair(elem,id)) into the multimap.
1087  provide_bcs.insert(std::make_pair(elem.key(),
1088  std::make_pair(&elem,
1089  elemset_id)));
1090 
1091  // Associate the name of this sideset with the ID we've
1092  // chosen. It's not necessary to do this for every
1093  // element in the set, but it's convenient to do it here
1094  // since we have all the necessary information...
1095  the_mesh.get_boundary_info().sideset_name(elemset_id) = it->first;
1096  }
1097  }
1098 
1099  // Loop over elements and try to assign boundary information
1100  for (auto & elem : the_mesh.active_element_ptr_range())
1101  if (elem->dim() == max_dim)
1102  for (auto sn : elem->side_index_range())
1103  {
1104  // This is a max-dimension element that may require BCs.
1105  // For each of its sides, including internal sides, we'll
1106  // see if a lower-dimensional element provides boundary
1107  // information for it. Note that we have not yet called
1108  // find_neighbors(), so we can't use elem->neighbor(sn) in
1109  // this algorithm...
1110  auto bounds = provide_bcs.equal_range (elem->key(sn));
1111 
1112  // Add boundary information for each side in the range.
1113  for (const auto & pr : as_range(bounds))
1114  {
1115  // We'll need to compare the lower dimensional element against the current side.
1116  std::unique_ptr<Elem> side (elem->build_side_ptr(sn));
1117 
1118  // Extract the relevant data. We don't need the key for anything.
1119  Elem * lower_dim_elem = pr.second.first;
1120  boundary_id_type bid = pr.second.second;
1121 
1122  // This was a hash, so it might not be perfect. Let's verify...
1123  if (*lower_dim_elem == *side)
1124  the_mesh.get_boundary_info().add_side(elem, sn, bid);
1125  }
1126  }
1127  }
1128 }
1129 
1130 
1131 
1133 {
1134  std::string dummy;
1135  while (true)
1136  {
1137  // We assume we are at the beginning of a line that may be
1138  // comments or may be data. We need to only discard the line if
1139  // it begins with **, but we must avoid calling std::getline()
1140  // since there's no way to put that back.
1141  if (_in.peek() == '*')
1142  {
1143  // The first character was a star, so actually read it from the stream.
1144  _in.get();
1145 
1146  // Peek at the next character...
1147  if (_in.peek() == '*')
1148  {
1149  // OK, second character was star also, by definition this
1150  // line must be a comment! Read the rest of the line and discard!
1151  std::getline(_in, dummy);
1152  }
1153  else
1154  {
1155  // The second character was _not_ a star, so put back the first star
1156  // we pulled out so that the line can be parsed correctly by somebody
1157  // else!
1158  _in.unget();
1159 
1160  // Finally, break out of the while loop, we are done parsing comments
1161  break;
1162  }
1163  }
1164  else
1165  {
1166  // First character was not *, so this line must be data! Break out of the
1167  // while loop!
1168  break;
1169  }
1170  }
1171 }
1172 
1173 
1174 
1176 {
1177  unsigned char max_dim = 0;
1178 
1179  unsigned char elem_dimensions_size = cast_int<unsigned char>
1180  (elems_of_dimension.size());
1181  // The elems_of_dimension array is 1-based in the UNV reader
1182  for (unsigned char i=1; i<elem_dimensions_size; ++i)
1183  if (elems_of_dimension[i])
1184  max_dim = i;
1185 
1186  return max_dim;
1187 }
1188 
1189 
1190 } // namespace
virtual void read(const std::string &name) override
Definition: abaqus_io.C:205
void read_sideset(std::string sideset_name, sideset_container_t &container)
Definition: abaqus_io.C:846
void process_and_discard_comments()
Definition: abaqus_io.C:1132
std::size_t n_boundary_conds() const
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2024
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
std::string & nodeset_name(boundary_id_type id)
bool _already_seen_part
Definition: abaqus_io.h:235
std::map< std::string, std::vector< dof_id_type > > container_t
Definition: abaqus_io.h:72
std::vector< bool > elems_of_dimension
Definition: mesh_input.h:97
sideset_container_t _sideset_ids
Definition: abaqus_io.h:199
virtual ~AbaqusIO()
Definition: abaqus_io.C:198
unsigned short int side
Definition: xdr_io.C:50
void read_ids(std::string set_name, container_t &container)
Definition: abaqus_io.C:764
The base class for all geometric element types.
Definition: elem.h:100
IterBase * end
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
bool build_sidesets_from_nodesets
Definition: abaqus_io.h:66
void assign_sideset_ids()
Definition: abaqus_io.C:984
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:131
void read_nodes(std::string nset_name)
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
Base class for Mesh.
Definition: mesh_base.h:77
std::map< dof_id_type, dof_id_type > _abaqus_to_libmesh_elem_mapping
Definition: abaqus_io.h:217
void read_elements(std::string upper, std::string elset_name)
Definition: abaqus_io.C:497
std::map< std::string, std::vector< std::pair< dof_id_type, unsigned > > > sideset_container_t
Definition: abaqus_io.h:79
std::string parse_label(std::string line, std::string label_name) const
Definition: abaqus_io.C:702
void add_node(const Node *node, const boundary_id_type id)
bool detect_generated_set(std::string upper) const
Definition: abaqus_io.C:747
int8_t boundary_id_type
Definition: id_types.h:51
virtual void delete_elem(Elem *e)=0
virtual SimpleRange< element_iterator > element_ptr_range()=0
dof_id_type id() const
Definition: dof_object.h:655
void build_side_list_from_node_list()
void generate_ids(std::string set_name, container_t &container)
Definition: abaqus_io.C:809
virtual Elem * add_elem(Elem *e)=0
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:245
virtual dof_id_type key(const unsigned int s) const =0
void assign_subdomain_ids()
Definition: abaqus_io.C:879
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
void assign_boundary_node_ids()
Definition: abaqus_io.C:948
container_t _nodeset_ids
Definition: abaqus_io.h:197
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:538
void set_mesh_dimension(unsigned char d)
Definition: mesh_base.h:213
virtual void clear()
Definition: mesh_base.C:260
std::string & sideset_name(boundary_id_type id)
std::string enum_to_string(const T e)
AbaqusIO(MeshBase &mesh)
Definition: abaqus_io.C:188
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::set< ElemType > _elem_types
Definition: abaqus_io.h:210
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
std::map< dof_id_type, dof_id_type > _abaqus_to_libmesh_node_mapping
Definition: abaqus_io.h:227
virtual unsigned short dim() const =0
unsigned char max_elem_dimension_seen()
Definition: abaqus_io.C:1175
void swap(Iterator &lhs, Iterator &rhs)
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:504
virtual const Node * node_ptr(const dof_id_type i) const =0
container_t _elemset_ids
Definition: abaqus_io.h:198
virtual ElemType type() const =0
A geometric point in (x,y,z) space.
Definition: point.h:38
virtual dof_id_type n_nodes() const =0
std::ifstream _in
Definition: abaqus_io.h:204
uint8_t dof_id_type
Definition: id_types.h:64