libMesh
unv_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 
19 // Local includes
20 #include "libmesh/libmesh_config.h"
21 #include "libmesh/libmesh_logging.h"
22 #include "libmesh/unv_io.h"
23 #include "libmesh/mesh_base.h"
24 #include "libmesh/edge_edge2.h"
25 #include "libmesh/face_quad4.h"
26 #include "libmesh/face_tri3.h"
27 #include "libmesh/face_tri6.h"
28 #include "libmesh/face_quad8.h"
29 #include "libmesh/face_quad9.h"
30 #include "libmesh/cell_tet4.h"
31 #include "libmesh/cell_hex8.h"
32 #include "libmesh/cell_hex20.h"
33 #include "libmesh/cell_tet10.h"
34 #include "libmesh/cell_prism6.h"
35 
36 // C++ includes
37 #include <iomanip>
38 #include <algorithm> // for std::sort
39 #include <fstream>
40 #include <ctype.h> // isspace
41 #include <sstream> // std::istringstream
42 #include <unordered_map>
43 
44 #ifdef LIBMESH_HAVE_GZSTREAM
45 # include "gzstream.h" // For reading/writing compressed streams
46 #endif
47 
48 
49 
50 namespace libMesh
51 {
52 
53 //-----------------------------------------------------------------------------
54 // UNVIO class static members
55 const std::string UNVIO::_nodes_dataset_label = "2411";
56 const std::string UNVIO::_elements_dataset_label = "2412";
57 const std::string UNVIO::_groups_dataset_label = "2467";
58 
59 
60 
61 // ------------------------------------------------------------
62 // UNVIO class members
63 
65  MeshInput<MeshBase> (mesh),
66  MeshOutput<MeshBase>(mesh),
67  _verbose (false)
68 {
69 }
70 
71 
72 
74  MeshOutput<MeshBase> (mesh),
75  _verbose (false)
76 {
77 }
78 
79 
80 
82 {
83 }
84 
85 
86 
87 bool & UNVIO::verbose ()
88 {
89  return _verbose;
90 }
91 
92 
93 
94 void UNVIO::read (const std::string & file_name)
95 {
96  if (file_name.rfind(".gz") < file_name.size())
97  {
98 #ifdef LIBMESH_HAVE_GZSTREAM
99 
100  igzstream in_stream (file_name.c_str());
101  this->read_implementation (in_stream);
102 
103 #else
104 
105  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
106 
107 #endif
108  return;
109  }
110 
111  else
112  {
113  std::ifstream in_stream (file_name.c_str());
114  this->read_implementation (in_stream);
115  return;
116  }
117 }
118 
119 
120 void UNVIO::read_implementation (std::istream & in_stream)
121 {
122  // Keep track of what kinds of elements this file contains
123  elems_of_dimension.clear();
124  elems_of_dimension.resize(4, false);
125 
126  {
127  if (!in_stream.good())
128  libmesh_error_msg("ERROR: Input file not good.");
129 
130  // Flags to be set when certain sections are encountered
131  bool
132  found_node = false,
133  found_elem = false,
134  found_group = false;
135 
136  // strings for reading the file line by line
137  std::string
138  old_line,
139  current_line;
140 
141  while (true)
142  {
143  // Save off the old_line. This will provide extra reliability
144  // for detecting the beginnings of the different sections.
145  old_line = current_line;
146 
147  // Try to read something. This may set EOF!
148  std::getline(in_stream, current_line);
149 
150  // If the stream is still "valid", parse the line
151  if (in_stream)
152  {
153  // UNV files always have some amount of leading
154  // whitespace, let's not rely on exactly how much... This
155  // command deletes it.
156  current_line.erase(std::remove_if(current_line.begin(), current_line.end(), isspace),
157  current_line.end());
158 
159  // Parse the nodes section
160  if (current_line == _nodes_dataset_label &&
161  old_line == "-1")
162  {
163  found_node = true;
164  this->nodes_in(in_stream);
165  }
166 
167  // Parse the elements section
168  else if (current_line == _elements_dataset_label &&
169  old_line == "-1")
170  {
171  // The current implementation requires the nodes to
172  // have been read before reaching the elements
173  // section.
174  if (!found_node)
175  libmesh_error_msg("ERROR: The Nodes section must come before the Elements section of the UNV file!");
176 
177  found_elem = true;
178  this->elements_in(in_stream);
179  }
180 
181  // Parse the groups section
182  else if (current_line == _groups_dataset_label &&
183  old_line == "-1")
184  {
185  // The current implementation requires the nodes and
186  // elements to have already been read before reaching
187  // the groups section.
188  if (!found_node || !found_elem)
189  libmesh_error_msg("ERROR: The Nodes and Elements sections must come before the Groups section of the UNV file!");
190 
191  found_group = true;
192  this->groups_in(in_stream);
193  }
194 
195  // We can stop reading once we've found the nodes, elements,
196  // and group sections.
197  if (found_node && found_elem && found_group)
198  break;
199 
200  // If we made it here, we're not done yet, so keep reading
201  continue;
202  }
203 
204  // if (!in_stream) check to see if EOF was set. If so, break out of while loop.
205  if (in_stream.eof())
206  break;
207 
208  // If !in_stream and !in_stream.eof(), stream is in a bad state!
209  libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
210  } // end while (true)
211 
212  // By now we better have found the datasets for nodes and elements,
213  // otherwise the unv file is bad!
214  if (!found_node)
215  libmesh_error_msg("ERROR: Could not find nodes!");
216 
217  if (!found_elem)
218  libmesh_error_msg("ERROR: Could not find elements!");
219  }
220 
221 
222 
223  {
224  // Set the mesh dimension to the largest element dimension encountered
225  MeshInput<MeshBase>::mesh().set_mesh_dimension(max_elem_dimension_seen());
226 
227 #if LIBMESH_DIM < 3
228  if (MeshInput<MeshBase>::mesh().mesh_dimension() > LIBMESH_DIM)
229  libmesh_error_msg("Cannot open dimension " \
230  << MeshInput<MeshBase>::mesh().mesh_dimension() \
231  << " mesh file when configured without " \
232  << MeshInput<MeshBase>::mesh().mesh_dimension() \
233  << "D support." );
234 #endif
235 
236  // Delete any lower-dimensional elements that might have been
237  // added to the mesh strictly for setting BCs.
238  {
239  // Grab reference to the Mesh
241 
242  unsigned char max_dim = this->max_elem_dimension_seen();
243 
244  for (const auto & elem : mesh.element_ptr_range())
245  if (elem->dim() < max_dim)
246  mesh.delete_elem(elem);
247  }
248 
249  if (this->verbose())
250  libMesh::out << " Finished." << std::endl << std::endl;
251  }
252 }
253 
254 
255 
256 
257 
258 void UNVIO::write (const std::string & file_name)
259 {
260  if (file_name.rfind(".gz") < file_name.size())
261  {
262 #ifdef LIBMESH_HAVE_GZSTREAM
263 
264  ogzstream out_stream(file_name.c_str());
265  this->write_implementation (out_stream);
266 
267 #else
268 
269  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
270 
271 #endif
272 
273  return;
274  }
275 
276  else
277  {
278  std::ofstream out_stream (file_name.c_str());
279  this->write_implementation (out_stream);
280  return;
281  }
282 }
283 
284 
285 
286 
287 void UNVIO::write_implementation (std::ostream & out_file)
288 {
289  if (!out_file.good())
290  libmesh_error_msg("ERROR: Output file not good.");
291 
292  // write the nodes, then the elements
293  this->nodes_out (out_file);
294  this->elements_out (out_file);
295 }
296 
297 
298 
299 
300 void UNVIO::nodes_in (std::istream & in_file)
301 {
302  LOG_SCOPE("nodes_in()","UNVIO");
303 
304  if (this->verbose())
305  libMesh::out << " Reading nodes" << std::endl;
306 
308 
309  // node label, we use an int here so we can read in a -1
310  int node_label;
311 
312  // always 3 coordinates in the UNV file, no matter
313  // what LIBMESH_DIM is.
314  Point xyz;
315 
316  // We'll just read the floating point values as strings even when
317  // there are no "D" characters in the file. This will make parsing
318  // the file a bit slower but the logic to do so will all be
319  // simplified and in one place...
320  std::string line;
321  std::istringstream coords_stream;
322 
323  // Continue reading nodes until there are none left
324  unsigned ctr = 0;
325  while (true)
326  {
327  // Read the node label
328  in_file >> node_label;
329 
330  // Break out of the while loop when we hit -1
331  if (node_label == -1)
332  break;
333 
334  // Read and discard the the rest of the node data on this line
335  // which we do not currently use:
336  // .) exp_coord_sys_num
337  // .) disp_coord_sys_num
338  // .) color
339  std::getline(in_file, line);
340 
341  // read the floating-point data
342  std::getline(in_file, line);
343 
344  // Replace any "D" characters used for exponents with "E"
345  size_t last_pos = 0;
346  while (true)
347  {
348  last_pos = line.find("D", last_pos);
349 
350  if (last_pos != std::string::npos)
351  line.replace(last_pos, 1, "E");
352  else
353  break;
354  }
355 
356  // Construct a stream object from the current line and then
357  // stream in the xyz values.
358  coords_stream.str(line);
359  coords_stream.clear(); // clear iostate bits! Important!
360  coords_stream >> xyz(0) >> xyz(1) >> xyz(2);
361 
362  // Add node to the Mesh, bump the counter.
363  Node * added_node = mesh.add_point(xyz, ctr++);
364 
365  // Maintain the mapping between UNV node ids and libmesh Node
366  // pointers.
367  _unv_node_id_to_libmesh_node_ptr[node_label] = added_node;
368  }
369 }
370 
371 
372 
374 {
375  unsigned char max_dim = 0;
376 
377  unsigned char elem_dimensions_size = cast_int<unsigned char>
378  (elems_of_dimension.size());
379  // The elems_of_dimension array is 1-based in the UNV reader
380  for (unsigned char i=1; i<elem_dimensions_size; ++i)
381  if (elems_of_dimension[i])
382  max_dim = i;
383 
384  return max_dim;
385 }
386 
387 
388 
389 bool UNVIO::need_D_to_e (std::string & number)
390 {
391  // find "D" in string, start looking at 6th element since dont expect a "D" earlier.
392  std::string::size_type position = number.find("D", 6);
393 
394  if (position != std::string::npos)
395  {
396  // replace "D" in string
397  number.replace(position, 1, "e");
398  return true;
399  }
400  else
401  return false;
402 }
403 
404 
405 
406 void UNVIO::groups_in (std::istream & in_file)
407 {
408  // Grab reference to the Mesh, so we can add boundary info data to it
410 
411  // Record the max and min element dimension seen while reading the file.
412  unsigned char max_dim = this->max_elem_dimension_seen();
413 
414  // Container which stores lower-dimensional elements (based on
415  // Elem::key()) for later assignment of boundary conditions. We
416  // could use e.g. an unordered_set with Elems as keys for this, but
417  // this turns out to be much slower on the search side, since we
418  // have to build an entire side in order to search, rather than just
419  // calling elem->key(side) to compute a key.
420  typedef std::unordered_multimap<dof_id_type, Elem *> map_type;
421  map_type provide_bcs;
422 
423  // Read groups until there aren't any more to read...
424  while (true)
425  {
426  // If we read a -1, it means there is nothing else to read in this section.
427  int group_number;
428  in_file >> group_number;
429 
430  if (group_number == -1)
431  break;
432 
433  // The first record consists of 8 entries:
434  // Field 1 -- group number (that we just read)
435  // Field 2 -- active constraint set no. for group
436  // Field 3 -- active restraint set no. for group
437  // Field 4 -- active load set no. for group
438  // Field 5 -- active dof set no. for group
439  // Field 6 -- active temperature set no. for group
440  // Field 7 -- active contact set no. for group
441  // Field 8 -- number of entities in group
442  // Only the first and last of these are relevant to us...
443  unsigned num_entities;
444  std::string group_name;
445  {
446  unsigned dummy;
447  in_file >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy
448  >> num_entities;
449 
450  // The second record has 1 field, the group name
451  in_file >> group_name;
452  }
453 
454  // The dimension of the elements in the group will determine
455  // whether this is a sideset group or a subdomain group.
456  bool
457  is_subdomain_group = false,
458  is_sideset_group = false;
459 
460  // Each subsequent line has 8 entries, there are two entity type
461  // codes and two tags per line that we need to read. In all my
462  // examples, the entity type code was always "8", which stands for
463  // "finite elements" but it's possible that we could eventually
464  // support other entity type codes...
465  // Field 1 -- entity type code
466  // Field 2 -- entity tag
467  // Field 3 -- entity node leaf id.
468  // Field 4 -- entity component/ ham id.
469  // Field 5 -- entity type code
470  // Field 6 -- entity tag
471  // Field 7 -- entity node leaf id.
472  // Field 8 -- entity component/ ham id.
473  {
474  unsigned entity_type_code, entity_tag, dummy;
475  for (unsigned entity=0; entity<num_entities; ++entity)
476  {
477  in_file >> entity_type_code >> entity_tag >> dummy >> dummy;
478 
479  if (entity_type_code != 8)
480  libMesh::err << "Warning, unrecognized entity type code = "
481  << entity_type_code
482  << " in group "
483  << group_name
484  << std::endl;
485 
486  // Try to find this ID in the map from UNV element ids to libmesh ids.
487  std::map<unsigned, unsigned>::iterator it =
488  _unv_elem_id_to_libmesh_elem_id.find(entity_tag);
489 
490  if (it != _unv_elem_id_to_libmesh_elem_id.end())
491  {
492  unsigned libmesh_elem_id = it->second;
493 
494  // Attempt to get a pointer to the elem listed in the group
495  Elem * group_elem = mesh.elem_ptr(libmesh_elem_id);
496 
497  // dim < max_dim means the Elem defines a boundary condition
498  if (group_elem->dim() < max_dim)
499  {
500  is_sideset_group = true;
501 
502  // We can only handle elements that are *exactly*
503  // one dimension lower than the max element
504  // dimension. Not sure if "edge" BCs in 3D
505  // actually make sense/are required...
506  if (group_elem->dim()+1 != max_dim)
507  libmesh_error_msg("ERROR: Expected boundary element of dimension " \
508  << max_dim-1 << " but got " << group_elem->dim());
509 
510  // Set the current group number as the lower-dimensional element's subdomain ID.
511  // We will use this later to set a boundary ID.
512  group_elem->subdomain_id() =
513  cast_int<subdomain_id_type>(group_number);
514 
515  // Store the lower-dimensional element in the provide_bcs container.
516  provide_bcs.insert(std::make_pair(group_elem->key(), group_elem));
517  }
518 
519  // dim == max_dim means this group defines a subdomain ID
520  else if (group_elem->dim() == max_dim)
521  {
522  is_subdomain_group = true;
523  group_elem->subdomain_id() =
524  cast_int<subdomain_id_type>(group_number);
525  }
526 
527  else
528  libmesh_error_msg("ERROR: Found an elem with dim=" \
529  << group_elem->dim() << " > " << "max_dim=" << +max_dim);
530  }
531  else
532  libMesh::err << "WARNING: UNV Element " << entity_tag << " not found while parsing groups." << std::endl;
533  } // end for (entity)
534  } // end scope
535 
536  // Associate this group_number with the group_name in the BoundaryInfo object.
537  if (is_sideset_group)
539  (cast_int<boundary_id_type>(group_number)) = group_name;
540 
541  if (is_subdomain_group)
542  mesh.subdomain_name
543  (cast_int<subdomain_id_type>(group_number)) = group_name;
544 
545  } // end while (true)
546 
547  // Loop over elements and try to assign boundary information
548  for (auto & elem : mesh.active_element_ptr_range())
549  if (elem->dim() == max_dim)
550  for (auto sn : elem->side_index_range())
551  for (const auto & pr : as_range(provide_bcs.equal_range (elem->key(sn))))
552  {
553  // This is a max-dimension element that may require BCs.
554  // For each of its sides, including internal sides, we'll
555  // see if a lower-dimensional element provides boundary
556  // information for it. Note that we have not yet called
557  // find_neighbors(), so we can't use elem->neighbor(sn) in
558  // this algorithm...
559 
560  // Build a side to confirm the hash mapped to the correct side.
561  UniquePtr<Elem> side (elem->build_side_ptr(sn));
562 
563  // Get a pointer to the lower-dimensional element
564  Elem * lower_dim_elem = pr.second;
565 
566  // This was a hash, so it might not be perfect. Let's verify...
567  if (*lower_dim_elem == *side)
568  mesh.get_boundary_info().add_side(elem, sn, lower_dim_elem->subdomain_id());
569  }
570 }
571 
572 
573 
574 void UNVIO::elements_in (std::istream & in_file)
575 {
576  LOG_SCOPE("elements_in()","UNVIO");
577 
578  if (this->verbose())
579  libMesh::out << " Reading elements" << std::endl;
580 
582 
583  // node label, we use an int here so we can read in a -1
584  int element_label;
585 
586  unsigned
587  n_nodes, // number of nodes on element
588  fe_descriptor_id, // FE descriptor id
589  phys_prop_tab_num, // physical property table number (not supported yet)
590  mat_prop_tab_num, // material property table number (not supported yet)
591  color; // color (not supported yet)
592 
593  // vector that temporarily holds the node labels defining element
594  std::vector<unsigned int> node_labels (21);
595 
596  // vector that assigns element nodes to their correct position
597  // for example:
598  // 44:plane stress | QUAD4
599  // linear quadrilateral |
600  // position in UNV-file | position in libmesh
601  // assign_elem_node[1] = 0
602  // assign_elem_node[2] = 3
603  // assign_elem_node[3] = 2
604  // assign_elem_node[4] = 1
605  //
606  // UNV is 1-based, we leave the 0th element of the vectors unused in order
607  // to prevent confusion, this way we can store elements with up to 20 nodes
608  unsigned int assign_elem_nodes[21];
609 
610  // Read elements until there are none left
611  unsigned ctr = 0;
612  while (true)
613  {
614  // read element label, break out when we read -1
615  in_file >> element_label;
616 
617  if (element_label == -1)
618  break;
619 
620  in_file >> fe_descriptor_id // read FE descriptor id
621  >> phys_prop_tab_num // (not supported yet)
622  >> mat_prop_tab_num // (not supported yet)
623  >> color // (not supported yet)
624  >> n_nodes; // read number of nodes on element
625 
626  // For "beam" type elements, the next three numbers are:
627  // .) beam orientation node number
628  // .) beam fore-end cross section number
629  // .) beam aft-end cross section number
630  // which we discard in libmesh. The "beam" type elements:
631  // 11 Rod
632  // 21 Linear beam
633  // 22 Tapered beam
634  // 23 Curved beam
635  // 24 Parabolic beam
636  // all have fe_descriptor_id < 25.
637  // http://www.sdrl.uc.edu/universal-file-formats-for-modal-analysis-testing-1/file-format-storehouse/unv_2412.htm
638  if (fe_descriptor_id < 25)
639  {
640  unsigned dummy;
641  in_file >> dummy >> dummy >> dummy;
642  }
643 
644  // read node labels (1-based)
645  for (unsigned int j=1; j<=n_nodes; j++)
646  in_file >> node_labels[j];
647 
648  // element pointer, to be allocated
649  Elem * elem = libmesh_nullptr;
650 
651  switch (fe_descriptor_id)
652  {
653  case 11: // Rod
654  {
655  elem = new Edge2;
656 
657  assign_elem_nodes[1]=0;
658  assign_elem_nodes[2]=1;
659  break;
660  }
661 
662  case 41: // Plane Stress Linear Triangle
663  case 91: // Thin Shell Linear Triangle
664  {
665  elem = new Tri3; // create new element
666 
667  assign_elem_nodes[1]=0;
668  assign_elem_nodes[2]=2;
669  assign_elem_nodes[3]=1;
670  break;
671  }
672 
673  case 42: // Plane Stress Quadratic Triangle
674  case 92: // Thin Shell Quadratic Triangle
675  {
676  elem = new Tri6; // create new element
677 
678  assign_elem_nodes[1]=0;
679  assign_elem_nodes[2]=5;
680  assign_elem_nodes[3]=2;
681  assign_elem_nodes[4]=4;
682  assign_elem_nodes[5]=1;
683  assign_elem_nodes[6]=3;
684  break;
685  }
686 
687  case 43: // Plane Stress Cubic Triangle
688  libmesh_error_msg("ERROR: UNV-element type 43: Plane Stress Cubic Triangle not supported.");
689 
690  case 44: // Plane Stress Linear Quadrilateral
691  case 94: // Thin Shell Linear Quadrilateral
692  {
693  elem = new Quad4; // create new element
694 
695  assign_elem_nodes[1]=0;
696  assign_elem_nodes[2]=3;
697  assign_elem_nodes[3]=2;
698  assign_elem_nodes[4]=1;
699  break;
700  }
701 
702  case 45: // Plane Stress Quadratic Quadrilateral
703  case 95: // Thin Shell Quadratic Quadrilateral
704  {
705  elem = new Quad8; // create new element
706 
707  assign_elem_nodes[1]=0;
708  assign_elem_nodes[2]=7;
709  assign_elem_nodes[3]=3;
710  assign_elem_nodes[4]=6;
711  assign_elem_nodes[5]=2;
712  assign_elem_nodes[6]=5;
713  assign_elem_nodes[7]=1;
714  assign_elem_nodes[8]=4;
715  break;
716  }
717 
718  case 300: // Thin Shell Quadratic Quadrilateral (nine nodes)
719  {
720  elem = new Quad9; // create new element
721 
722  assign_elem_nodes[1]=0;
723  assign_elem_nodes[2]=7;
724  assign_elem_nodes[3]=3;
725  assign_elem_nodes[4]=6;
726  assign_elem_nodes[5]=2;
727  assign_elem_nodes[6]=5;
728  assign_elem_nodes[7]=1;
729  assign_elem_nodes[8]=4;
730  assign_elem_nodes[9]=8;
731  break;
732  }
733 
734  case 46: // Plane Stress Cubic Quadrilateral
735  libmesh_error_msg("ERROR: UNV-element type 46: Plane Stress Cubic Quadrilateral not supported.");
736 
737  case 111: // Solid Linear Tetrahedron
738  {
739  elem = new Tet4; // create new element
740 
741  assign_elem_nodes[1]=0;
742  assign_elem_nodes[2]=1;
743  assign_elem_nodes[3]=2;
744  assign_elem_nodes[4]=3;
745  break;
746  }
747 
748  case 112: // Solid Linear Prism
749  {
750  elem = new Prism6; // create new element
751 
752  assign_elem_nodes[1]=0;
753  assign_elem_nodes[2]=1;
754  assign_elem_nodes[3]=2;
755  assign_elem_nodes[4]=3;
756  assign_elem_nodes[5]=4;
757  assign_elem_nodes[6]=5;
758  break;
759  }
760 
761  case 115: // Solid Linear Brick
762  {
763  elem = new Hex8; // create new element
764 
765  assign_elem_nodes[1]=0;
766  assign_elem_nodes[2]=4;
767  assign_elem_nodes[3]=5;
768  assign_elem_nodes[4]=1;
769  assign_elem_nodes[5]=3;
770  assign_elem_nodes[6]=7;
771  assign_elem_nodes[7]=6;
772  assign_elem_nodes[8]=2;
773  break;
774  }
775 
776  case 116: // Solid Quadratic Brick
777  {
778  elem = new Hex20; // create new element
779 
780  assign_elem_nodes[1]=0;
781  assign_elem_nodes[2]=12;
782  assign_elem_nodes[3]=4;
783  assign_elem_nodes[4]=16;
784  assign_elem_nodes[5]=5;
785  assign_elem_nodes[6]=13;
786  assign_elem_nodes[7]=1;
787  assign_elem_nodes[8]=8;
788 
789  assign_elem_nodes[9]=11;
790  assign_elem_nodes[10]=19;
791  assign_elem_nodes[11]=17;
792  assign_elem_nodes[12]=9;
793 
794  assign_elem_nodes[13]=3;
795  assign_elem_nodes[14]=15;
796  assign_elem_nodes[15]=7;
797  assign_elem_nodes[16]=18;
798  assign_elem_nodes[17]=6;
799  assign_elem_nodes[18]=14;
800  assign_elem_nodes[19]=2;
801  assign_elem_nodes[20]=10;
802  break;
803  }
804 
805  case 117: // Solid Cubic Brick
806  libmesh_error_msg("Error: UNV-element type 117: Solid Cubic Brick not supported.");
807 
808  case 118: // Solid Quadratic Tetrahedron
809  {
810  elem = new Tet10; // create new element
811 
812  assign_elem_nodes[1]=0;
813  assign_elem_nodes[2]=4;
814  assign_elem_nodes[3]=1;
815  assign_elem_nodes[4]=5;
816  assign_elem_nodes[5]=2;
817  assign_elem_nodes[6]=6;
818  assign_elem_nodes[7]=7;
819  assign_elem_nodes[8]=8;
820  assign_elem_nodes[9]=9;
821  assign_elem_nodes[10]=3;
822  break;
823  }
824 
825  default: // Unrecognized element type
826  libmesh_error_msg("ERROR: UNV-element type " << fe_descriptor_id << " not supported.");
827  }
828 
829  // nodes are being stored in element
830  for (dof_id_type j=1; j<=n_nodes; j++)
831  {
832  // Map the UNV node ID to the libmesh node ID
833  std::map<dof_id_type, Node *>::iterator it =
834  _unv_node_id_to_libmesh_node_ptr.find(node_labels[j]);
835 
836  if (it != _unv_node_id_to_libmesh_node_ptr.end())
837  elem->set_node(assign_elem_nodes[j]) = it->second;
838  else
839  libmesh_error_msg("ERROR: UNV node " << node_labels[j] << " not found!");
840  }
841 
842  elems_of_dimension[elem->dim()] = true;
843 
844  // Set the element's ID
845  elem->set_id(ctr);
846 
847  // Maintain a map from the libmesh (0-based) numbering to the
848  // UNV numbering.
849  _unv_elem_id_to_libmesh_elem_id[element_label] = ctr;
850 
851  // Add the element to the Mesh
852  mesh.add_elem(elem);
853 
854  // Increment the counter for the next iteration
855  ctr++;
856  } // end while (true)
857 }
858 
859 
860 
861 
862 
863 
864 void UNVIO::nodes_out (std::ostream & out_file)
865 {
866  if (this->verbose())
867  libMesh::out << " Writing " << MeshOutput<MeshBase>::mesh().n_nodes() << " nodes" << std::endl;
868 
869  // Write beginning of dataset
870  out_file << " -1\n"
871  << " "
873  << '\n';
874 
875  unsigned int
876  exp_coord_sys_dummy = 0, // export coordinate sys. (not supported yet)
877  disp_coord_sys_dummy = 0, // displacement coordinate sys. (not supp. yet)
878  color_dummy = 0; // color(not supported yet)
879 
880  // A reference to the parent class's mesh
882 
883  // Use scientific notation with capital E and 16 digits for printing out the coordinates
884  out_file << std::scientific << std::setprecision(16) << std::uppercase;
885 
886  for (const auto & current_node : mesh.node_ptr_range())
887  {
888  dof_id_type node_id = current_node->id();
889 
890  out_file << std::setw(10) << node_id
891  << std::setw(10) << exp_coord_sys_dummy
892  << std::setw(10) << disp_coord_sys_dummy
893  << std::setw(10) << color_dummy
894  << '\n';
895 
896  // The coordinates - always write out three coords regardless of LIBMESH_DIM
897  Real x = (*current_node)(0);
898 
899 #if LIBMESH_DIM > 1
900  Real y = (*current_node)(1);
901 #else
902  Real y = 0.;
903 #endif
904 
905 #if LIBMESH_DIM > 2
906  Real z = (*current_node)(2);
907 #else
908  Real z = 0.;
909 #endif
910 
911  out_file << std::setw(25) << x
912  << std::setw(25) << y
913  << std::setw(25) << z
914  << '\n';
915  }
916 
917  // Write end of dataset
918  out_file << " -1\n";
919 }
920 
921 
922 
923 
924 
925 
926 void UNVIO::elements_out(std::ostream & out_file)
927 {
928  if (this->verbose())
929  libMesh::out << " Writing elements" << std::endl;
930 
931  // Write beginning of dataset
932  out_file << " -1\n"
933  << " "
935  << "\n";
936 
937  unsigned
938  fe_descriptor_id = 0, // FE descriptor id
939  phys_prop_tab_dummy = 2, // physical property (not supported yet)
940  mat_prop_tab_dummy = 1, // material property (not supported yet)
941  color_dummy = 7; // color (not supported yet)
942 
943  // vector that assigns element nodes to their correct position
944  // currently only elements with up to 20 nodes
945  //
946  // Example:
947  // QUAD4 | 44:plane stress
948  // | linear quad
949  // position in libMesh | UNV numbering
950  // (note: 0-based) | (note: 1-based)
951  //
952  // assign_elem_node[0] = 0
953  // assign_elem_node[1] = 3
954  // assign_elem_node[2] = 2
955  // assign_elem_node[3] = 1
956  unsigned int assign_elem_nodes[20];
957 
958  unsigned int n_elem_written=0;
959 
960  // A reference to the parent class's mesh
962 
963  for (const auto & elem : mesh.element_ptr_range())
964  {
965  switch (elem->type())
966  {
967 
968  case TRI3:
969  {
970  fe_descriptor_id = 41; // Plane Stress Linear Triangle
971  assign_elem_nodes[0] = 0;
972  assign_elem_nodes[1] = 2;
973  assign_elem_nodes[2] = 1;
974  break;
975  }
976 
977  case TRI6:
978  {
979  fe_descriptor_id = 42; // Plane Stress Quadratic Triangle
980  assign_elem_nodes[0] = 0;
981  assign_elem_nodes[1] = 5;
982  assign_elem_nodes[2] = 2;
983  assign_elem_nodes[3] = 4;
984  assign_elem_nodes[4] = 1;
985  assign_elem_nodes[5] = 3;
986  break;
987  }
988 
989  case QUAD4:
990  {
991  fe_descriptor_id = 44; // Plane Stress Linear Quadrilateral
992  assign_elem_nodes[0] = 0;
993  assign_elem_nodes[1] = 3;
994  assign_elem_nodes[2] = 2;
995  assign_elem_nodes[3] = 1;
996  break;
997  }
998 
999  case QUAD8:
1000  {
1001  fe_descriptor_id = 45; // Plane Stress Quadratic Quadrilateral
1002  assign_elem_nodes[0] = 0;
1003  assign_elem_nodes[1] = 7;
1004  assign_elem_nodes[2] = 3;
1005  assign_elem_nodes[3] = 6;
1006  assign_elem_nodes[4] = 2;
1007  assign_elem_nodes[5] = 5;
1008  assign_elem_nodes[6] = 1;
1009  assign_elem_nodes[7] = 4;
1010  break;
1011  }
1012 
1013  case QUAD9:
1014  {
1015  fe_descriptor_id = 300; // Plane Stress Quadratic Quadrilateral
1016  assign_elem_nodes[0] = 0;
1017  assign_elem_nodes[1] = 7;
1018  assign_elem_nodes[2] = 3;
1019  assign_elem_nodes[3] = 6;
1020  assign_elem_nodes[4] = 2;
1021  assign_elem_nodes[5] = 5;
1022  assign_elem_nodes[6] = 1;
1023  assign_elem_nodes[7] = 4;
1024  assign_elem_nodes[8] = 8;
1025  break;
1026  }
1027 
1028  case TET4:
1029  {
1030  fe_descriptor_id = 111; // Solid Linear Tetrahedron
1031  assign_elem_nodes[0] = 0;
1032  assign_elem_nodes[1] = 1;
1033  assign_elem_nodes[2] = 2;
1034  assign_elem_nodes[3] = 3;
1035  break;
1036  }
1037 
1038  case PRISM6:
1039  {
1040  fe_descriptor_id = 112; // Solid Linear Prism
1041  assign_elem_nodes[0] = 0;
1042  assign_elem_nodes[1] = 1;
1043  assign_elem_nodes[2] = 2;
1044  assign_elem_nodes[3] = 3;
1045  assign_elem_nodes[4] = 4;
1046  assign_elem_nodes[5] = 5;
1047  break;
1048  }
1049 
1050  case HEX8:
1051  {
1052  fe_descriptor_id = 115; // Solid Linear Brick
1053  assign_elem_nodes[0] = 0;
1054  assign_elem_nodes[1] = 4;
1055  assign_elem_nodes[2] = 5;
1056  assign_elem_nodes[3] = 1;
1057  assign_elem_nodes[4] = 3;
1058  assign_elem_nodes[5] = 7;
1059  assign_elem_nodes[6] = 6;
1060  assign_elem_nodes[7] = 2;
1061  break;
1062  }
1063 
1064  case HEX20:
1065  {
1066  fe_descriptor_id = 116; // Solid Quadratic Brick
1067  assign_elem_nodes[ 0] = 0;
1068  assign_elem_nodes[ 1] = 12;
1069  assign_elem_nodes[ 2] = 4;
1070  assign_elem_nodes[ 3] = 16;
1071  assign_elem_nodes[ 4] = 5;
1072  assign_elem_nodes[ 5] = 13;
1073  assign_elem_nodes[ 6] = 1;
1074  assign_elem_nodes[ 7] = 8;
1075 
1076  assign_elem_nodes[ 8] = 11;
1077  assign_elem_nodes[ 9] = 19;
1078  assign_elem_nodes[10] = 17;
1079  assign_elem_nodes[11] = 9;
1080 
1081  assign_elem_nodes[12] = 3;
1082  assign_elem_nodes[13] = 15;
1083  assign_elem_nodes[14] = 7;
1084  assign_elem_nodes[15] = 18;
1085  assign_elem_nodes[16] = 6;
1086  assign_elem_nodes[17] = 14;
1087  assign_elem_nodes[18] = 2;
1088  assign_elem_nodes[19] = 10;
1089 
1090 
1091  break;
1092  }
1093 
1094  case TET10:
1095  {
1096  fe_descriptor_id = 118; // Solid Quadratic Tetrahedron
1097  assign_elem_nodes[0] = 0;
1098  assign_elem_nodes[1] = 4;
1099  assign_elem_nodes[2] = 1;
1100  assign_elem_nodes[3] = 5;
1101  assign_elem_nodes[4] = 2;
1102  assign_elem_nodes[5] = 6;
1103  assign_elem_nodes[6] = 7;
1104  assign_elem_nodes[7] = 8;
1105  assign_elem_nodes[8] = 9;
1106  assign_elem_nodes[9] = 3;
1107  break;
1108  }
1109 
1110  default:
1111  libmesh_error_msg("ERROR: Element type = " << elem->type() << " not supported in " << "UNVIO!");
1112  }
1113 
1114  dof_id_type elem_id = elem->id();
1115 
1116  out_file << std::setw(10) << elem_id // element ID
1117  << std::setw(10) << fe_descriptor_id // type of element
1118  << std::setw(10) << phys_prop_tab_dummy // not supported
1119  << std::setw(10) << mat_prop_tab_dummy // not supported
1120  << std::setw(10) << color_dummy // not supported
1121  << std::setw(10) << elem->n_nodes() // No. of nodes per element
1122  << '\n';
1123 
1124  for (auto j : elem->node_index_range())
1125  {
1126  // assign_elem_nodes[j]-th node: i.e., j loops over the
1127  // libMesh numbering, and assign_elem_nodes[j] over the
1128  // UNV numbering.
1129  const Node * node_in_unv_order = elem->node_ptr(assign_elem_nodes[j]);
1130 
1131  // new record after 8 id entries
1132  if (j==8 || j==16)
1133  out_file << '\n';
1134 
1135  // Write foreign label for this node
1136  dof_id_type node_id = node_in_unv_order->id();
1137 
1138  out_file << std::setw(10) << node_id;
1139  }
1140 
1141  out_file << '\n';
1142 
1143  n_elem_written++;
1144  }
1145 
1146  if (this->verbose())
1147  libMesh::out << " Finished writing " << n_elem_written << " elements" << std::endl;
1148 
1149  // Write end of dataset
1150  out_file << " -1\n";
1151 }
1152 
1153 
1154 
1155 void UNVIO::read_dataset(std::string file_name)
1156 {
1157  std::ifstream in_stream(file_name.c_str());
1158 
1159  if (!in_stream.good())
1160  libmesh_error_msg("Error opening UNV data file.");
1161 
1162  std::string olds, news, dummy;
1163 
1164  while (true)
1165  {
1166  in_stream >> olds >> news;
1167 
1168  // A "-1" followed by a number means the beginning of a dataset.
1169  while (((olds != "-1") || (news == "-1")) && !in_stream.eof())
1170  {
1171  olds = news;
1172  in_stream >> news;
1173  }
1174 
1175  if (in_stream.eof())
1176  break;
1177 
1178  // We only support reading datasets in the "2414" format.
1179  if (news == "2414")
1180  {
1181  // Ignore the rest of the current line and the next two
1182  // lines that contain analysis dataset label and name.
1183  for (unsigned int i=0; i<3; i++)
1184  std::getline(in_stream, dummy);
1185 
1186  // Read the dataset location, where
1187  // 1: Data at nodes
1188  // 2: Data on elements
1189  // Currently only data on nodes is supported.
1190  unsigned int dataset_location;
1191  in_stream >> dataset_location;
1192 
1193  // Currently only nodal datasets are supported.
1194  if (dataset_location != 1)
1195  libmesh_error_msg("ERROR: Currently only Data at nodes is supported.");
1196 
1197  // Ignore the rest of this line and the next five records.
1198  for (unsigned int i=0; i<6; i++)
1199  std::getline(in_stream, dummy);
1200 
1201  // These data are all of no interest to us...
1202  unsigned int
1203  model_type,
1204  analysis_type,
1205  data_characteristic,
1206  result_type;
1207 
1208  // The type of data (complex, real, float, double etc, see
1209  // below).
1210  unsigned int data_type;
1211 
1212  // The number of floating-point values per entity.
1213  unsigned int num_vals_per_node;
1214 
1215  in_stream >> model_type // not used here
1216  >> analysis_type // not used here
1217  >> data_characteristic // not used here
1218  >> result_type // not used here
1219  >> data_type
1220  >> num_vals_per_node;
1221 
1222  // Ignore the rest of record 9, and records 10-13.
1223  for (unsigned int i=0; i<5; i++)
1224  std::getline(in_stream, dummy);
1225 
1226  // Now get the foreign (aka UNV node) node number and
1227  // the respective nodal data.
1228  int f_n_id;
1229  std::vector<Number> values;
1230 
1231  while (true)
1232  {
1233  in_stream >> f_n_id;
1234 
1235  // If -1 then we have reached the end of the dataset.
1236  if (f_n_id == -1)
1237  break;
1238 
1239  // Resize the values vector (usually data in three
1240  // principle directions, i.e. num_vals_per_node = 3).
1241  values.resize(num_vals_per_node);
1242 
1243  // Read the values for the respective node.
1244  for (unsigned int data_cnt=0; data_cnt<num_vals_per_node; data_cnt++)
1245  {
1246  // Check what data type we are reading.
1247  // 2,4: Real
1248  // 5,6: Complex
1249  // other data types are not supported yet.
1250  // As again, these floats may also be written
1251  // using a 'D' instead of an 'e'.
1252  if (data_type == 2 || data_type == 4)
1253  {
1254  std::string buf;
1255  in_stream >> buf;
1256  this->need_D_to_e(buf);
1257 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1258  values[data_cnt] = Complex(std::atof(buf.c_str()), 0.);
1259 #else
1260  values[data_cnt] = std::atof(buf.c_str());
1261 #endif
1262  }
1263 
1264  else if (data_type == 5 || data_type == 6)
1265  {
1266 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1267  Real re_val, im_val;
1268 
1269  std::string buf;
1270  in_stream >> buf;
1271 
1272  if (this->need_D_to_e(buf))
1273  {
1274  re_val = std::atof(buf.c_str());
1275  in_stream >> buf;
1276  this->need_D_to_e(buf);
1277  im_val = std::atof(buf.c_str());
1278  }
1279  else
1280  {
1281  re_val = std::atof(buf.c_str());
1282  in_stream >> im_val;
1283  }
1284 
1285  values[data_cnt] = Complex(re_val,im_val);
1286 #else
1287 
1288  libmesh_error_msg("ERROR: Complex data only supported when libMesh is configured with --enable-complex!");
1289 #endif
1290  }
1291 
1292  else
1293  libmesh_error_msg("ERROR: Data type not supported.");
1294 
1295  } // end loop data_cnt
1296 
1297  // Get a pointer to the Node associated with the UNV node id.
1298  std::map<dof_id_type, Node *>::const_iterator it =
1299  _unv_node_id_to_libmesh_node_ptr.find(f_n_id);
1300 
1301  if (it == _unv_node_id_to_libmesh_node_ptr.end())
1302  libmesh_error_msg("UNV node id " << f_n_id << " was not found.");
1303 
1304  // Store the nodal values in our map which uses the
1305  // libMesh Node* as the key. We use operator[] here
1306  // because we want to create an empty vector if the
1307  // entry does not already exist.
1308  _node_data[it->second] = values;
1309  } // end while (true)
1310  } // end if (news == "2414")
1311  } // end while (true)
1312 }
1313 
1314 
1315 
1316 const std::vector<Number> *
1317 UNVIO::get_data (Node * node) const
1318 {
1319  std::map<Node *, std::vector<Number>>::const_iterator
1320  it = _node_data.find(node);
1321 
1322  if (it == _node_data.end())
1323  return libmesh_nullptr;
1324  else
1325  return &(it->second);
1326 }
1327 
1328 
1329 } // namespace libMesh
OStreamProxy err
void write_implementation(std::ostream &out_stream)
The actual implementation of the write function.
Definition: unv_io.C:287
bool _verbose
should be be verbose?
Definition: unv_io.h:182
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:54
UNVIO(MeshBase &mesh)
Constructor.
Definition: unv_io.C:64
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
const std::vector< Number > * get_data(Node *node) const
Definition: unv_io.C:1317
bool & verbose()
Set the flag indicating if we should be verbose.
Definition: unv_io.C:87
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:49
std::map< dof_id_type, Node * > _unv_node_id_to_libmesh_node_ptr
Maps UNV node IDs to libMesh Node*s.
Definition: unv_io.h:188
static const std::string _elements_dataset_label
label for the element dataset
Definition: unv_io.h:198
virtual void read(const std::string &) libmesh_override
This method implements reading a mesh from a specified file.
Definition: unv_io.C:94
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
void nodes_in(std::istream &in_file)
Read nodes from file.
Definition: unv_io.C:300
void read_implementation(std::istream &in_stream)
The actual implementation of the read function.
Definition: unv_io.C:120
void read_dataset(std::string file_name)
Read a UNV data file containing a dataset of type "2414".
Definition: unv_io.C:1155
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
std::map< Node *, std::vector< Number > > _node_data
Map from libMesh Node* to data at that node, as read in by the read_dataset() function.
Definition: unv_io.h:214
const class libmesh_nullptr_t libmesh_nullptr
The QUAD8 is an element in 2D composed of 8 nodes.
Definition: face_quad8.h:49
virtual void write(const std::string &) libmesh_override
This method implements writing a mesh to a specified file.
Definition: unv_io.C:258
MPI_Datatype data_type
Data types for communication.
Definition: parallel.h:166
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
const MT & mesh() const
Definition: mesh_output.h:216
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
unsigned char max_elem_dimension_seen()
Definition: unv_io.C:373
The libMesh namespace provides an interface to certain functionality in the library.
void nodes_out(std::ostream &out_file)
Outputs nodes to the file out_file.
Definition: unv_io.C:864
static const std::string _groups_dataset_label
label for the groups dataset
Definition: unv_io.h:203
void elements_out(std::ostream &out_file)
Outputs the element data to the file out_file.
Definition: unv_io.C:926
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.
This is the MeshBase class.
Definition: mesh_base.h:68
dof_id_type & set_id()
Definition: dof_object.h:641
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
The Tri6 is an element in 2D composed of 6 nodes.
Definition: face_tri6.h:49
const dof_id_type n_nodes
Definition: tecplot_io.C:67
PetscErrorCode Vec x
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
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual SimpleRange< node_iterator > node_ptr_range()=0
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
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:576
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
The QUAD9 is an element in 2D composed of 9 nodes.
Definition: face_quad9.h:49
std::map< unsigned, unsigned > _unv_elem_id_to_libmesh_elem_id
Map UNV element IDs to libmesh element IDs.
Definition: unv_io.h:208
std::string & sideset_name(boundary_id_type id)
std::complex< Real > Complex
bool need_D_to_e(std::string &number)
Replaces "1.1111D+00" with "1.1111e+00" if necessary.
Definition: unv_io.C:389
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
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
The Edge2 is an element in 1D composed of 2 nodes.
Definition: edge_edge2.h:43
virtual ~UNVIO()
Destructor.
Definition: unv_io.C:81
dof_id_type id() const
Definition: dof_object.h:632
static const std::string _nodes_dataset_label
label for the node dataset
Definition: unv_io.h:193
virtual const Elem * elem_ptr(const dof_id_type i) const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void elements_in(std::istream &in_file)
Method reads elements and stores them in std::vector<Elem *> _elements in the same order as they come...
Definition: unv_io.C:574
uint8_t dof_id_type
Definition: id_types.h:64
void groups_in(std::istream &in_file)
Reads the "groups" section of the file.
Definition: unv_io.C:406