abaqus_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
23 // Local includes
24 #include "libmesh/abaqus_io.h"
25 #include "libmesh/point.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/string_to_enum.h"
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/utility.h"
30 
31 // Anonymous namespace to hold mapping Data for Abaqus/libMesh element types
32 namespace
33 {
38  struct ElementDefinition
39  {
40  // Maps (zero-based!) Abaqus local node numbers to libmesh local node numbers
41  std::vector<unsigned> abaqus_zero_based_node_id_to_libmesh_node_id;
42 
43  // Maps (zero-based!) Abaqus side numbers to libmesh side numbers
44  std::vector<unsigned> abaqus_zero_based_side_id_to_libmesh_side_id;
45  };
46 
50  std::map<ElemType, ElementDefinition> eletypes;
51 
55  void add_eletype_entry(ElemType libmesh_elem_type,
56  const unsigned* node_map,
57  unsigned node_map_size,
58  const unsigned* side_map,
59  unsigned side_map_size)
60  {
61  // If map entry does not exist, this will create it
62  ElementDefinition& map_entry = eletypes[libmesh_elem_type];
63 
64 
65  // Use the "swap trick" from Scott Meyer's "Effective STL" to swap
66  // an unnamed temporary vector into the map_entry's vector. Note:
67  // the vector(iter, iter) constructor is used.
68  std::vector<unsigned>(node_map,
69  node_map+node_map_size).swap(map_entry.abaqus_zero_based_node_id_to_libmesh_node_id);
70 
71  std::vector<unsigned>(side_map,
72  side_map+side_map_size).swap(map_entry.abaqus_zero_based_side_id_to_libmesh_side_id);
73  }
74 
75 
79  void init_eletypes ()
80  {
81  // This should happen only once. The first time this method is
82  // called the eletypes data struture will be empty, and we will
83  // fill it. Any subsequent calls will find an initialized
84  // eletypes map and will do nothing.
85  if (eletypes.empty())
86  {
87  {
88  // TRI3
89  const unsigned int node_map[] = {0,1,2}; // identity
90  const unsigned int side_map[] = {0,1,2}; // identity
91  add_eletype_entry(TRI3, node_map, 3, side_map, 3);
92  }
93 
94  {
95  // QUAD4
96  const unsigned int node_map[] = {0,1,2,3}; // identity
97  const unsigned int side_map[] = {0,1,2,3}; // identity
98  add_eletype_entry(QUAD4, node_map, 4, side_map, 4);
99  }
100 
101  {
102  // TET4
103  const unsigned int node_map[] = {0,1,2,3}; // identity
104  const unsigned int side_map[] = {0,1,2,3}; // identity
105  add_eletype_entry(TET4, node_map, 4, side_map, 4);
106  }
107 
108  {
109  // TET10
110  const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,9}; // identity
111  const unsigned int side_map[] = {0,1,2,3}; // identity
112  add_eletype_entry(TET10, node_map, 10, side_map, 4);
113  }
114 
115  {
116  // HEX8
117  const unsigned int node_map[] = {0,1,2,3,4,5,6,7}; // identity
118  const unsigned int side_map[] = {0,5,1,2,3,4}; // inverse = 0,2,3,4,5,1
119  add_eletype_entry(HEX8, node_map, 8, side_map, 6);
120  }
121 
122  {
123  // HEX20
124  const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15}; // map is its own inverse
125  const unsigned int side_map[] = {0,5,1,2,3,4}; // inverse = 0,2,3,4,5,1
126  add_eletype_entry(HEX20, node_map, 20, side_map, 6);
127  }
128 
129  {
130  // HEX27
131  const unsigned int node_map[] = {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}; // inverse = ...,21,23,24,25,26,22,20
132  const unsigned int side_map[] = {0,5,1,2,3,4}; // inverse = 0,2,3,4,5,1
133  add_eletype_entry(HEX27, node_map, 27, side_map, 6);
134  }
135 
136  {
137  // PRISM6
138  const unsigned int node_map[] = {0,1,2,3,4,5}; // identity
139  const unsigned int side_map[] = {0,4,1,2,3}; // inverse = 0,2,3,4,1
140  add_eletype_entry(PRISM6, node_map, 6, side_map, 5);
141  }
142 
143  {
144  // PRISM15
145  const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,12,13,14,9,10,11}; // map is its own inverse
146  const unsigned int side_map[] = {0,4,1,2,3}; // inverse = 0,2,3,4,1
147  add_eletype_entry(PRISM15, node_map, 15, side_map, 5);
148  }
149 
150  {
151  // PRISM18
152  const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,12,13,14,9,10,11,15,16,17}; // map is its own inverse
153  const unsigned int side_map[] = {0,4,1,2,3}; // inverse = 0,2,3,4,1
154  add_eletype_entry(PRISM18, node_map, 18, side_map, 5);
155  }
156 
157 
158 
159  } // if (eletypes.empty())
160  } // init_eletypes()
161 } // anonymous namespace
162 
163 
164 
165 
166 
167 namespace libMesh
168 {
169 
171  MeshInput<MeshBase> (mesh_in),
172  build_sidesets_from_nodesets(false),
173  _already_seen_part(false)
174  {
175  }
176 
177 
178 
179 
181  {
182  }
183 
184 
185 
186 
187  void AbaqusIO::read (const std::string& fname)
188  {
189  // Get a reference to the mesh we are reading
190  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
191 
192  // Clear any existing mesh data
193  the_mesh.clear();
194 
195  // Open stream for reading
196  _in.open(fname.c_str());
197  libmesh_assert(_in.good());
198 
199  // Read file line-by-line... this is based on a set of different
200  // test input files. I have not looked at the full input file
201  // specs for Abaqus.
202  std::string s;
203  while (true)
204  {
205  // Try to read something. This may set EOF!
206  std::getline(_in, s);
207 
208  if (_in)
209  {
210  // Process s...
211  //
212  // There are many sections in Abaqus files, we read some
213  // but others are just ignored... Some sections may occur
214  // more than once. For example for a hybrid grid, you
215  // will have multiple *Element sections...
216 
217  // Some Abaqus files use all upper-case for section names,
218  // so we will just convert s to uppercase
219  std::string upper(s);
220  std::transform(upper.begin(), upper.end(), upper.begin(), ::toupper);
221 
222  // 0.) Look for the "*Part" section
223  if (upper.find("*PART") == 0)
224  {
225  // libMesh::out << "Found parts section!" << std::endl;
226 
227  if (_already_seen_part)
228  {
229  libMesh::err << "We currently don't support reading Abaqus files with multiple PART sections" << std::endl;
230  libmesh_error();
231  }
232 
233  _already_seen_part = true;
234  }
235 
236  // 1.) Look for the "*Nodes" section
237  if (upper.find("*NODE") == 0)
238  {
239  // Process any lines of comments that may be present
241 
242  // Read a block of nodes
243  this->read_nodes();
244  }
245 
246 
247 
248  // 2.) Look for the "*Element" section
249  else if (upper.find("*ELEMENT,") == 0)
250  {
251  // Process any lines of comments that may be present
253 
254  // Read a block of elements
255  this->read_elements(upper);
256  }
257 
258 
259 
260  // 3.) Look for a Nodeset section
261  else if (upper.find("*NSET") == 0)
262  {
263  std::string nset_name = this->parse_label(s, "nset");
264 
265  // I haven't seen an unnamed elset yet, but let's detect it
266  // just in case...
267  if (nset_name == "")
268  {
269  libMesh::err << "Unnamed nset encountered!" << std::endl;
270  libmesh_error();
271  }
272 
273  // Process any lines of comments that may be present
275 
276  // Read the IDs, storing them in _nodeset_ids
277  this->read_ids(nset_name, _nodeset_ids);
278  } // *Nodeset
279 
280 
281 
282  // 4.) Look for an Elset section
283  else if (upper.find("*ELSET") == 0)
284  {
285  std::string elset_name = this->parse_label(s, "elset");
286 
287  // I haven't seen an unnamed elset yet, but let's detect it
288  // just in case...
289  if (elset_name == "")
290  {
291  libMesh::err << "Unnamed elset encountered!" << std::endl;
292  libmesh_error();
293  }
294 
295  // Debugging
296  // libMesh::out << "Processing ELSET: " << elset_name << std::endl;
297 
298  // Process any lines of comments that may be present
300 
301  // Read the IDs, storing them in _elemset_ids
302  this->read_ids(elset_name, _elemset_ids);
303  } // *Elset
304 
305 
306 
307  // 5.) Look for a Surface section. Need to be a little
308  // careful, since there are also "surface interaction"
309  // sections we don't want to read here.
310  else if (upper.find("*SURFACE,") == 0)
311  {
312  // libMesh::out << "Found SURFACE section: " << s << std::endl;
313 
314  // Get the name from the Name=Foo label. This will be the map key.
315  std::string sideset_name = this->parse_label(s, "name");
316 
317  // Print name of section we just found
318  // libMesh::out << "Found surface section named: " << sideset_name << std::endl;
319 
320  // Process any lines of comments that may be present
322 
323  // Read the sideset IDs
324  this->read_sideset(sideset_name, _sideset_ids);
325 
326  // Debugging: print status of most recently read sideset
327  // libMesh::out << "Read " << _sideset_ids[sideset_name].size() << " sides in " << sideset_name << std::endl;
328  }
329 
330  continue;
331  } // if (_in)
332 
333  // If !file, check to see if EOF was set. If so, break out
334  // of while loop.
335  if (_in.eof())
336  break;
337 
338  // If !in and !in.eof(), stream is in a bad state!
339  libMesh::err << "Stream is bad!\n";
340  libMesh::err << "Perhaps the file: " << fname << " does not exist?" << std::endl;
341  libmesh_error();
342  } // while
343 
344 
345  //
346  // We've read everything we can from the file at this point. Now
347  // do some more processing.
348  //
349  libMesh::out << "Mesh contains "
350  << the_mesh.n_elem()
351  << " elements, and "
352  << the_mesh.n_nodes()
353  << " nodes." << std::endl;
354 
355  // TODO: Remove these or write a function to do it?
356 // {
357 // container_t::iterator it=_nodeset_ids.begin();
358 // for (; it != _nodeset_ids.end(); ++it)
359 // {
360 // libMesh::out << "Node set '" << (*it).first << "' contains " << (*it).second.size() << " ID(s)." << std::endl;
361 // }
362 // }
363 //
364 // {
365 // container_t::iterator it=_elemset_ids.begin();
366 // for (; it != _elemset_ids.end(); ++it)
367 // {
368 // libMesh::out << "Elem set '" << (*it).first << "' contains " << (*it).second.size() << " ID(s)." << std::endl;
369 // }
370 // }
371 
372 
373  // Set element IDs based on the element sets.
374  this->assign_subdomain_ids();
375 
376  // Assign nodeset values to the BoundaryInfo object
377  this->assign_boundary_node_ids();
378 
379  // Assign sideset values in the BoundaryInfo object
380  this->assign_sideset_ids();
381 
382  // Abaqus files only contain nodesets by default. To be useful in
383  // applying most types of BCs in libmesh, we will definitely need
384  // sidesets. So we can call the new BoundaryInfo function which
385  // generates sidesets from nodesets.
387  the_mesh.boundary_info->build_side_list_from_node_list();
388 
389  } // read()
390 
391 
392 
393 
394 
395 
396 
398  {
399  // Get a reference to the mesh we are reading
400  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
401 
402  // Debugging: print node count
403  // libMesh::out << "Before read_nodes(), mesh contains "
404  // << the_mesh.n_elem()
405  // << " elements, and "
406  // << the_mesh.n_nodes()
407  // << " nodes." << std::endl;
408 
409  // In the input file I have, Abaqus neither tells what
410  // the mesh dimension is nor how many nodes it has...
411 
412  // The node line format is:
413  // id, x, y, z
414  // and you do have to parse out the commas.
415  // The z-coordinate will only be present for 3D meshes
416 
417  // Temporary variables for parsing lines of text
418  dof_id_type abaqus_node_id=0;
419  Real x=0, y=0, z=0;
420  char c;
421  std::string dummy;
422 
423  // Defines the sequential node numbering used by libmesh
424  dof_id_type libmesh_node_id = 0;
425 
426  // We will read nodes until the next line begins with *, since that will be the
427  // next section.
428  // TODO: Is Abaqus guaranteed to start the line with '*' or can there be leading white space?
429  while (_in.peek() != '*' && _in.peek() != EOF)
430  {
431  // Re-Initialize variables to be read in from file
432  abaqus_node_id=0;
433  x = y = z = 0.;
434 
435  // Note: we assume *at least* 2D points here, should we worry about
436  // trying to read 1D Abaqus meshes?
437  _in >> abaqus_node_id >> c >> x >> c >> y;
438 
439  // Peek at the next character. If it is a comma, then there is another
440  // value to read!
441  if (_in.peek() == ',')
442  _in >> c >> z;
443 
444  // Debugging: Print what we just read in.
445  // libMesh::out << "node_id=" << node_id
446  // << ", x=" << x
447  // << ", y=" << y
448  // << ", z=" << z
449  // << std::endl;
450 
451  // Read (and discard) the rest of the line, including the newline.
452  // This is required so that our 'peek()' at the beginning of this
453  // loop doesn't read the newline character, for example.
454  std::getline(_in, dummy);
455 
456  // Set up the abaqus -> libmesh node mapping. This is usually just the
457  // "off-by-one" map.
458  _abaqus_to_libmesh_node_mapping[abaqus_node_id] = libmesh_node_id;
459 
460  // Add the point to the mesh using libmesh's numbering,
461  // and post-increment the libmesh node counter.
462  the_mesh.add_point(Point(x,y,z), libmesh_node_id++);
463  } // while
464 
465  // Debugging: print node count. Note: in serial mesh, this count may
466  // be off by one, since Abaqus uses one-based numbering, and libmesh
467  // just reports the length of its _nodes vector for the number of nodes.
468  // libMesh::out << "After read_nodes(), mesh contains "
469  // << the_mesh.n_elem()
470  // << " elements, and "
471  // << the_mesh.n_nodes()
472  // << " nodes." << std::endl;
473 
474  } // read_nodes()
475 
476 
477 
478 
479 
480  void AbaqusIO::read_elements(std::string upper)
481  {
482  // Some *Element sections also specify an Elset name on the same line.
483  // Look for one here.
484  std::string elset_name = this->parse_label(upper, "elset");
485 
486  // Get a reference to the mesh we are reading
487  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
488 
489  // initialize the eletypes map (eletypes is a file-global variable)
490  init_eletypes();
491 
492  ElemType elem_type = INVALID_ELEM;
493  unsigned n_nodes_per_elem = 0;
494 
495  // Within s, we should have "type=XXXX"
496  if (upper.find("CPE4") != std::string::npos ||
497  upper.find("CPS4") != std::string::npos)
498  {
499  elem_type = QUAD4;
500  n_nodes_per_elem = 4;
501  the_mesh.set_mesh_dimension(2);
502  }
503  else if (upper.find("CPS3") != std::string::npos)
504  {
505  elem_type = TRI3;
506  n_nodes_per_elem = 3;
507  the_mesh.set_mesh_dimension(2);
508  }
509  else if (upper.find("C3D8") != std::string::npos)
510  {
511  elem_type = HEX8;
512  n_nodes_per_elem = 8;
513  the_mesh.set_mesh_dimension(3);
514  }
515  else if (upper.find("C3D4") != std::string::npos)
516  {
517  elem_type = TET4;
518  n_nodes_per_elem = 4;
519  the_mesh.set_mesh_dimension(3);
520  }
521  else if (upper.find("C3D20") != std::string::npos)
522  {
523  elem_type = HEX20;
524  n_nodes_per_elem = 20;
525  the_mesh.set_mesh_dimension(3);
526  }
527  else if (upper.find("C3D6") != std::string::npos)
528  {
529  elem_type = PRISM6;
530  n_nodes_per_elem = 6;
531  the_mesh.set_mesh_dimension(3);
532  }
533  else if (upper.find("C3D15") != std::string::npos)
534  {
535  elem_type = PRISM15;
536  n_nodes_per_elem = 15;
537  the_mesh.set_mesh_dimension(3);
538  }
539  else if (upper.find("C3D10") != std::string::npos)
540  {
541  elem_type = TET10;
542  n_nodes_per_elem = 10;
543  the_mesh.set_mesh_dimension(3);
544  }
545  else
546  {
547  libMesh::err << "Unrecognized element type: " << upper << std::endl;
548  libmesh_error();
549  }
550 
551  // Insert the elem type we detected into the set of all elem types for this mesh
552  _elem_types.insert(elem_type);
553 
554  // For reading in line endings
555  std::string dummy;
556 
557  // Grab a reference to the element definition for this element type
558  const ElementDefinition& eledef = eletypes[elem_type];
559 
560  // If the element definition was not found, the call above would have
561  // created one with an uninitialized struct. Check for that here...
562  if (eledef.abaqus_zero_based_node_id_to_libmesh_node_id.size() == 0)
563  {
564  // libMesh::err << "No Abaqus->LibMesh mapping information for ElemType "
565  // << Utility::enum_to_string(elem_type)
566  // << "!"
567  // << std::endl;
568  libmesh_error();
569  }
570 
571  // We will read elements until the next line begins with *, since that will be the
572  // next section.
573  while (_in.peek() != '*' && _in.peek() != EOF)
574  {
575  // Read the element ID, it is the first number on each line. It is
576  // followed by a comma, so read that also. We will need this ID later
577  // when we try to assign subdomain IDs
578  dof_id_type abaqus_elem_id = 0;
579  char c;
580  _in >> abaqus_elem_id >> c;
581 
582  // Debugging:
583  // libMesh::out << "Reading data for element " << abaqus_elem_id << std::endl;
584 
585  // Add an element of the appropriate type to the Mesh.
586  Elem* elem = the_mesh.add_elem(Elem::build(elem_type).release());
587 
588  // Associate the ID returned from libmesh with the abaqus element ID
589  //_libmesh_to_abaqus_elem_mapping[elem->id()] = abaqus_elem_id;
590  _abaqus_to_libmesh_elem_mapping[abaqus_elem_id] = elem->id();
591 
592  // The count of the total number of IDs read for the current element.
593  unsigned id_count=0;
594 
595  // Continue reading line-by-line until we have read enough nodes for this element
596  while (id_count < n_nodes_per_elem)
597  {
598  // Read entire line (up to carriage return) of comma-separated values
599  std::string csv_line;
600  std::getline(_in, csv_line);
601 
602  // Create a stream object out of the current line
603  std::stringstream line_stream(csv_line);
604 
605  // Process the comma-separated values
606  std::string cell;
607  while (std::getline(line_stream, cell, ','))
608  {
609  // FIXME: factor out this strtol stuff into a utility function.
610  char* endptr;
611  long abaqus_global_node_id = std::strtol(cell.c_str(), &endptr, /*base=*/10);
612 
613  if (abaqus_global_node_id!=0 || cell.c_str() != endptr)
614  {
615  // Use the global node number mapping to determine the corresponding libmesh global node id
616  dof_id_type libmesh_global_node_id = _abaqus_to_libmesh_node_mapping[abaqus_global_node_id];
617 
618  // Grab the node pointer from the mesh for this ID
619  Node* node = the_mesh.node_ptr(libmesh_global_node_id);
620 
621  // Debugging:
622  // libMesh::out << "Assigning global node id: " << abaqus_global_node_id
623  // << "(Abaqus), " << node->id() << "(LibMesh)" << std::endl;
624 
625  // If node_ptr() returns NULL, it may mean we have not yet read the
626  // *Nodes section, though I assumed that always came before the *Elements section...
627  if (node == NULL)
628  {
629  libMesh::err << "Error! Mesh returned NULL Node pointer.\n";
630  libMesh::err << "Either no node exists with ID " << libmesh_global_node_id
631  << " or perhaps this input file has *Elements defined before *Nodes?" << std::endl;
632  libmesh_error();
633  }
634 
635  // Note: id_count is the zero-based abaqus (elem local) node index. We therefore map
636  // it to a libmesh elem local node index using the element definition map
637  unsigned libmesh_elem_local_node_id =
638  eledef.abaqus_zero_based_node_id_to_libmesh_node_id[id_count];
639 
640  // Set this node pointer within the element.
641  elem->set_node(libmesh_elem_local_node_id) = node;
642 
643  // Debugging:
644  // libMesh::out << "Setting elem " << elem->id()
645  // << ", local node " << libmesh_elem_local_node_id
646  // << " to global node " << node->id() << std::endl;
647 
648  // Increment the count of IDs read for this element
649  id_count++;
650  } // end if strtol success
651  } // end while getline(',')
652  } // end while (id_count)
653 
654  // Ensure that we read *exactly* as many nodes as we were expecting to, no more.
655  if (id_count != n_nodes_per_elem)
656  {
657  libMesh::err << "Error: Needed to read "
658  << n_nodes_per_elem
659  << " nodes, but read "
660  << id_count
661  << " instead!" << std::endl;
662  libmesh_error();
663  }
664 
665  // If we are recording Elset IDs, add this element to the correct set for later processing.
666  // Make sure to add it with the Abaqus ID, not the libmesh one!
667  if (elset_name != "")
668  {
669  // Debugging:
670  // libMesh::out << "Adding Elem " << abaqus_elem_id << " to Elmset " << elset_name << std::endl;
671  _elemset_ids[elset_name].push_back(abaqus_elem_id);
672  }
673  } // end while (peek)
674  } // read_elements()
675 
676 
677 
678 
679  std::string AbaqusIO::parse_label(std::string line, std::string label_name)
680  {
681  // Do all string comparisons in upper-case
682  std::string upper_line(line), upper_label_name(label_name);
683  std::transform(upper_line.begin(), upper_line.end(), upper_line.begin(), ::toupper);
684  std::transform(upper_label_name.begin(), upper_label_name.end(), upper_label_name.begin(), ::toupper);
685 
686  // Get index of start of "label="
687  size_t label_index = upper_line.find(upper_label_name + "=");
688 
689  if (label_index != std::string::npos)
690  {
691  // Location of the first comma following "label="
692  size_t comma_index = upper_line.find(",", label_index);
693 
694  // Construct iterators from which to build the sub-string.
695  // Note the +1 is to skip past the "=" which follows the label name
696  std::string::iterator
697  beg = line.begin() + label_name.size() + 1 + label_index,
698  end = (comma_index == std::string::npos) ? line.end() : line.begin()+comma_index;
699 
700  return std::string(beg, end);
701  }
702 
703  // The label index was not found, return the empty string
704  return std::string("");
705  } // parse_label()
706 
707 
708 
709 
710  void AbaqusIO::read_ids(std::string set_name, container_t& container)
711  {
712  // Debugging
713  // libMesh::out << "Reading ids for set: " << set_name << std::endl;
714 
715  // Grab a reference to a vector that will hold all the IDs
716  std::vector<dof_id_type>& id_storage = container[set_name];
717 
718  // Read until the start of another section is detected, or EOF is encountered
719  while (_in.peek() != '*' && _in.peek() != EOF)
720  {
721  // Read entire comma-separated line into a string
722  std::string csv_line;
723  std::getline(_in, csv_line);
724 
725  // On that line, use std::getline again to parse each
726  // comma-separated entry.
727  std::string cell;
728  std::stringstream line_stream(csv_line);
729  while (std::getline(line_stream, cell, ','))
730  {
731  // If no conversion can be performed by strtol, 0 is returned.
732  //
733  // If endptr is not NULL, strtol() stores the address of the
734  // first invalid character in *endptr. If there were no
735  // digits at all, however, strtol() stores the original
736  // value of str in *endptr.
737  char* endptr;
738 
739  // FIXME - this needs to be updated for 64-bit inputs
740  long id = std::strtol(cell.c_str(), &endptr, /*base=*/10);
741 
742  // Note that lists of comma-separated values in abaqus also
743  // *end* with a comma, so the last call to getline on a given
744  // line will get an empty string, which we must detect.
745  if (id!=0 || cell.c_str() != endptr)
746  {
747  // Debugging
748  // libMesh::out << "Read id: " << id << std::endl;
749 
750  // 'cell' is now a string with an integer id in it
751  id_storage.push_back( id );
752  }
753  }
754  }
755 
756  // Status message
757  // libMesh::out << "Read " << id_storage.size() << " ID(s) for the set " << set_name << std::endl;
758  } // read_ids()
759 
760 
761 
762 
763  void AbaqusIO::read_sideset(std::string sideset_name, sideset_container_t& container)
764  {
765  // Grab a reference to a vector that will hold all the IDs
766  std::vector<std::pair<dof_id_type, unsigned> >& id_storage = container[sideset_name];
767 
768  // Variables for storing values read in from file
769  dof_id_type elem_id=0;
770  unsigned side_id=0;
771  char c;
772  std::string dummy;
773 
774  // Read until the start of another section is detected, or EOF is encountered
775  while (_in.peek() != '*' && _in.peek() != EOF)
776  {
777  // The strings are of the form: "391, S2"
778 
779  // Read the element ID and the leading comma
780  _in >> elem_id >> c;
781 
782  // Read another character (the 'S') and finally the side ID
783  _in >> c >> side_id;
784 
785  // Debugging: print status
786  // libMesh::out << "Read elem_id=" << elem_id << ", side_id=" << side_id << std::endl;
787 
788  // Store this pair of data in the vector
789  id_storage.push_back( std::make_pair(elem_id, side_id) );
790 
791  // Extract remaining characters on line including newline
792  std::getline(_in, dummy);
793  } // while
794  }
795 
796 
797 
798 
800  {
801  // Get a reference to the mesh we are reading
802  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
803 
804  // The number of elemsets we've found while reading
805  std::size_t n_elemsets = _elemset_ids.size();
806 
807  // Fill in a temporary map with (ElemType, index) pairs based on the _elem_types set. This
808  // will allow us to easily look up this index in the loop below.
809  std::map<ElemType, unsigned> elem_types_map;
810  {
811  unsigned ctr=0;
812  for (std::set<ElemType>::iterator it=_elem_types.begin(); it!=_elem_types.end(); ++it)
813  elem_types_map[*it] = ctr++;
814  }
815 
816  // Loop over each Elemset and assign subdomain IDs to Mesh elements
817  {
818  // The elemset_id counter assigns a logical numbering to the _elemset_ids keys
819  container_t::iterator it=_elemset_ids.begin();
820  for (unsigned elemset_id=0; it != _elemset_ids.end(); ++it, ++elemset_id)
821  {
822  // Grab a reference to the vector of IDs
823  std::vector<dof_id_type>& id_vector = (*it).second;
824 
825  // Loop over this vector
826  for (std::size_t i=0; i<id_vector.size(); ++i)
827  {
828  // Map the id_vector[i]'th element ID (Abaqus numbering) to LibMesh numbering
829  dof_id_type libmesh_elem_id = _abaqus_to_libmesh_elem_mapping[ id_vector[i] ];
830 
831  // Get pointer to that element
832  Elem* elem = the_mesh.elem(libmesh_elem_id);
833 
834  if (elem == NULL)
835  {
836  libMesh::err << "Mesh returned NULL pointer for Elem " << libmesh_elem_id << std::endl;
837  libmesh_error();
838  }
839 
840  // Compute the proper subdomain ID, based on the formula in the
841  // documentation for this function.
842  subdomain_id_type computed_id = elemset_id + (elem_types_map[elem->type()] * n_elemsets);
843 
844  // Assign this ID to the element in question
845  elem->subdomain_id() = computed_id;
846  }
847  }
848  }
849  } // assign_subdomain_ids()
850 
851 
852 
853 
855  {
856  // Get a reference to the mesh we are reading
857  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
858 
859  // Iterate over the container of nodesets
860  container_t::iterator it=_nodeset_ids.begin();
861  for (unsigned current_id=0; it != _nodeset_ids.end(); ++it, ++current_id)
862  {
863  libMesh::out << "Assigning node boundary ID " << current_id << " to nodeset '"
864  << (*it).first
865  << "'." << std::endl;
866 
867  // Get a reference to the current vector of nodeset ID values
868  std::vector<dof_id_type>& nodeset_ids = (*it).second;
869 
870  for (std::size_t i=0; i<nodeset_ids.size(); ++i)
871  {
872  // Map the Abaqus global node ID to the libmesh node ID
873  dof_id_type libmesh_global_node_id = _abaqus_to_libmesh_node_mapping[nodeset_ids[i]];
874 
875  // Get node pointer from the mesh
876  Node* node = the_mesh.node_ptr(libmesh_global_node_id);
877 
878  if (node == NULL)
879  {
880  libMesh::err << "Error! Mesh returned NULL node pointer!" << std::endl;
881  libmesh_error();
882  }
883 
884  // Add this node with the current_id (which is determined by the
885  // alphabetical ordering of the map) to the BoundaryInfo object
886  the_mesh.boundary_info->add_node(node, current_id);
887  }
888  }
889 
890  } // assign_boundary_node_ids()
891 
892 
893 
894 
896  {
897  // Get a reference to the mesh we are reading
898  MeshBase& the_mesh = MeshInput<MeshBase>::mesh();
899 
900  // initialize the eletypes map (eletypes is a file-global variable)
901  init_eletypes();
902 
903  // Iterate over the container of sidesets
904  sideset_container_t::iterator it=_sideset_ids.begin();
905  for (unsigned current_id=0; it != _sideset_ids.end(); ++it, ++current_id)
906  {
907  libMesh::out << "Assigning sideset ID " << current_id << " to sideset '"
908  << (*it).first
909  << "'." << std::endl;
910 
911  // Get a reference to the current vector of nodeset ID values
912  std::vector<std::pair<dof_id_type,unsigned> >& sideset_ids = (*it).second;
913 
914  for (std::size_t i=0; i<sideset_ids.size(); ++i)
915  {
916  // sideset_ids is a vector of pairs (elem id, side id). Pull them out
917  // now to make the code below more readable.
918  dof_id_type abaqus_elem_id = sideset_ids[i].first;
919  unsigned abaqus_side_number = sideset_ids[i].second;
920 
921  // Map the Abaqus element ID to LibMesh numbering
922  dof_id_type libmesh_elem_id = _abaqus_to_libmesh_elem_mapping[ abaqus_elem_id ];
923 
924  // Get pointer to that element
925  Elem* elem = the_mesh.elem(libmesh_elem_id);
926 
927  // Check that the pointer returned from the Mesh is non-NULL
928  if (elem == NULL)
929  {
930  libMesh::err << "Mesh returned NULL pointer for Elem " << libmesh_elem_id << std::endl;
931  libmesh_error();
932  }
933 
934  // Grab a reference to the element definition for this element type
935  const ElementDefinition& eledef = eletypes[elem->type()];
936 
937  // If the element definition was not found, the call above would have
938  // created one with an uninitialized struct. Check for that here...
939  if (eledef.abaqus_zero_based_side_id_to_libmesh_side_id.size() == 0)
940  {
941  libMesh::err << "No Abaqus->LibMesh mapping information for ElemType " << Utility::enum_to_string(elem->type()) << "!" << std::endl;
942  libmesh_error();
943  }
944 
945  // Add this node with the current_id (which is determined by the
946  // alphabetical ordering of the map). Side numbers in Abaqus are 1-based,
947  // so we subtract 1 here before passing the abaqus side number to the
948  // mapping array
949  the_mesh.boundary_info->add_side(elem,
950  eledef.abaqus_zero_based_side_id_to_libmesh_side_id[abaqus_side_number-1],
951  current_id);
952  }
953  }
954  } // assign_sideset_ids()
955 
956 
957 
959  {
960  std::string dummy;
961  while (true)
962  {
963  // We assume we are at the beginning of a line that may be
964  // comments or may be data. We need to only discard the line if
965  // it begins with **, but we must avoid calling std::getline()
966  // since there's no way to put that back.
967  if (_in.peek() == '*')
968  {
969  // The first character was a star, so actually read it from the stream.
970  _in.get();
971 
972  // Peek at the next character...
973  if (_in.peek() == '*')
974  {
975  // OK, second character was star also, by definition this
976  // line must be a comment! Read the rest of the line and discard!
977  std::getline(_in, dummy);
978 
979  // Debugging:
980  // libMesh::out << "Read comment line: " << dummy << std::endl;
981  }
982  else
983  {
984  // The second character was _not_ a star, so put back the first star
985  // we pulled out so that the line can be parsed correctly by somebody
986  // else!
987  _in.unget();
988 
989  // Finally, break out of the while loop, we are done parsing comments
990  break;
991  }
992  }
993  else
994  {
995  // First character was not *, so this line must be data! Break out of the
996  // while loop!
997  break;
998  }
999  }
1000 
1001  } // process_and_discard_comments()
1002 
1003 
1004 } // namespace

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo