unv_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 
19 // Local includes
20 #include "libmesh/libmesh_config.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 "libmesh/ignore_warnings.h" // shadowing in gzstream.h
46 # include "gzstream.h" // For reading/writing compressed streams
48 #endif
49 
50 
51 
52 namespace libMesh
53 {
54 
55 //-----------------------------------------------------------------------------
56 // UNVIO class static members
57 const std::string UNVIO::_nodes_dataset_label = "2411";
58 const std::string UNVIO::_elements_dataset_label = "2412";
59 const std::string UNVIO::_groups_dataset_label = "2467";
60 
61 
62 
63 // ------------------------------------------------------------
64 // UNVIO class members
65 
69  _verbose (false)
70 {
71 }
72 
73 
74 
77  _verbose (false)
78 {
79 }
80 
81 
82 
84 {
85 }
86 
87 
88 
89 bool & UNVIO::verbose ()
90 {
91  return _verbose;
92 }
93 
94 
95 
96 void UNVIO::read (const std::string & file_name)
97 {
98  if (file_name.rfind(".gz") < file_name.size())
99  {
100 #ifdef LIBMESH_HAVE_GZSTREAM
101 
102  igzstream in_stream (file_name.c_str());
103  this->read_implementation (in_stream);
104 
105 #else
106 
107  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
108 
109 #endif
110  return;
111  }
112 
113  else
114  {
115  std::ifstream in_stream (file_name.c_str());
116  this->read_implementation (in_stream);
117  return;
118  }
119 }
120 
121 
122 void UNVIO::read_implementation (std::istream & in_stream)
123 {
124  // Keep track of what kinds of elements this file contains
125  elems_of_dimension.clear();
126  elems_of_dimension.resize(4, false);
127 
128  {
129  if (!in_stream.good())
130  libmesh_error_msg("ERROR: Input file not good.");
131 
132  // Flags to be set when certain sections are encountered
133  bool
134  found_node = false,
135  found_elem = false,
136  found_group = false;
137 
138  // strings for reading the file line by line
139  std::string
140  old_line,
141  current_line;
142 
143  while (true)
144  {
145  // Save off the old_line. This will provide extra reliability
146  // for detecting the beginnings of the different sections.
147  old_line = current_line;
148 
149  // Try to read something. This may set EOF!
150  std::getline(in_stream, current_line);
151 
152  // If the stream is still "valid", parse the line
153  if (in_stream)
154  {
155  // UNV files always have some amount of leading
156  // whitespace, let's not rely on exactly how much... This
157  // command deletes it.
158  current_line.erase(std::remove_if(current_line.begin(), current_line.end(), isspace),
159  current_line.end());
160 
161  // Parse the nodes section
162  if (current_line == _nodes_dataset_label &&
163  old_line == "-1")
164  {
165  found_node = true;
166  this->nodes_in(in_stream);
167  }
168 
169  // Parse the elements section
170  else if (current_line == _elements_dataset_label &&
171  old_line == "-1")
172  {
173  // The current implementation requires the nodes to
174  // have been read before reaching the elements
175  // section.
176  if (!found_node)
177  libmesh_error_msg("ERROR: The Nodes section must come before the Elements section of the UNV file!");
178 
179  found_elem = true;
180  this->elements_in(in_stream);
181  }
182 
183  // Parse the groups section
184  else if (current_line == _groups_dataset_label &&
185  old_line == "-1")
186  {
187  // The current implementation requires the nodes and
188  // elements to have already been read before reaching
189  // the groups section.
190  if (!found_node || !found_elem)
191  libmesh_error_msg("ERROR: The Nodes and Elements sections must come before the Groups section of the UNV file!");
192 
193  found_group = true;
194  this->groups_in(in_stream);
195  }
196 
197  // We can stop reading once we've found the nodes, elements,
198  // and group sections.
199  if (found_node && found_elem && found_group)
200  break;
201 
202  // If we made it here, we're not done yet, so keep reading
203  continue;
204  }
205 
206  // if (!in_stream) check to see if EOF was set. If so, break out of while loop.
207  if (in_stream.eof())
208  break;
209 
210  // If !in_stream and !in_stream.eof(), stream is in a bad state!
211  libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
212  } // end while (true)
213 
214  // By now we better have found the datasets for nodes and elements,
215  // otherwise the unv file is bad!
216  if (!found_node)
217  libmesh_error_msg("ERROR: Could not find nodes!");
218 
219  if (!found_elem)
220  libmesh_error_msg("ERROR: Could not find elements!");
221  }
222 
223 
224 
225  {
226  // Set the mesh dimension to the largest element dimension encountered
227  MeshInput<MeshBase>::mesh().set_mesh_dimension(max_elem_dimension_seen());
228 
229 #if LIBMESH_DIM < 3
230  if (MeshInput<MeshBase>::mesh().mesh_dimension() > LIBMESH_DIM)
231  libmesh_error_msg("Cannot open dimension " \
232  << MeshInput<MeshBase>::mesh().mesh_dimension() \
233  << " mesh file when configured without " \
234  << MeshInput<MeshBase>::mesh().mesh_dimension() \
235  << "D support." );
236 #endif
237 
238  // Delete any lower-dimensional elements that might have been
239  // added to the mesh strictly for setting BCs.
240  {
241  // Grab reference to the Mesh
243 
244  unsigned char max_dim = this->max_elem_dimension_seen();
245 
246  for (const auto & elem : mesh.element_ptr_range())
247  if (elem->dim() < max_dim)
248  mesh.delete_elem(elem);
249  }
250 
251  if (this->verbose())
252  libMesh::out << " Finished." << std::endl << std::endl;
253  }
254 }
255 
256 
257 
258 
259 
260 void UNVIO::write (const std::string & file_name)
261 {
262  if (file_name.rfind(".gz") < file_name.size())
263  {
264 #ifdef LIBMESH_HAVE_GZSTREAM
265 
266  ogzstream out_stream(file_name.c_str());
267  this->write_implementation (out_stream);
268 
269 #else
270 
271  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
272 
273 #endif
274 
275  return;
276  }
277 
278  else
279  {
280  std::ofstream out_stream (file_name.c_str());
281  this->write_implementation (out_stream);
282  return;
283  }
284 }
285 
286 
287 
288 
289 void UNVIO::write_implementation (std::ostream & out_file)
290 {
291  if (!out_file.good())
292  libmesh_error_msg("ERROR: Output file not good.");
293 
294  // write the nodes, then the elements
295  this->nodes_out (out_file);
296  this->elements_out (out_file);
297 }
298 
299 
300 
301 
302 void UNVIO::nodes_in (std::istream & in_file)
303 {
304  LOG_SCOPE("nodes_in()","UNVIO");
305 
306  if (this->verbose())
307  libMesh::out << " Reading nodes" << std::endl;
308 
310 
311  // node label, we use an int here so we can read in a -1
312  int node_label;
313 
314  // always 3 coordinates in the UNV file, no matter
315  // what LIBMESH_DIM is.
316  Point xyz;
317 
318  // We'll just read the floating point values as strings even when
319  // there are no "D" characters in the file. This will make parsing
320  // the file a bit slower but the logic to do so will all be
321  // simplified and in one place...
322  std::string line;
323  std::istringstream coords_stream;
324 
325  // Continue reading nodes until there are none left
326  unsigned ctr = 0;
327  while (true)
328  {
329  // Read the node label
330  in_file >> node_label;
331 
332  // Break out of the while loop when we hit -1
333  if (node_label == -1)
334  break;
335 
336  // Read and discard the the rest of the node data on this line
337  // which we do not currently use:
338  // .) exp_coord_sys_num
339  // .) disp_coord_sys_num
340  // .) color
341  std::getline(in_file, line);
342 
343  // read the floating-point data
344  std::getline(in_file, line);
345 
346  // Replace any "D" characters used for exponents with "E"
347  size_t last_pos = 0;
348  while (true)
349  {
350  last_pos = line.find("D", last_pos);
351 
352  if (last_pos != std::string::npos)
353  line.replace(last_pos, 1, "E");
354  else
355  break;
356  }
357 
358  // Construct a stream object from the current line and then
359  // stream in the xyz values.
360  coords_stream.str(line);
361  coords_stream.clear(); // clear iostate bits! Important!
362  coords_stream >> xyz(0) >> xyz(1) >> xyz(2);
363 
364  // Add node to the Mesh, bump the counter.
365  Node * added_node = mesh.add_point(xyz, ctr++);
366 
367  // Maintain the mapping between UNV node ids and libmesh Node
368  // pointers.
369  _unv_node_id_to_libmesh_node_ptr[node_label] = added_node;
370  }
371 }
372 
373 
374 
376 {
377  unsigned char max_dim = 0;
378 
379  unsigned char elem_dimensions_size = cast_int<unsigned char>
380  (elems_of_dimension.size());
381  // The elems_of_dimension array is 1-based in the UNV reader
382  for (unsigned char i=1; i<elem_dimensions_size; ++i)
383  if (elems_of_dimension[i])
384  max_dim = i;
385 
386  return max_dim;
387 }
388 
389 
390 
391 bool UNVIO::need_D_to_e (std::string & number)
392 {
393  // find "D" in string, start looking at 6th element since dont expect a "D" earlier.
394  std::string::size_type position = number.find("D", 6);
395 
396  if (position != std::string::npos)
397  {
398  // replace "D" in string
399  number.replace(position, 1, "e");
400  return true;
401  }
402  else
403  return false;
404 }
405 
406 
407 
408 void UNVIO::groups_in (std::istream & in_file)
409 {
410  // Grab reference to the Mesh, so we can add boundary info data to it
412 
413  // Record the max and min element dimension seen while reading the file.
414  unsigned char max_dim = this->max_elem_dimension_seen();
415 
416  // Container which stores lower-dimensional elements (based on
417  // Elem::key()) for later assignment of boundary conditions. We
418  // could use e.g. an unordered_set with Elems as keys for this, but
419  // this turns out to be much slower on the search side, since we
420  // have to build an entire side in order to search, rather than just
421  // calling elem->key(side) to compute a key.
422  typedef std::unordered_multimap<dof_id_type, Elem *> map_type;
423  map_type provide_bcs;
424 
425  // Read groups until there aren't any more to read...
426  while (true)
427  {
428  // If we read a -1, it means there is nothing else to read in this section.
429  int group_number;
430  in_file >> group_number;
431 
432  if (group_number == -1)
433  break;
434 
435  // The first record consists of 8 entries:
436  // Field 1 -- group number (that we just read)
437  // Field 2 -- active constraint set no. for group
438  // Field 3 -- active restraint set no. for group
439  // Field 4 -- active load set no. for group
440  // Field 5 -- active dof set no. for group
441  // Field 6 -- active temperature set no. for group
442  // Field 7 -- active contact set no. for group
443  // Field 8 -- number of entities in group
444  // Only the first and last of these are relevant to us...
445  unsigned num_entities;
446  std::string group_name;
447  {
448  unsigned dummy;
449  in_file >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy
450  >> num_entities;
451 
452  // The second record has 1 field, the group name
453  in_file >> group_name;
454  }
455 
456  // The dimension of the elements in the group will determine
457  // whether this is a sideset group or a subdomain group.
458  bool
459  is_subdomain_group = false,
460  is_sideset_group = false;
461 
462  // Each subsequent line has 8 entries, there are two entity type
463  // codes and two tags per line that we need to read. In all my
464  // examples, the entity type code was always "8", which stands for
465  // "finite elements" but it's possible that we could eventually
466  // support other entity type codes...
467  // Field 1 -- entity type code
468  // Field 2 -- entity tag
469  // Field 3 -- entity node leaf id.
470  // Field 4 -- entity component/ ham id.
471  // Field 5 -- entity type code
472  // Field 6 -- entity tag
473  // Field 7 -- entity node leaf id.
474  // Field 8 -- entity component/ ham id.
475  {
476  unsigned entity_type_code, entity_tag, dummy;
477  for (unsigned entity=0; entity<num_entities; ++entity)
478  {
479  in_file >> entity_type_code >> entity_tag >> dummy >> dummy;
480 
481  if (entity_type_code != 8)
482  libMesh::err << "Warning, unrecognized entity type code = "
483  << entity_type_code
484  << " in group "
485  << group_name
486  << std::endl;
487 
488  // Try to find this ID in the map from UNV element ids to libmesh ids.
489  auto it = _unv_elem_id_to_libmesh_elem_id.find(entity_tag);
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)
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  std::unique_ptr<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 = 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  auto it = _unv_node_id_to_libmesh_node_ptr.find(node_labels[j]);
834  if (it != _unv_node_id_to_libmesh_node_ptr.end())
835  elem->set_node(assign_elem_nodes[j]) = it->second;
836  else
837  libmesh_error_msg("ERROR: UNV node " << node_labels[j] << " not found!");
838  }
839 
840  elems_of_dimension[elem->dim()] = true;
841 
842  // Set the element's ID
843  elem->set_id(ctr);
844 
845  // Maintain a map from the libmesh (0-based) numbering to the
846  // UNV numbering.
847  _unv_elem_id_to_libmesh_elem_id[element_label] = ctr;
848 
849  // Add the element to the Mesh
850  mesh.add_elem(elem);
851 
852  // Increment the counter for the next iteration
853  ctr++;
854  } // end while (true)
855 }
856 
857 
858 
859 
860 
861 
862 void UNVIO::nodes_out (std::ostream & out_file)
863 {
864  if (this->verbose())
865  libMesh::out << " Writing " << MeshOutput<MeshBase>::mesh().n_nodes() << " nodes" << std::endl;
866 
867  // Write beginning of dataset
868  out_file << " -1\n"
869  << " "
871  << '\n';
872 
873  unsigned int
874  exp_coord_sys_dummy = 0, // export coordinate sys. (not supported yet)
875  disp_coord_sys_dummy = 0, // displacement coordinate sys. (not supp. yet)
876  color_dummy = 0; // color(not supported yet)
877 
878  // A reference to the parent class's mesh
880 
881  // Use scientific notation with capital E and 16 digits for printing out the coordinates
882  out_file << std::scientific << std::setprecision(16) << std::uppercase;
883 
884  for (const auto & current_node : mesh.node_ptr_range())
885  {
886  dof_id_type node_id = current_node->id();
887 
888  out_file << std::setw(10) << node_id
889  << std::setw(10) << exp_coord_sys_dummy
890  << std::setw(10) << disp_coord_sys_dummy
891  << std::setw(10) << color_dummy
892  << '\n';
893 
894  // The coordinates - always write out three coords regardless of LIBMESH_DIM
895  Real x = (*current_node)(0);
896 
897 #if LIBMESH_DIM > 1
898  Real y = (*current_node)(1);
899 #else
900  Real y = 0.;
901 #endif
902 
903 #if LIBMESH_DIM > 2
904  Real z = (*current_node)(2);
905 #else
906  Real z = 0.;
907 #endif
908 
909  out_file << std::setw(25) << x
910  << std::setw(25) << y
911  << std::setw(25) << z
912  << '\n';
913  }
914 
915  // Write end of dataset
916  out_file << " -1\n";
917 }
918 
919 
920 
921 
922 
923 
924 void UNVIO::elements_out(std::ostream & out_file)
925 {
926  if (this->verbose())
927  libMesh::out << " Writing elements" << std::endl;
928 
929  // Write beginning of dataset
930  out_file << " -1\n"
931  << " "
933  << "\n";
934 
935  unsigned
936  fe_descriptor_id = 0, // FE descriptor id
937  phys_prop_tab_dummy = 2, // physical property (not supported yet)
938  mat_prop_tab_dummy = 1, // material property (not supported yet)
939  color_dummy = 7; // color (not supported yet)
940 
941  // vector that assigns element nodes to their correct position
942  // currently only elements with up to 20 nodes
943  //
944  // Example:
945  // QUAD4 | 44:plane stress
946  // | linear quad
947  // position in libMesh | UNV numbering
948  // (note: 0-based) | (note: 1-based)
949  //
950  // assign_elem_node[0] = 0
951  // assign_elem_node[1] = 3
952  // assign_elem_node[2] = 2
953  // assign_elem_node[3] = 1
954  unsigned int assign_elem_nodes[20];
955 
956  unsigned int n_elem_written=0;
957 
958  // A reference to the parent class's mesh
960 
961  for (const auto & elem : mesh.element_ptr_range())
962  {
963  switch (elem->type())
964  {
965 
966  case TRI3:
967  {
968  fe_descriptor_id = 41; // Plane Stress Linear Triangle
969  assign_elem_nodes[0] = 0;
970  assign_elem_nodes[1] = 2;
971  assign_elem_nodes[2] = 1;
972  break;
973  }
974 
975  case TRI6:
976  {
977  fe_descriptor_id = 42; // Plane Stress Quadratic Triangle
978  assign_elem_nodes[0] = 0;
979  assign_elem_nodes[1] = 5;
980  assign_elem_nodes[2] = 2;
981  assign_elem_nodes[3] = 4;
982  assign_elem_nodes[4] = 1;
983  assign_elem_nodes[5] = 3;
984  break;
985  }
986 
987  case QUAD4:
988  {
989  fe_descriptor_id = 44; // Plane Stress Linear Quadrilateral
990  assign_elem_nodes[0] = 0;
991  assign_elem_nodes[1] = 3;
992  assign_elem_nodes[2] = 2;
993  assign_elem_nodes[3] = 1;
994  break;
995  }
996 
997  case QUAD8:
998  {
999  fe_descriptor_id = 45; // Plane Stress Quadratic Quadrilateral
1000  assign_elem_nodes[0] = 0;
1001  assign_elem_nodes[1] = 7;
1002  assign_elem_nodes[2] = 3;
1003  assign_elem_nodes[3] = 6;
1004  assign_elem_nodes[4] = 2;
1005  assign_elem_nodes[5] = 5;
1006  assign_elem_nodes[6] = 1;
1007  assign_elem_nodes[7] = 4;
1008  break;
1009  }
1010 
1011  case QUAD9:
1012  {
1013  fe_descriptor_id = 300; // Plane Stress Quadratic Quadrilateral
1014  assign_elem_nodes[0] = 0;
1015  assign_elem_nodes[1] = 7;
1016  assign_elem_nodes[2] = 3;
1017  assign_elem_nodes[3] = 6;
1018  assign_elem_nodes[4] = 2;
1019  assign_elem_nodes[5] = 5;
1020  assign_elem_nodes[6] = 1;
1021  assign_elem_nodes[7] = 4;
1022  assign_elem_nodes[8] = 8;
1023  break;
1024  }
1025 
1026  case TET4:
1027  {
1028  fe_descriptor_id = 111; // Solid Linear Tetrahedron
1029  assign_elem_nodes[0] = 0;
1030  assign_elem_nodes[1] = 1;
1031  assign_elem_nodes[2] = 2;
1032  assign_elem_nodes[3] = 3;
1033  break;
1034  }
1035 
1036  case PRISM6:
1037  {
1038  fe_descriptor_id = 112; // Solid Linear Prism
1039  assign_elem_nodes[0] = 0;
1040  assign_elem_nodes[1] = 1;
1041  assign_elem_nodes[2] = 2;
1042  assign_elem_nodes[3] = 3;
1043  assign_elem_nodes[4] = 4;
1044  assign_elem_nodes[5] = 5;
1045  break;
1046  }
1047 
1048  case HEX8:
1049  {
1050  fe_descriptor_id = 115; // Solid Linear Brick
1051  assign_elem_nodes[0] = 0;
1052  assign_elem_nodes[1] = 4;
1053  assign_elem_nodes[2] = 5;
1054  assign_elem_nodes[3] = 1;
1055  assign_elem_nodes[4] = 3;
1056  assign_elem_nodes[5] = 7;
1057  assign_elem_nodes[6] = 6;
1058  assign_elem_nodes[7] = 2;
1059  break;
1060  }
1061 
1062  case HEX20:
1063  {
1064  fe_descriptor_id = 116; // Solid Quadratic Brick
1065  assign_elem_nodes[ 0] = 0;
1066  assign_elem_nodes[ 1] = 12;
1067  assign_elem_nodes[ 2] = 4;
1068  assign_elem_nodes[ 3] = 16;
1069  assign_elem_nodes[ 4] = 5;
1070  assign_elem_nodes[ 5] = 13;
1071  assign_elem_nodes[ 6] = 1;
1072  assign_elem_nodes[ 7] = 8;
1073 
1074  assign_elem_nodes[ 8] = 11;
1075  assign_elem_nodes[ 9] = 19;
1076  assign_elem_nodes[10] = 17;
1077  assign_elem_nodes[11] = 9;
1078 
1079  assign_elem_nodes[12] = 3;
1080  assign_elem_nodes[13] = 15;
1081  assign_elem_nodes[14] = 7;
1082  assign_elem_nodes[15] = 18;
1083  assign_elem_nodes[16] = 6;
1084  assign_elem_nodes[17] = 14;
1085  assign_elem_nodes[18] = 2;
1086  assign_elem_nodes[19] = 10;
1087 
1088 
1089  break;
1090  }
1091 
1092  case TET10:
1093  {
1094  fe_descriptor_id = 118; // Solid Quadratic Tetrahedron
1095  assign_elem_nodes[0] = 0;
1096  assign_elem_nodes[1] = 4;
1097  assign_elem_nodes[2] = 1;
1098  assign_elem_nodes[3] = 5;
1099  assign_elem_nodes[4] = 2;
1100  assign_elem_nodes[5] = 6;
1101  assign_elem_nodes[6] = 7;
1102  assign_elem_nodes[7] = 8;
1103  assign_elem_nodes[8] = 9;
1104  assign_elem_nodes[9] = 3;
1105  break;
1106  }
1107 
1108  default:
1109  libmesh_error_msg("ERROR: Element type = " << elem->type() << " not supported in " << "UNVIO!");
1110  }
1111 
1112  dof_id_type elem_id = elem->id();
1113 
1114  out_file << std::setw(10) << elem_id // element ID
1115  << std::setw(10) << fe_descriptor_id // type of element
1116  << std::setw(10) << phys_prop_tab_dummy // not supported
1117  << std::setw(10) << mat_prop_tab_dummy // not supported
1118  << std::setw(10) << color_dummy // not supported
1119  << std::setw(10) << elem->n_nodes() // No. of nodes per element
1120  << '\n';
1121 
1122  for (auto j : elem->node_index_range())
1123  {
1124  // assign_elem_nodes[j]-th node: i.e., j loops over the
1125  // libMesh numbering, and assign_elem_nodes[j] over the
1126  // UNV numbering.
1127  const Node * node_in_unv_order = elem->node_ptr(assign_elem_nodes[j]);
1128 
1129  // new record after 8 id entries
1130  if (j==8 || j==16)
1131  out_file << '\n';
1132 
1133  // Write foreign label for this node
1134  dof_id_type node_id = node_in_unv_order->id();
1135 
1136  out_file << std::setw(10) << node_id;
1137  }
1138 
1139  out_file << '\n';
1140 
1141  n_elem_written++;
1142  }
1143 
1144  if (this->verbose())
1145  libMesh::out << " Finished writing " << n_elem_written << " elements" << std::endl;
1146 
1147  // Write end of dataset
1148  out_file << " -1\n";
1149 }
1150 
1151 
1152 
1153 void UNVIO::read_dataset(std::string file_name)
1154 {
1155  std::ifstream in_stream(file_name.c_str());
1156 
1157  if (!in_stream.good())
1158  libmesh_error_msg("Error opening UNV data file.");
1159 
1160  std::string olds, news, dummy;
1161 
1162  while (true)
1163  {
1164  in_stream >> olds >> news;
1165 
1166  // A "-1" followed by a number means the beginning of a dataset.
1167  while (((olds != "-1") || (news == "-1")) && !in_stream.eof())
1168  {
1169  olds = news;
1170  in_stream >> news;
1171  }
1172 
1173  if (in_stream.eof())
1174  break;
1175 
1176  // We only support reading datasets in the "2414" format.
1177  if (news == "2414")
1178  {
1179  // Ignore the rest of the current line and the next two
1180  // lines that contain analysis dataset label and name.
1181  for (unsigned int i=0; i<3; i++)
1182  std::getline(in_stream, dummy);
1183 
1184  // Read the dataset location, where
1185  // 1: Data at nodes
1186  // 2: Data on elements
1187  // Currently only data on nodes is supported.
1188  unsigned int dataset_location;
1189  in_stream >> dataset_location;
1190 
1191  // Currently only nodal datasets are supported.
1192  if (dataset_location != 1)
1193  libmesh_error_msg("ERROR: Currently only Data at nodes is supported.");
1194 
1195  // Ignore the rest of this line and the next five records.
1196  for (unsigned int i=0; i<6; i++)
1197  std::getline(in_stream, dummy);
1198 
1199  // These data are all of no interest to us...
1200  unsigned int
1201  model_type,
1202  analysis_type,
1203  data_characteristic,
1204  result_type;
1205 
1206  // The type of data (complex, real, float, double etc, see
1207  // below).
1208  unsigned int data_type;
1209 
1210  // The number of floating-point values per entity.
1211  unsigned int num_vals_per_node;
1212 
1213  in_stream >> model_type // not used here
1214  >> analysis_type // not used here
1215  >> data_characteristic // not used here
1216  >> result_type // not used here
1217  >> data_type
1218  >> num_vals_per_node;
1219 
1220  // Ignore the rest of record 9, and records 10-13.
1221  for (unsigned int i=0; i<5; i++)
1222  std::getline(in_stream, dummy);
1223 
1224  // Now get the foreign (aka UNV node) node number and
1225  // the respective nodal data.
1226  int f_n_id;
1227  std::vector<Number> values;
1228 
1229  while (true)
1230  {
1231  in_stream >> f_n_id;
1232 
1233  // If -1 then we have reached the end of the dataset.
1234  if (f_n_id == -1)
1235  break;
1236 
1237  // Resize the values vector (usually data in three
1238  // principle directions, i.e. num_vals_per_node = 3).
1239  values.resize(num_vals_per_node);
1240 
1241  // Read the values for the respective node.
1242  for (unsigned int data_cnt=0; data_cnt<num_vals_per_node; data_cnt++)
1243  {
1244  // Check what data type we are reading.
1245  // 2,4: Real
1246  // 5,6: Complex
1247  // other data types are not supported yet.
1248  // As again, these floats may also be written
1249  // using a 'D' instead of an 'e'.
1250  if (data_type == 2 || data_type == 4)
1251  {
1252  std::string buf;
1253  in_stream >> buf;
1254  this->need_D_to_e(buf);
1255 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1256  values[data_cnt] = Complex(std::atof(buf.c_str()), 0.);
1257 #else
1258  values[data_cnt] = std::atof(buf.c_str());
1259 #endif
1260  }
1261 
1262  else if (data_type == 5 || data_type == 6)
1263  {
1264 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1265  Real re_val, im_val;
1266 
1267  std::string buf;
1268  in_stream >> buf;
1269 
1270  if (this->need_D_to_e(buf))
1271  {
1272  re_val = std::atof(buf.c_str());
1273  in_stream >> buf;
1274  this->need_D_to_e(buf);
1275  im_val = std::atof(buf.c_str());
1276  }
1277  else
1278  {
1279  re_val = std::atof(buf.c_str());
1280  in_stream >> im_val;
1281  }
1282 
1283  values[data_cnt] = Complex(re_val,im_val);
1284 #else
1285 
1286  libmesh_error_msg("ERROR: Complex data only supported when libMesh is configured with --enable-complex!");
1287 #endif
1288  }
1289 
1290  else
1291  libmesh_error_msg("ERROR: Data type not supported.");
1292 
1293  } // end loop data_cnt
1294 
1295  // Get a pointer to the Node associated with the UNV node id.
1296  auto it = _unv_node_id_to_libmesh_node_ptr.find(f_n_id);
1297 
1298  if (it == _unv_node_id_to_libmesh_node_ptr.end())
1299  libmesh_error_msg("UNV node id " << f_n_id << " was not found.");
1300 
1301  // Store the nodal values in our map which uses the
1302  // libMesh Node* as the key. We use operator[] here
1303  // because we want to create an empty vector if the
1304  // entry does not already exist.
1305  _node_data[it->second] = values;
1306  } // end while (true)
1307  } // end if (news == "2414")
1308  } // end while (true)
1309 }
1310 
1311 
1312 
1313 const std::vector<Number> *
1314 UNVIO::get_data (Node * node) const
1315 {
1316  auto it = _node_data.find(node);
1317 
1318  if (it == _node_data.end())
1319  return nullptr;
1320  else
1321  return &(it->second);
1322 }
1323 
1324 
1325 } // namespace libMesh
const MT & mesh() const
Definition: mesh_output.h:234
void write_implementation(std::ostream &out_stream)
Definition: unv_io.C:289
bool _verbose
Definition: unv_io.h:182
A 2D triangular element with 3 nodes.
Definition: face_tri3.h:56
UNVIO(MeshBase &mesh)
Definition: unv_io.C:66
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
bool & verbose()
Definition: unv_io.C:89
A 2D quadrilateral element with 4 nodes.
Definition: face_quad4.h:51
A 3D hexahedral element with 8 nodes.
Definition: cell_hex8.h:53
std::map< dof_id_type, Node * > _unv_node_id_to_libmesh_node_ptr
Definition: unv_io.h:188
static const std::string _elements_dataset_label
Definition: unv_io.h:198
std::vector< bool > elems_of_dimension
Definition: mesh_input.h:97
A 3D prismatic element with 6 nodes.
Definition: cell_prism6.h:52
void nodes_in(std::istream &in_file)
Definition: unv_io.C:302
void read_implementation(std::istream &in_stream)
Definition: unv_io.C:122
A 3D hexahedral element with 20 nodes.
Definition: cell_hex20.h:68
void read_dataset(std::string file_name)
Definition: unv_io.C:1153
unsigned short int side
Definition: xdr_io.C:50
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
std::map< Node *, std::vector< Number > > _node_data
Definition: unv_io.h:214
A 2D quadrilateral element with 8 nodes.
Definition: face_quad8.h:51
MPI_Datatype data_type
Definition: data_type.h:46
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
unsigned char max_elem_dimension_seen()
Definition: unv_io.C:375
void nodes_out(std::ostream &out_file)
Definition: unv_io.C:862
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:131
static const std::string _groups_dataset_label
Definition: unv_io.h:203
void elements_out(std::ostream &out_file)
Definition: unv_io.C:924
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
dof_id_type & set_id()
Definition: dof_object.h:664
A 2D triangular element with 6 nodes.
Definition: face_tri6.h:56
virtual void write(const std::string &) override
Definition: unv_io.C:260
const dof_id_type n_nodes
Definition: tecplot_io.C:68
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
virtual Elem * add_elem(Elem *e)=0
virtual dof_id_type key(const unsigned int s) const =0
virtual SimpleRange< node_iterator > node_ptr_range()=0
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:538
OStreamProxy err(std::cerr)
A 2D quadrilateral element with 9 nodes.
Definition: face_quad9.h:51
std::map< unsigned, unsigned > _unv_elem_id_to_libmesh_elem_id
Definition: unv_io.h:208
std::string & sideset_name(boundary_id_type id)
virtual const Elem * elem_ptr(const dof_id_type i) const =0
std::complex< Real > Complex
bool need_D_to_e(std::string &number)
Definition: unv_io.C:391
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
A 3D tetrahedral element with 10 nodes.
Definition: cell_tet10.h:60
virtual unsigned short dim() const =0
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
A 1D geometric element with 2 nodes.
Definition: edge_edge2.h:43
const std::vector< Number > * get_data(Node *node) const
Definition: unv_io.C:1314
virtual ~UNVIO()
Definition: unv_io.C:83
A 3D tetrahedral element with 4 nodes.
Definition: cell_tet4.h:53
static const std::string _nodes_dataset_label
Definition: unv_io.h:193
OStreamProxy out(std::cout)
virtual void read(const std::string &) override
Definition: unv_io.C:96
A geometric point in (x,y,z) space.
Definition: point.h:38
void elements_in(std::istream &in_file)
Definition: unv_io.C:574
uint8_t dof_id_type
Definition: id_types.h:64
void groups_in(std::istream &in_file)
Definition: unv_io.C:408