Reference documentation for deal.II version Git ac8d010384 2020-11-27 19:49:05 +0100
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_in.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
19 #include <deal.II/base/utilities.h>
20 
21 #include <deal.II/grid/grid_in.h>
24 #include <deal.II/grid/tria.h>
25 
26 #include <boost/archive/binary_iarchive.hpp>
27 #include <boost/io/ios_state.hpp>
28 #include <boost/property_tree/ptree.hpp>
29 #include <boost/property_tree/xml_parser.hpp>
30 #include <boost/serialization/serialization.hpp>
31 
32 #include <algorithm>
33 #include <cctype>
34 #include <fstream>
35 #include <functional>
36 #include <map>
37 
38 #ifdef DEAL_II_WITH_ASSIMP
39 # include <assimp/Importer.hpp> // C++ importer interface
40 # include <assimp/postprocess.h> // Post processing flags
41 # include <assimp/scene.h> // Output data structure
42 #endif
43 
44 #ifdef DEAL_II_TRILINOS_WITH_SEACAS
45 # include <exodusII.h>
46 #endif
47 
48 
50 
51 
52 namespace
53 {
62  template <int spacedim>
63  void
64  assign_1d_boundary_ids(
65  const std::map<unsigned int, types::boundary_id> &boundary_ids,
67  {
68  if (boundary_ids.size() > 0)
69  for (const auto &cell : triangulation.active_cell_iterators())
70  for (unsigned int f : GeometryInfo<1>::face_indices())
71  if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
72  {
74  cell->at_boundary(f),
75  ExcMessage(
76  "You are trying to prescribe boundary ids on the face "
77  "of a 1d cell (i.e., on a vertex), but this face is not actually at "
78  "the boundary of the mesh. This is not allowed."));
79  cell->face(f)->set_boundary_id(
80  boundary_ids.find(cell->vertex_index(f))->second);
81  }
82  }
83 
84 
85  template <int dim, int spacedim>
86  void
87  assign_1d_boundary_ids(const std::map<unsigned int, types::boundary_id> &,
89  {
90  // we shouldn't get here since boundary ids are not assigned to
91  // vertices except in 1d
92  Assert(dim != 1, ExcInternalError());
93  }
94 } // namespace
95 
96 template <int dim, int spacedim>
98  : tria(nullptr, typeid(*this).name())
100 {}
101 
102 
103 
104 template <int dim, int spacedim>
106  : tria(&t, typeid(*this).name())
108 {}
109 
110 
111 
112 template <int dim, int spacedim>
113 void
115 {
116  tria = &t;
117 }
118 
119 
120 
121 template <int dim, int spacedim>
122 void
124 {
125  std::string line;
126 
127  // verify that the first, third and fourth lines match
128  // expectations. the second line of the file may essentially be
129  // anything the author of the file chose to identify what's in
130  // there, so we just ensure that we can read it
131  {
132  std::string text[4];
133  text[0] = "# vtk DataFile Version 3.0";
134  text[1] = "****";
135  text[2] = "ASCII";
136  text[3] = "DATASET UNSTRUCTURED_GRID";
137 
138  for (unsigned int i = 0; i < 4; ++i)
139  {
140  getline(in, line);
141  if (i != 1)
142  AssertThrow(
143  line.compare(text[i]) == 0,
144  ExcMessage(
145  std::string(
146  "While reading VTK file, failed to find a header line with text <") +
147  text[i] + ">"));
148  }
149  }
150 
152 
153  std::vector<Point<spacedim>> vertices;
154  std::vector<CellData<dim>> cells;
155  SubCellData subcelldata;
156 
157  std::string keyword;
158 
159  in >> keyword;
160 
162 
163  if (keyword == "POINTS")
164  {
165  unsigned int n_vertices;
166  in >> n_vertices;
167 
168  in >> keyword; // float, double, int, char, etc.
169 
170  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
171  {
172  // VTK format always specifies vertex coordinates with 3 components
173  Point<3> x;
174  in >> x(0) >> x(1) >> x(2);
175 
176  vertices.emplace_back();
177  for (unsigned int d = 0; d < spacedim; ++d)
178  vertices.back()(d) = x(d);
179  }
180  }
181 
182  else
183  AssertThrow(false,
184  ExcMessage(
185  "While reading VTK file, failed to find POINTS section"));
186 
187  in >> keyword;
188 
189  unsigned int n_geometric_objects = 0;
190  unsigned int n_ints;
191 
192  bool is_quad_or_hex_mesh = false;
193  bool is_tria_or_tet_mesh = false;
194 
195  if (keyword == "CELLS")
196  {
197  // jump to the `CELL_TYPES` section and read in cell types
198  std::vector<unsigned int> cell_types;
199  {
200  std::streampos oldpos = in.tellg();
201 
202 
203  while (in >> keyword)
204  if (keyword == "CELL_TYPES")
205  {
206  in >> n_ints;
207 
208  cell_types.resize(n_ints);
209 
210  for (unsigned int i = 0; i < n_ints; ++i)
211  in >> cell_types[i];
212 
213  break;
214  }
215 
216  in.seekg(oldpos);
217  }
218 
219  in >> n_geometric_objects;
220  in >> n_ints; // Ignore this, since we don't need it.
221 
222  if (dim == 3)
223  {
224  for (unsigned int count = 0; count < n_geometric_objects; count++)
225  {
226  unsigned int type;
227  in >> type;
228 
229  if (cell_types[count] == 10 || cell_types[count] == 12)
230  {
231  if (cell_types[count] == 10)
232  is_tria_or_tet_mesh = true;
233  if (cell_types[count] == 12)
234  is_quad_or_hex_mesh = true;
235 
236  // we assume that the file contains first all cells,
237  // and only then any faces or lines
238  AssertThrow(subcelldata.boundary_quads.size() == 0 &&
239  subcelldata.boundary_lines.size() == 0,
241 
242  cells.emplace_back(type);
243 
244  for (unsigned int j = 0; j < type; j++) // loop to feed data
245  in >> cells.back().vertices[j];
246 
247  cells.back().material_id = 0;
248  }
249 
250  else if (cell_types[count] == 5 || cell_types[count] == 9)
251  {
252  if (cell_types[count] == 5)
253  is_tria_or_tet_mesh = true;
254  if (cell_types[count] == 9)
255  is_quad_or_hex_mesh = true;
256 
257  // we assume that the file contains first all cells,
258  // then all faces, and finally all lines
259  AssertThrow(subcelldata.boundary_lines.size() == 0,
261 
262  subcelldata.boundary_quads.emplace_back(type);
263 
264  for (unsigned int j = 0; j < type;
265  j++) // loop to feed the data to the boundary
266  in >> subcelldata.boundary_quads.back().vertices[j];
267 
268  subcelldata.boundary_quads.back().material_id = 0;
269  }
270  else if (cell_types[count] == 3)
271  {
272  subcelldata.boundary_lines.emplace_back(type);
273 
274  for (unsigned int j = 0; j < type;
275  j++) // loop to feed the data to the boundary
276  in >> subcelldata.boundary_lines.back().vertices[j];
277 
278  subcelldata.boundary_lines.back().material_id = 0;
279  }
280 
281  else
282  AssertThrow(
283  false,
284  ExcMessage(
285  "While reading VTK file, unknown file type encountered"));
286  }
287  }
288  else if (dim == 2)
289  {
290  for (unsigned int count = 0; count < n_geometric_objects; count++)
291  {
292  unsigned int type;
293  in >> type;
294 
295  if (cell_types[count] == 5 || cell_types[count] == 9)
296  {
297  // we assume that the file contains first all cells,
298  // and only then any faces
299  AssertThrow(subcelldata.boundary_lines.size() == 0,
301 
302  if (cell_types[count] == 5)
303  is_tria_or_tet_mesh = true;
304  if (cell_types[count] == 9)
305  is_quad_or_hex_mesh = true;
306 
307  cells.emplace_back(type);
308 
309  for (unsigned int j = 0; j < type; j++) // loop to feed data
310  in >> cells.back().vertices[j];
311 
312  cells.back().material_id = 0;
313  }
314 
315  else if (cell_types[count] == 3)
316  {
317  // If this is encountered, the pointer comes out of the loop
318  // and starts processing boundaries.
319  subcelldata.boundary_lines.emplace_back(type);
320 
321  for (unsigned int j = 0; j < type;
322  j++) // loop to feed the data to the boundary
323  {
324  in >> subcelldata.boundary_lines.back().vertices[j];
325  }
326 
327  subcelldata.boundary_lines.back().material_id = 0;
328  }
329 
330  else
331  AssertThrow(
332  false,
333  ExcMessage(
334  "While reading VTK file, unknown cell type encountered"));
335  }
336  }
337  else if (dim == 1)
338  {
339  for (unsigned int count = 0; count < n_geometric_objects; count++)
340  {
341  unsigned int type;
342  in >> type;
343 
344  AssertThrow(
345  cell_types[count] == 3 && type == 2,
346  ExcMessage(
347  "While reading VTK file, unknown cell type encountered"));
348  cells.emplace_back(type);
349 
350  for (unsigned int j = 0; j < type; j++) // loop to feed data
351  in >> cells.back().vertices[j];
352 
353  cells.back().material_id = 0;
354  }
355  }
356  else
357  AssertThrow(false,
358  ExcMessage(
359  "While reading VTK file, failed to find CELLS section"));
360 
363 
364  in >> keyword;
365 
366  AssertThrow(
367  keyword == "CELL_TYPES",
368  ExcMessage(std::string(
369  "While reading VTK file, missing CELL_TYPES section. Found <" +
370  keyword + "> instead.")));
371 
372  in >> n_ints;
373  AssertThrow(
374  n_ints == n_geometric_objects,
375  ExcMessage("The VTK reader found a CELL_DATA statement "
376  "that lists a total of " +
377  Utilities::int_to_string(n_ints) +
378  " cell data objects, but this needs to "
379  "equal the number of cells (which is " +
380  Utilities::int_to_string(cells.size()) +
381  ") plus the number of quads (" +
382  Utilities::int_to_string(subcelldata.boundary_quads.size()) +
383  " in 3d or the number of lines (" +
384  Utilities::int_to_string(subcelldata.boundary_lines.size()) +
385  ") in 2d."));
386 
387  int tmp_int;
388  for (unsigned int i = 0; i < n_ints; ++i)
389  in >> tmp_int;
390 
391  // Ignore everything up to CELL_DATA
392  while (in >> keyword)
393  if (keyword == "CELL_DATA")
394  {
395  unsigned int n_ids;
396  in >> n_ids;
397 
398  AssertThrow(n_ids == n_geometric_objects,
399  ExcMessage("The VTK reader found a CELL_DATA statement "
400  "that lists a total of " +
401  Utilities::int_to_string(n_ids) +
402  " cell data objects, but this needs to "
403  "equal the number of cells (which is " +
404  Utilities::int_to_string(cells.size()) +
405  ") plus the number of quads (" +
407  subcelldata.boundary_quads.size()) +
408  " in 3d or the number of lines (" +
410  subcelldata.boundary_lines.size()) +
411  ") in 2d."));
412 
413  const std::vector<std::string> data_sets{"MaterialID",
414  "ManifoldID"};
415 
416  for (unsigned int i = 0; i < data_sets.size(); ++i)
417  {
418  // Ignore everything until we get to a SCALARS data set
419  while (in >> keyword)
420  if (keyword == "SCALARS")
421  {
422  // Now see if we know about this type of data set,
423  // if not, just ignore everything till the next SCALARS
424  // keyword
425  std::string field_name;
426  in >> field_name;
427  if (std::find(data_sets.begin(),
428  data_sets.end(),
429  field_name) == data_sets.end())
430  // The data set here is not one of the ones we know, so
431  // keep ignoring everything until the next SCALARS
432  // keyword.
433  continue;
434 
435  // Now we got somewhere. Proceed from here, assert
436  // that the type of the table is int, and ignore the
437  // rest of the line.
438  // SCALARS MaterialID int 1
439  // (the last number is optional)
440  std::string line;
441  std::getline(in, line);
442  AssertThrow(
443  line.substr(1,
444  std::min(static_cast<std::size_t>(3),
445  line.size() - 1)) == "int",
446  ExcMessage(
447  "While reading VTK file, material- and manifold IDs can only have type 'int'."));
448 
449  in >> keyword;
450  AssertThrow(
451  keyword == "LOOKUP_TABLE",
452  ExcMessage(
453  "While reading VTK file, missing keyword 'LOOKUP_TABLE'."));
454 
455  in >> keyword;
456  AssertThrow(
457  keyword == "default",
458  ExcMessage(
459  "While reading VTK file, missing keyword 'default'."));
460 
461  // read material or manifold ids first for all cells,
462  // then for all faces, and finally for all lines. the
463  // assumption that cells come before all faces and
464  // lines has been verified above via an assertion, so
465  // the order used in the following blocks makes sense
466  for (unsigned int i = 0; i < cells.size(); i++)
467  {
468  int id;
469  in >> id;
470  if (field_name == "MaterialID")
471  cells[i].material_id =
472  static_cast<types::material_id>(id);
473  else if (field_name == "ManifoldID")
474  cells[i].manifold_id =
475  static_cast<types::manifold_id>(id);
476  else
477  Assert(false, ExcInternalError());
478  }
479 
480  if (dim == 3)
481  {
482  for (auto &boundary_quad : subcelldata.boundary_quads)
483  {
484  int id;
485  in >> id;
486  if (field_name == "MaterialID")
487  boundary_quad.material_id =
488  static_cast<types::material_id>(id);
489  else if (field_name == "ManifoldID")
490  boundary_quad.manifold_id =
491  static_cast<types::manifold_id>(id);
492  else
493  Assert(false, ExcInternalError());
494  }
495  for (auto &boundary_line : subcelldata.boundary_lines)
496  {
497  int id;
498  in >> id;
499  if (field_name == "MaterialID")
500  boundary_line.material_id =
501  static_cast<types::material_id>(id);
502  else if (field_name == "ManifoldID")
503  boundary_line.manifold_id =
504  static_cast<types::manifold_id>(id);
505  else
506  Assert(false, ExcInternalError());
507  }
508  }
509  else if (dim == 2)
510  {
511  for (auto &boundary_line : subcelldata.boundary_lines)
512  {
513  int id;
514  in >> id;
515  if (field_name == "MaterialID")
516  boundary_line.material_id =
517  static_cast<types::material_id>(id);
518  else if (field_name == "ManifoldID")
519  boundary_line.manifold_id =
520  static_cast<types::manifold_id>(id);
521  else
522  Assert(false, ExcInternalError());
523  }
524  }
525  }
526  }
527  }
528 
529  Assert(subcelldata.check_consistency(dim), ExcInternalError());
530 
531 
532  // make sure that only either simplex or hypercube cells are available
533  //
534  // TODO: the functions below (GridTools::delete_unused_vertices(),
535  // GridTools::invert_all_cells_of_negative_grid(),
536  // GridReordering::reorder_cells(),
537  // Triangulation::create_triangulation_compatibility()) need to be
538  // revisited for simplex meshes
539  AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
541 
542  if (dim == 1 || is_quad_or_hex_mesh)
543  {
544  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
545 
546  if (dim == spacedim)
548  vertices, cells);
549 
551  tria->create_triangulation_compatibility(vertices,
552  cells,
553  subcelldata);
554 
555  return;
556  }
557  else
558  {
559  tria->create_triangulation(vertices, cells, subcelldata);
560  }
561  }
562  else
563  AssertThrow(false,
564  ExcMessage(
565  "While reading VTK file, failed to find CELLS section"));
566 }
567 
568 
569 
570 template <int dim, int spacedim>
571 void
573 {
574  namespace pt = boost::property_tree;
575  pt::ptree tree;
576  pt::read_xml(in, tree);
577  auto section = tree.get_optional<std::string>("VTKFile.dealiiData");
578 
579  AssertThrow(section,
580  ExcMessage(
581  "While reading a VTU file, failed to find dealiiData section. "
582  "Notice that we can only read grid files in .vtu format that "
583  "were created by the deal.II library, using a call to "
584  "GridOut::write_vtu(), where the flag "
585  "GridOutFlags::Vtu::serialize_triangulation is set to true."));
586 
587  const auto decoded =
588  Utilities::decode_base64({section->begin(), section->end() - 1});
589  const auto string_archive =
590  Utilities::decompress({decoded.begin(), decoded.end()});
591  std::istringstream in_stream(string_archive);
592  boost::archive::binary_iarchive ia(in_stream);
593  tria->load(ia, 0);
594 }
595 
596 
597 template <int dim, int spacedim>
598 void
600 {
601  Assert(tria != nullptr, ExcNoTriangulationSelected());
602  Assert((dim == 2) || (dim == 3), ExcNotImplemented());
603 
604  AssertThrow(in, ExcIO());
605  skip_comment_lines(in, '#'); // skip comments (if any) at beginning of file
606 
607  int tmp;
608 
609  AssertThrow(in, ExcIO());
610  in >> tmp;
611  AssertThrow(in, ExcIO());
612  in >> tmp;
613 
614  // section 2411 describes vertices: see
615  // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2411
616  AssertThrow(tmp == 2411, ExcUnknownSectionType(tmp));
617 
618  std::vector<Point<spacedim>> vertices; // vector of vertex coordinates
619  std::map<int, int>
620  vertex_indices; // # vert in unv (key) ---> # vert in deal.II (value)
621 
622  int no_vertex = 0; // deal.II
623 
624  while (tmp != -1) // we do until reach end of 2411
625  {
626  int no; // unv
627  int dummy;
628  double x[3];
629 
630  AssertThrow(in, ExcIO());
631  in >> no;
632 
633  tmp = no;
634  if (tmp == -1)
635  break;
636 
637  in >> dummy >> dummy >> dummy;
638 
639  AssertThrow(in, ExcIO());
640  in >> x[0] >> x[1] >> x[2];
641 
642  vertices.emplace_back();
643 
644  for (unsigned int d = 0; d < spacedim; d++)
645  vertices.back()(d) = x[d];
646 
647  vertex_indices[no] = no_vertex;
648 
649  no_vertex++;
650  }
651 
652  AssertThrow(in, ExcIO());
653  in >> tmp;
654  AssertThrow(in, ExcIO());
655  in >> tmp;
656 
657  // section 2412 describes elements: see
658  // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2412
659  AssertThrow(tmp == 2412, ExcUnknownSectionType(tmp));
660 
661  std::vector<CellData<dim>> cells; // vector of cells
662  SubCellData subcelldata;
663 
664  std::map<int, int>
665  cell_indices; // # cell in unv (key) ---> # cell in deal.II (value)
666  std::map<int, int>
667  line_indices; // # line in unv (key) ---> # line in deal.II (value)
668  std::map<int, int>
669  quad_indices; // # quad in unv (key) ---> # quad in deal.II (value)
670 
671  int no_cell = 0; // deal.II
672  int no_line = 0; // deal.II
673  int no_quad = 0; // deal.II
674 
675  while (tmp != -1) // we do until reach end of 2412
676  {
677  int no; // unv
678  int type;
679  int dummy;
680 
681  AssertThrow(in, ExcIO());
682  in >> no;
683 
684  tmp = no;
685  if (tmp == -1)
686  break;
687 
688  in >> type >> dummy >> dummy >> dummy >> dummy;
689 
690  AssertThrow((type == 11) || (type == 44) || (type == 94) || (type == 115),
691  ExcUnknownElementType(type));
692 
693  if ((((type == 44) || (type == 94)) && (dim == 2)) ||
694  ((type == 115) && (dim == 3))) // cell
695  {
696  cells.emplace_back();
697 
698  AssertThrow(in, ExcIO());
699  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
700  in >> cells.back().vertices[v];
701 
702  cells.back().material_id = 0;
703 
704  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
705  cells.back().vertices[v] = vertex_indices[cells.back().vertices[v]];
706 
707  cell_indices[no] = no_cell;
708 
709  no_cell++;
710  }
711  else if (((type == 11) && (dim == 2)) ||
712  ((type == 11) && (dim == 3))) // boundary line
713  {
714  AssertThrow(in, ExcIO());
715  in >> dummy >> dummy >> dummy;
716 
717  subcelldata.boundary_lines.emplace_back();
718 
719  AssertThrow(in, ExcIO());
720  for (unsigned int &vertex :
721  subcelldata.boundary_lines.back().vertices)
722  in >> vertex;
723 
724  subcelldata.boundary_lines.back().material_id = 0;
725 
726  for (unsigned int &vertex :
727  subcelldata.boundary_lines.back().vertices)
728  vertex = vertex_indices[vertex];
729 
730  line_indices[no] = no_line;
731 
732  no_line++;
733  }
734  else if (((type == 44) || (type == 94)) && (dim == 3)) // boundary quad
735  {
736  subcelldata.boundary_quads.emplace_back();
737 
738  AssertThrow(in, ExcIO());
739  for (unsigned int &vertex :
740  subcelldata.boundary_quads.back().vertices)
741  in >> vertex;
742 
743  subcelldata.boundary_quads.back().material_id = 0;
744 
745  for (unsigned int &vertex :
746  subcelldata.boundary_quads.back().vertices)
747  vertex = vertex_indices[vertex];
748 
749  quad_indices[no] = no_quad;
750 
751  no_quad++;
752  }
753  else
754  AssertThrow(false,
755  ExcMessage("Unknown element label <" +
757  "> when running in dim=" +
759  }
760 
761  // note that so far all materials and bcs are explicitly set to 0
762  // if we do not need more info on materials and bcs - this is end of file
763  // if we do - section 2467 or 2477 comes
764 
765  in >> tmp; // tmp can be either -1 or end-of-file
766 
767  if (!in.eof())
768  {
769  AssertThrow(in, ExcIO());
770  in >> tmp;
771 
772  // section 2467 (2477) describes (materials - first and bcs - second) or
773  // (bcs - first and materials - second) - sequence depends on which
774  // group is created first: see
775  // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2467
776  AssertThrow((tmp == 2467) || (tmp == 2477), ExcUnknownSectionType(tmp));
777 
778  while (tmp != -1) // we do until reach end of 2467 or 2477
779  {
780  int n_entities; // number of entities in group
781  int id; // id is either material or bc
782  int no; // unv
783  int dummy;
784 
785  AssertThrow(in, ExcIO());
786  in >> dummy;
787 
788  tmp = dummy;
789  if (tmp == -1)
790  break;
791 
792  in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
793  n_entities;
794 
795  AssertThrow(in, ExcIO());
796  in >> id;
797 
798  const unsigned int n_lines =
799  (n_entities % 2 == 0) ? (n_entities / 2) : ((n_entities + 1) / 2);
800 
801  for (unsigned int line = 0; line < n_lines; line++)
802  {
803  unsigned int n_fragments;
804 
805  if (line == n_lines - 1)
806  n_fragments = (n_entities % 2 == 0) ? (2) : (1);
807  else
808  n_fragments = 2;
809 
810  for (unsigned int no_fragment = 0; no_fragment < n_fragments;
811  no_fragment++)
812  {
813  AssertThrow(in, ExcIO());
814  in >> dummy >> no >> dummy >> dummy;
815 
816  if (cell_indices.count(no) > 0) // cell - material
817  cells[cell_indices[no]].material_id = id;
818 
819  if (line_indices.count(no) > 0) // boundary line - bc
820  subcelldata.boundary_lines[line_indices[no]].material_id =
821  id;
822 
823  if (quad_indices.count(no) > 0) // boundary quad - bc
824  subcelldata.boundary_quads[quad_indices[no]].material_id =
825  id;
826  }
827  }
828  }
829  }
830 
831  Assert(subcelldata.check_consistency(dim), ExcInternalError());
832 
833  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
834 
835  if (dim == spacedim)
837  cells);
838 
840 
841  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
842 }
843 
844 
845 
846 template <int dim, int spacedim>
847 void
849  const bool apply_all_indicators_to_manifolds)
850 {
851  Assert(tria != nullptr, ExcNoTriangulationSelected());
852  AssertThrow(in, ExcIO());
853 
854  // skip comments at start of file
855  skip_comment_lines(in, '#');
856 
857 
858  unsigned int n_vertices;
859  unsigned int n_cells;
860  int dummy;
861 
862  in >> n_vertices >> n_cells >> dummy // number of data vectors
863  >> dummy // cell data
864  >> dummy; // model data
865  AssertThrow(in, ExcIO());
866 
867  // set up array of vertices
868  std::vector<Point<spacedim>> vertices(n_vertices);
869  // set up mapping between numbering
870  // in ucd-file (key) and in the
871  // vertices vector
872  std::map<int, int> vertex_indices;
873 
874  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
875  {
876  int vertex_number;
877  double x[3];
878 
879  // read vertex
880  AssertThrow(in, ExcIO());
881  in >> vertex_number >> x[0] >> x[1] >> x[2];
882 
883  // store vertex
884  for (unsigned int d = 0; d < spacedim; ++d)
885  vertices[vertex](d) = x[d];
886  // store mapping; note that
887  // vertices_indices[i] is automatically
888  // created upon first usage
889  vertex_indices[vertex_number] = vertex;
890  }
891 
892  // set up array of cells
893  std::vector<CellData<dim>> cells;
894  SubCellData subcelldata;
895 
896  for (unsigned int cell = 0; cell < n_cells; ++cell)
897  {
898  // note that since in the input
899  // file we found the number of
900  // cells at the top, there
901  // should still be input here,
902  // so check this:
903  AssertThrow(in, ExcIO());
904 
905  std::string cell_type;
906 
907  // we use an unsigned int because we
908  // fill this variable through an read-in process
909  unsigned int material_id;
910 
911  in >> dummy // cell number
912  >> material_id;
913  in >> cell_type;
914 
915  if (((cell_type == "line") && (dim == 1)) ||
916  ((cell_type == "quad") && (dim == 2)) ||
917  ((cell_type == "hex") && (dim == 3)))
918  // found a cell
919  {
920  // allocate and read indices
921  cells.emplace_back();
922  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
923  in >> cells.back().vertices[i];
924 
925  // to make sure that the cast won't fail
927  ExcIndexRange(material_id,
928  0,
930  // we use only material_ids in the range from 0 to
931  // numbers::invalid_material_id-1
933 
934  if (apply_all_indicators_to_manifolds)
935  cells.back().manifold_id =
936  static_cast<types::manifold_id>(material_id);
937  cells.back().material_id = material_id;
938 
939  // transform from ucd to
940  // consecutive numbering
941  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
942  if (vertex_indices.find(cells.back().vertices[i]) !=
943  vertex_indices.end())
944  // vertex with this index exists
945  cells.back().vertices[i] =
946  vertex_indices[cells.back().vertices[i]];
947  else
948  {
949  // no such vertex index
950  AssertThrow(false,
952  cells.back().vertices[i]));
953 
954  cells.back().vertices[i] = numbers::invalid_unsigned_int;
955  }
956  }
957  else if ((cell_type == "line") && ((dim == 2) || (dim == 3)))
958  // boundary info
959  {
960  subcelldata.boundary_lines.emplace_back();
961  in >> subcelldata.boundary_lines.back().vertices[0] >>
962  subcelldata.boundary_lines.back().vertices[1];
963 
964  // to make sure that the cast won't fail
966  ExcIndexRange(material_id,
967  0,
969  // we use only boundary_ids in the range from 0 to
970  // numbers::internal_face_boundary_id-1
972 
973  // Make sure to set both manifold id and boundary id appropriately in
974  // both cases:
975  // numbers::internal_face_boundary_id and numbers::flat_manifold_id
976  // are ignored in Triangulation::create_triangulation.
977  if (apply_all_indicators_to_manifolds)
978  {
979  subcelldata.boundary_lines.back().boundary_id =
981  subcelldata.boundary_lines.back().manifold_id =
982  static_cast<types::manifold_id>(material_id);
983  }
984  else
985  {
986  subcelldata.boundary_lines.back().boundary_id =
987  static_cast<types::boundary_id>(material_id);
988  subcelldata.boundary_lines.back().manifold_id =
990  }
991 
992  // transform from ucd to
993  // consecutive numbering
994  for (unsigned int &vertex :
995  subcelldata.boundary_lines.back().vertices)
996  if (vertex_indices.find(vertex) != vertex_indices.end())
997  // vertex with this index exists
998  vertex = vertex_indices[vertex];
999  else
1000  {
1001  // no such vertex index
1002  AssertThrow(false, ExcInvalidVertexIndex(cell, vertex));
1004  }
1005  }
1006  else if ((cell_type == "quad") && (dim == 3))
1007  // boundary info
1008  {
1009  subcelldata.boundary_quads.emplace_back();
1010  in >> subcelldata.boundary_quads.back().vertices[0] >>
1011  subcelldata.boundary_quads.back().vertices[1] >>
1012  subcelldata.boundary_quads.back().vertices[2] >>
1013  subcelldata.boundary_quads.back().vertices[3];
1014 
1015  // to make sure that the cast won't fail
1017  ExcIndexRange(material_id,
1018  0,
1020  // we use only boundary_ids in the range from 0 to
1021  // numbers::internal_face_boundary_id-1
1023 
1024  // Make sure to set both manifold id and boundary id appropriately in
1025  // both cases:
1026  // numbers::internal_face_boundary_id and numbers::flat_manifold_id
1027  // are ignored in Triangulation::create_triangulation.
1028  if (apply_all_indicators_to_manifolds)
1029  {
1030  subcelldata.boundary_quads.back().boundary_id =
1032  subcelldata.boundary_quads.back().manifold_id =
1033  static_cast<types::manifold_id>(material_id);
1034  }
1035  else
1036  {
1037  subcelldata.boundary_quads.back().boundary_id =
1038  static_cast<types::boundary_id>(material_id);
1039  subcelldata.boundary_quads.back().manifold_id =
1041  }
1042 
1043  // transform from ucd to
1044  // consecutive numbering
1045  for (unsigned int &vertex :
1046  subcelldata.boundary_quads.back().vertices)
1047  if (vertex_indices.find(vertex) != vertex_indices.end())
1048  // vertex with this index exists
1049  vertex = vertex_indices[vertex];
1050  else
1051  {
1052  // no such vertex index
1053  Assert(false, ExcInvalidVertexIndex(cell, vertex));
1055  }
1056  }
1057  else
1058  // cannot read this
1059  AssertThrow(false, ExcUnknownIdentifier(cell_type));
1060  }
1061 
1062 
1063  // check that no forbidden arrays are used
1064  Assert(subcelldata.check_consistency(dim), ExcInternalError());
1065 
1066  AssertThrow(in, ExcIO());
1067 
1068  // do some clean-up on vertices...
1069  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1070  // ... and cells
1071  if (dim == spacedim)
1073  cells);
1075  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1076 }
1077 
1078 namespace
1079 {
1080  template <int dim, int spacedim>
1081  class Abaqus_to_UCD
1082  {
1083  public:
1084  Abaqus_to_UCD();
1085 
1086  void
1087  read_in_abaqus(std::istream &in);
1088  void
1089  write_out_avs_ucd(std::ostream &out) const;
1090 
1091  private:
1092  const double tolerance;
1093 
1094  std::vector<double>
1095  get_global_node_numbers(const int face_cell_no,
1096  const int face_cell_face_no) const;
1097 
1098  // NL: Stored as [ global node-id (int), x-coord, y-coord, z-coord ]
1099  std::vector<std::vector<double>> node_list;
1100  // CL: Stored as [ material-id (int), node1, node2, node3, node4, node5,
1101  // node6, node7, node8 ]
1102  std::vector<std::vector<double>> cell_list;
1103  // FL: Stored as [ sideset-id (int), node1, node2, node3, node4 ]
1104  std::vector<std::vector<double>> face_list;
1105  // ELSET: Stored as [ (std::string) elset_name = (std::vector) of cells
1106  // numbers]
1107  std::map<std::string, std::vector<int>> elsets_list;
1108  };
1109 } // namespace
1110 
1111 template <int dim, int spacedim>
1112 void
1114  const bool apply_all_indicators_to_manifolds)
1115 {
1116  Assert(tria != nullptr, ExcNoTriangulationSelected());
1117  // This implementation has only been verified for:
1118  // - 2d grids with codimension 0
1119  // - 3d grids with codimension 0
1120  // - 3d grids with codimension 1
1121  Assert((spacedim == 2 && dim == spacedim) ||
1122  (spacedim == 3 && (dim == spacedim || dim == spacedim - 1)),
1123  ExcNotImplemented());
1124  AssertThrow(in, ExcIO());
1125 
1126  // Read in the Abaqus file into an intermediate object
1127  // that is to be passed along to the UCD reader
1128  Abaqus_to_UCD<dim, spacedim> abaqus_to_ucd;
1129  abaqus_to_ucd.read_in_abaqus(in);
1130 
1131  std::stringstream in_ucd;
1132  abaqus_to_ucd.write_out_avs_ucd(in_ucd);
1133 
1134  // This next call is wrapped in a try-catch for the following reason:
1135  // It ensures that if the Abaqus mesh is read in correctly but produces
1136  // an erroneous result then the user is alerted to the source of the problem
1137  // and doesn't think that they've somehow called the wrong function.
1138  try
1139  {
1140  read_ucd(in_ucd, apply_all_indicators_to_manifolds);
1141  }
1142  catch (std::exception &exc)
1143  {
1144  std::cerr << "Exception on processing internal UCD data: " << std::endl
1145  << exc.what() << std::endl;
1146 
1147  AssertThrow(
1148  false,
1149  ExcMessage(
1150  "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1151  "More information is provided in an error message printed above. "
1152  "Are you sure that your ABAQUS mesh file conforms with the requirements "
1153  "listed in the documentation?"));
1154  }
1155  catch (...)
1156  {
1157  AssertThrow(
1158  false,
1159  ExcMessage(
1160  "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1161  "Are you sure that your ABAQUS mesh file conforms with the requirements "
1162  "listed in the documentation?"));
1163  }
1164 }
1165 
1166 
1167 template <int dim, int spacedim>
1168 void
1170 {
1171  Assert(tria != nullptr, ExcNoTriangulationSelected());
1172  Assert(dim == 2, ExcNotImplemented());
1173 
1174  AssertThrow(in, ExcIO());
1175 
1176  // skip comments at start of file
1177  skip_comment_lines(in, '#');
1178 
1179  // first read in identifier string
1180  std::string line;
1181  getline(in, line);
1182 
1183  AssertThrow(line == "MeshVersionFormatted 0", ExcInvalidDBMESHInput(line));
1184 
1185  skip_empty_lines(in);
1186 
1187  // next read dimension
1188  getline(in, line);
1189  AssertThrow(line == "Dimension", ExcInvalidDBMESHInput(line));
1190  unsigned int dimension;
1191  in >> dimension;
1192  AssertThrow(dimension == dim, ExcDBMESHWrongDimension(dimension));
1193  skip_empty_lines(in);
1194 
1195  // now there are a lot of fields of
1196  // which we don't know the exact
1197  // meaning and which are far from
1198  // being properly documented in the
1199  // manual. we skip everything until
1200  // we find a comment line with the
1201  // string "# END". at some point in
1202  // the future, someone may have the
1203  // knowledge to parse and interpret
1204  // the other fields in between as
1205  // well...
1206  while (getline(in, line), line.find("# END") == std::string::npos)
1207  ;
1208  skip_empty_lines(in);
1209 
1210 
1211  // now read vertices
1212  getline(in, line);
1213  AssertThrow(line == "Vertices", ExcInvalidDBMESHInput(line));
1214 
1215  unsigned int n_vertices;
1216  double dummy;
1217 
1218  in >> n_vertices;
1219  std::vector<Point<spacedim>> vertices(n_vertices);
1220  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1221  {
1222  // read vertex coordinates
1223  for (unsigned int d = 0; d < dim; ++d)
1224  in >> vertices[vertex][d];
1225  // read Ref phi_i, whatever that may be
1226  in >> dummy;
1227  }
1229 
1230  skip_empty_lines(in);
1231 
1232  // read edges. we ignore them at
1233  // present, so just read them and
1234  // discard the input
1235  getline(in, line);
1236  AssertThrow(line == "Edges", ExcInvalidDBMESHInput(line));
1237 
1238  unsigned int n_edges;
1239  in >> n_edges;
1240  for (unsigned int edge = 0; edge < n_edges; ++edge)
1241  {
1242  // read vertex indices
1243  in >> dummy >> dummy;
1244  // read Ref phi_i, whatever that may be
1245  in >> dummy;
1246  }
1248 
1249  skip_empty_lines(in);
1250 
1251 
1252 
1253  // read cracked edges (whatever
1254  // that may be). we ignore them at
1255  // present, so just read them and
1256  // discard the input
1257  getline(in, line);
1258  AssertThrow(line == "CrackedEdges", ExcInvalidDBMESHInput(line));
1259 
1260  in >> n_edges;
1261  for (unsigned int edge = 0; edge < n_edges; ++edge)
1262  {
1263  // read vertex indices
1264  in >> dummy >> dummy;
1265  // read Ref phi_i, whatever that may be
1266  in >> dummy;
1267  }
1269 
1270  skip_empty_lines(in);
1271 
1272 
1273  // now read cells.
1274  // set up array of cells
1275  getline(in, line);
1276  AssertThrow(line == "Quadrilaterals", ExcInvalidDBMESHInput(line));
1277 
1278  std::vector<CellData<dim>> cells;
1279  SubCellData subcelldata;
1280  unsigned int n_cells;
1281  in >> n_cells;
1282  for (unsigned int cell = 0; cell < n_cells; ++cell)
1283  {
1284  // read in vertex numbers. they
1285  // are 1-based, so subtract one
1286  cells.emplace_back();
1287  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1288  {
1289  in >> cells.back().vertices[i];
1290 
1291  AssertThrow((cells.back().vertices[i] >= 1) &&
1292  (static_cast<unsigned int>(cells.back().vertices[i]) <=
1293  vertices.size()),
1294  ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
1295 
1296  --cells.back().vertices[i];
1297  }
1298 
1299  // read and discard Ref phi_i
1300  in >> dummy;
1301  }
1303 
1304  skip_empty_lines(in);
1305 
1306 
1307  // then there are again a whole lot
1308  // of fields of which I have no
1309  // clue what they mean. skip them
1310  // all and leave the interpretation
1311  // to other implementors...
1312  while (getline(in, line), ((line.find("End") == std::string::npos) && (in)))
1313  ;
1314  // ok, so we are not at the end of
1315  // the file, that's it, mostly
1316 
1317 
1318  // check that no forbidden arrays are used
1319  Assert(subcelldata.check_consistency(dim), ExcInternalError());
1320 
1321  AssertThrow(in, ExcIO());
1322 
1323  // do some clean-up on vertices...
1324  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1325  // ...and cells
1327  cells);
1329  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1330 }
1331 
1332 
1333 
1334 template <int dim, int spacedim>
1335 void
1337 {
1338  Assert(false, ExcNotImplemented());
1339 }
1340 
1341 
1342 
1343 template <>
1344 void
1345 GridIn<2>::read_xda(std::istream &in)
1346 {
1347  Assert(tria != nullptr, ExcNoTriangulationSelected());
1348  AssertThrow(in, ExcIO());
1349 
1350  std::string line;
1351  // skip comments at start of file
1352  getline(in, line);
1353 
1354 
1355  unsigned int n_vertices;
1356  unsigned int n_cells;
1357 
1358  // read cells, throw away rest of line
1359  in >> n_cells;
1360  getline(in, line);
1361 
1362  in >> n_vertices;
1363  getline(in, line);
1364 
1365  // ignore following 8 lines
1366  for (unsigned int i = 0; i < 8; ++i)
1367  getline(in, line);
1368 
1369  // set up array of cells
1370  std::vector<CellData<2>> cells(n_cells);
1371  SubCellData subcelldata;
1372 
1373  for (unsigned int cell = 0; cell < n_cells; ++cell)
1374  {
1375  // note that since in the input
1376  // file we found the number of
1377  // cells at the top, there
1378  // should still be input here,
1379  // so check this:
1380  AssertThrow(in, ExcIO());
1382 
1383  for (unsigned int &vertex : cells[cell].vertices)
1384  in >> vertex;
1385  }
1386 
1387 
1388 
1389  // set up array of vertices
1390  std::vector<Point<2>> vertices(n_vertices);
1391  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1392  {
1393  double x[3];
1394 
1395  // read vertex
1396  in >> x[0] >> x[1] >> x[2];
1397 
1398  // store vertex
1399  for (unsigned int d = 0; d < 2; ++d)
1400  vertices[vertex](d) = x[d];
1401  }
1402  AssertThrow(in, ExcIO());
1403 
1404  // do some clean-up on vertices...
1405  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1406  // ... and cells
1409  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1410 }
1411 
1412 
1413 
1414 template <>
1415 void
1416 GridIn<3>::read_xda(std::istream &in)
1417 {
1418  Assert(tria != nullptr, ExcNoTriangulationSelected());
1419  AssertThrow(in, ExcIO());
1420 
1421  static const unsigned int xda_to_dealII_map[] = {0, 1, 5, 4, 3, 2, 6, 7};
1422 
1423  std::string line;
1424  // skip comments at start of file
1425  getline(in, line);
1426 
1427 
1428  unsigned int n_vertices;
1429  unsigned int n_cells;
1430 
1431  // read cells, throw away rest of line
1432  in >> n_cells;
1433  getline(in, line);
1434 
1435  in >> n_vertices;
1436  getline(in, line);
1437 
1438  // ignore following 8 lines
1439  for (unsigned int i = 0; i < 8; ++i)
1440  getline(in, line);
1441 
1442  // set up array of cells
1443  std::vector<CellData<3>> cells(n_cells);
1444  SubCellData subcelldata;
1445 
1446  for (unsigned int cell = 0; cell < n_cells; ++cell)
1447  {
1448  // note that since in the input
1449  // file we found the number of
1450  // cells at the top, there
1451  // should still be input here,
1452  // so check this:
1453  AssertThrow(in, ExcIO());
1455 
1456  unsigned int xda_ordered_nodes[8];
1457 
1458  for (unsigned int &xda_ordered_node : xda_ordered_nodes)
1459  in >> xda_ordered_node;
1460 
1461  for (unsigned int i = 0; i < 8; i++)
1462  cells[cell].vertices[i] = xda_ordered_nodes[xda_to_dealII_map[i]];
1463  }
1464 
1465 
1466 
1467  // set up array of vertices
1468  std::vector<Point<3>> vertices(n_vertices);
1469  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1470  {
1471  double x[3];
1472 
1473  // read vertex
1474  in >> x[0] >> x[1] >> x[2];
1475 
1476  // store vertex
1477  for (unsigned int d = 0; d < 3; ++d)
1478  vertices[vertex](d) = x[d];
1479  }
1480  AssertThrow(in, ExcIO());
1481 
1482  // do some clean-up on vertices...
1483  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1484  // ... and cells
1487  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1488 }
1489 
1490 
1491 
1492 template <int dim, int spacedim>
1493 void
1495 {
1496  Assert(tria != nullptr, ExcNoTriangulationSelected());
1497  AssertThrow(in, ExcIO());
1498 
1499  unsigned int n_vertices;
1500  unsigned int n_cells;
1501  unsigned int dummy;
1502  std::string line;
1503  // This array stores maps from the 'entities' to the 'physical tags' for
1504  // points, curves, surfaces and volumes. We use this information later to
1505  // assign boundary ids.
1506  std::array<std::map<int, int>, 4> tag_maps;
1507 
1508  in >> line;
1509 
1510  // first determine file format
1511  unsigned int gmsh_file_format = 0;
1512  if (line == "@f$NOD")
1513  gmsh_file_format = 10;
1514  else if (line == "@f$MeshFormat")
1515  gmsh_file_format = 20;
1516  else
1517  AssertThrow(false, ExcInvalidGMSHInput(line));
1518 
1519  // if file format is 2.0 or greater then we also have to read the rest of the
1520  // header
1521  if (gmsh_file_format == 20)
1522  {
1523  double version;
1524  unsigned int file_type, data_size;
1525 
1526  in >> version >> file_type >> data_size;
1527 
1528  Assert((version >= 2.0) && (version <= 4.1), ExcNotImplemented());
1529  gmsh_file_format = static_cast<unsigned int>(version * 10);
1530 
1531  Assert(file_type == 0, ExcNotImplemented());
1532  Assert(data_size == sizeof(double), ExcNotImplemented());
1533 
1534  // read the end of the header and the first line of the nodes description
1535  // to synch ourselves with the format 1 handling above
1536  in >> line;
1537  AssertThrow(line == "@f$EndMeshFormat", ExcInvalidGMSHInput(line));
1538 
1539  in >> line;
1540  // if the next block is of kind @f$PhysicalNames, ignore it
1541  if (line == "@f$PhysicalNames")
1542  {
1543  do
1544  {
1545  in >> line;
1546  }
1547  while (line != "@f$EndPhysicalNames");
1548  in >> line;
1549  }
1550 
1551  // if the next block is of kind @f$Entities, parse it
1552  if (line == "@f$Entities")
1553  {
1554  unsigned long n_points, n_curves, n_surfaces, n_volumes;
1555 
1556  in >> n_points >> n_curves >> n_surfaces >> n_volumes;
1557  for (unsigned int i = 0; i < n_points; ++i)
1558  {
1559  // parse point ids
1560  int tag;
1561  unsigned int n_physicals;
1562  double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1563  box_max_z;
1564 
1565  // we only care for 'tag' as key for tag_maps[0]
1566  if (gmsh_file_format > 40)
1567  {
1568  in >> tag >> box_min_x >> box_min_y >> box_min_z >>
1569  n_physicals;
1570  box_max_x = box_min_x;
1571  box_max_y = box_min_y;
1572  box_max_z = box_min_z;
1573  }
1574  else
1575  {
1576  in >> tag >> box_min_x >> box_min_y >> box_min_z >>
1577  box_max_x >> box_max_y >> box_max_z >> n_physicals;
1578  }
1579  // if there is a physical tag, we will use it as boundary id below
1580  AssertThrow(n_physicals < 2,
1581  ExcMessage("More than one tag is not supported!"));
1582  // if there is no physical tag, use 0 as default
1583  int physical_tag = 0;
1584  for (unsigned int j = 0; j < n_physicals; ++j)
1585  in >> physical_tag;
1586  tag_maps[0][tag] = physical_tag;
1587  }
1588  for (unsigned int i = 0; i < n_curves; ++i)
1589  {
1590  // parse curve ids
1591  int tag;
1592  unsigned int n_physicals;
1593  double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1594  box_max_z;
1595 
1596  // we only care for 'tag' as key for tag_maps[1]
1597  in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1598  box_max_y >> box_max_z >> n_physicals;
1599  // if there is a physical tag, we will use it as boundary id below
1600  AssertThrow(n_physicals < 2,
1601  ExcMessage("More than one tag is not supported!"));
1602  // if there is no physical tag, use 0 as default
1603  int physical_tag = 0;
1604  for (unsigned int j = 0; j < n_physicals; ++j)
1605  in >> physical_tag;
1606  tag_maps[1][tag] = physical_tag;
1607  // we don't care about the points associated to a curve, but have
1608  // to parse them anyway because their format is unstructured
1609  in >> n_points;
1610  for (unsigned int j = 0; j < n_points; ++j)
1611  in >> tag;
1612  }
1613 
1614  for (unsigned int i = 0; i < n_surfaces; ++i)
1615  {
1616  // parse surface ids
1617  int tag;
1618  unsigned int n_physicals;
1619  double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1620  box_max_z;
1621 
1622  // we only care for 'tag' as key for tag_maps[2]
1623  in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1624  box_max_y >> box_max_z >> n_physicals;
1625  // if there is a physical tag, we will use it as boundary id below
1626  AssertThrow(n_physicals < 2,
1627  ExcMessage("More than one tag is not supported!"));
1628  // if there is no physical tag, use 0 as default
1629  int physical_tag = 0;
1630  for (unsigned int j = 0; j < n_physicals; ++j)
1631  in >> physical_tag;
1632  tag_maps[2][tag] = physical_tag;
1633  // we don't care about the curves associated to a surface, but
1634  // have to parse them anyway because their format is unstructured
1635  in >> n_curves;
1636  for (unsigned int j = 0; j < n_curves; ++j)
1637  in >> tag;
1638  }
1639  for (unsigned int i = 0; i < n_volumes; ++i)
1640  {
1641  // parse volume ids
1642  int tag;
1643  unsigned int n_physicals;
1644  double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1645  box_max_z;
1646 
1647  // we only care for 'tag' as key for tag_maps[3]
1648  in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1649  box_max_y >> box_max_z >> n_physicals;
1650  // if there is a physical tag, we will use it as boundary id below
1651  AssertThrow(n_physicals < 2,
1652  ExcMessage("More than one tag is not supported!"));
1653  // if there is no physical tag, use 0 as default
1654  int physical_tag = 0;
1655  for (unsigned int j = 0; j < n_physicals; ++j)
1656  in >> physical_tag;
1657  tag_maps[3][tag] = physical_tag;
1658  // we don't care about the surfaces associated to a volume, but
1659  // have to parse them anyway because their format is unstructured
1660  in >> n_surfaces;
1661  for (unsigned int j = 0; j < n_surfaces; ++j)
1662  in >> tag;
1663  }
1664  in >> line;
1665  AssertThrow(line == "@f$EndEntities", ExcInvalidGMSHInput(line));
1666  in >> line;
1667  }
1668 
1669  // if the next block is of kind @f$PartitionedEntities, ignore it
1670  if (line == "@f$PartitionedEntities")
1671  {
1672  do
1673  {
1674  in >> line;
1675  }
1676  while (line != "@f$EndPartitionedEntities");
1677  in >> line;
1678  }
1679 
1680  // but the next thing should,
1681  // in any case, be the list of
1682  // nodes:
1683  AssertThrow(line == "@f$Nodes", ExcInvalidGMSHInput(line));
1684  }
1685 
1686  // now read the nodes list
1687  int n_entity_blocks = 1;
1688  if (gmsh_file_format > 40)
1689  {
1690  int min_node_tag;
1691  int max_node_tag;
1692  in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag;
1693  }
1694  else if (gmsh_file_format == 40)
1695  {
1696  in >> n_entity_blocks >> n_vertices;
1697  }
1698  else
1699  in >> n_vertices;
1700  std::vector<Point<spacedim>> vertices(n_vertices);
1701  // set up mapping between numbering
1702  // in msh-file (nod) and in the
1703  // vertices vector
1704  std::map<int, int> vertex_indices;
1705 
1706  {
1707  unsigned int global_vertex = 0;
1708  for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
1709  {
1710  int parametric;
1711  unsigned long numNodes;
1712 
1713  if (gmsh_file_format < 40)
1714  {
1715  numNodes = n_vertices;
1716  parametric = 0;
1717  }
1718  else
1719  {
1720  // for gmsh_file_format 4.1 the order of tag and dim is reversed,
1721  // but we are ignoring both anyway.
1722  int tagEntity, dimEntity;
1723  in >> tagEntity >> dimEntity >> parametric >> numNodes;
1724  }
1725 
1726  std::vector<int> vertex_numbers;
1727  int vertex_number;
1728  if (gmsh_file_format > 40)
1729  for (unsigned long vertex_per_entity = 0;
1730  vertex_per_entity < numNodes;
1731  ++vertex_per_entity)
1732  {
1733  in >> vertex_number;
1734  vertex_numbers.push_back(vertex_number);
1735  }
1736 
1737  for (unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes;
1738  ++vertex_per_entity, ++global_vertex)
1739  {
1740  int vertex_number;
1741  double x[3];
1742 
1743  // read vertex
1744  if (gmsh_file_format > 40)
1745  {
1746  vertex_number = vertex_numbers[vertex_per_entity];
1747  in >> x[0] >> x[1] >> x[2];
1748  }
1749  else
1750  in >> vertex_number >> x[0] >> x[1] >> x[2];
1751 
1752  for (unsigned int d = 0; d < spacedim; ++d)
1753  vertices[global_vertex](d) = x[d];
1754  // store mapping
1755  vertex_indices[vertex_number] = global_vertex;
1756 
1757  // ignore parametric coordinates
1758  if (parametric != 0)
1759  {
1760  double u = 0.;
1761  double v = 0.;
1762  in >> u >> v;
1763  (void)u;
1764  (void)v;
1765  }
1766  }
1767  }
1768  AssertDimension(global_vertex, n_vertices);
1769  }
1770 
1771  // Assert we reached the end of the block
1772  in >> line;
1773  static const std::string end_nodes_marker[] = {"@f$ENDNOD", "@f$EndNodes"};
1774  AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1],
1775  ExcInvalidGMSHInput(line));
1776 
1777  // Now read in next bit
1778  in >> line;
1779  static const std::string begin_elements_marker[] = {"@f$ELM", "@f$Elements"};
1780  AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1],
1781  ExcInvalidGMSHInput(line));
1782 
1783  // now read the cell list
1784  if (gmsh_file_format > 40)
1785  {
1786  int min_node_tag;
1787  int max_node_tag;
1788  in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag;
1789  }
1790  else if (gmsh_file_format == 40)
1791  {
1792  in >> n_entity_blocks >> n_cells;
1793  }
1794  else
1795  {
1796  n_entity_blocks = 1;
1797  in >> n_cells;
1798  }
1799 
1800  // set up array of cells and subcells (faces). In 1d, there is currently no
1801  // standard way in deal.II to pass boundary indicators attached to individual
1802  // vertices, so do this by hand via the boundary_ids_1d array
1803  std::vector<CellData<dim>> cells;
1804  SubCellData subcelldata;
1805  std::map<unsigned int, types::boundary_id> boundary_ids_1d;
1806  bool is_quad_or_hex_mesh = false;
1807  bool is_tria_or_tet_mesh = false;
1808 
1809  {
1810  unsigned int global_cell = 0;
1811  for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
1812  {
1813  unsigned int material_id;
1814  unsigned long numElements;
1815  int cell_type;
1816 
1817  if (gmsh_file_format < 40)
1818  {
1819  material_id = 0;
1820  cell_type = 0;
1821  numElements = n_cells;
1822  }
1823  else if (gmsh_file_format == 40)
1824  {
1825  int tagEntity, dimEntity;
1826  in >> tagEntity >> dimEntity >> cell_type >> numElements;
1827  material_id = tag_maps[dimEntity][tagEntity];
1828  }
1829  else
1830  {
1831  // for gmsh_file_format 4.1 the order of tag and dim is reversed,
1832  int tagEntity, dimEntity;
1833  in >> dimEntity >> tagEntity >> cell_type >> numElements;
1834  material_id = tag_maps[dimEntity][tagEntity];
1835  }
1836 
1837  for (unsigned int cell_per_entity = 0; cell_per_entity < numElements;
1838  ++cell_per_entity, ++global_cell)
1839  {
1840  // note that since in the input
1841  // file we found the number of
1842  // cells at the top, there
1843  // should still be input here,
1844  // so check this:
1845  AssertThrow(in, ExcIO());
1846 
1847  unsigned int nod_num;
1848 
1849  /*
1850  For file format version 1, the format of each cell is as follows:
1851  elm-number elm-type reg-phys reg-elem number-of-nodes
1852  node-number-list
1853 
1854  However, for version 2, the format reads like this:
1855  elm-number elm-type number-of-tags < tag > ... node-number-list
1856 
1857  For version 4, we have:
1858  tag(int) numVert(int) ...
1859 
1860  In the following, we will ignore the element number (we simply
1861  enumerate them in the order in which we read them, and we will
1862  take reg-phys (version 1) or the first tag (version 2, if any tag
1863  is given at all) as material id. For version 4, we already read
1864  the material and the cell type in above.
1865  */
1866 
1867  unsigned int elm_number = 0;
1868  if (gmsh_file_format < 40)
1869  {
1870  in >> elm_number // ELM-NUMBER
1871  >> cell_type; // ELM-TYPE
1872  }
1873 
1874  if (gmsh_file_format < 20)
1875  {
1876  in >> material_id // REG-PHYS
1877  >> dummy // reg_elm
1878  >> nod_num;
1879  }
1880  else if (gmsh_file_format < 40)
1881  {
1882  // read the tags; ignore all but the first one which we will
1883  // interpret as the material_id (for cells) or boundary_id
1884  // (for faces)
1885  unsigned int n_tags;
1886  in >> n_tags;
1887  if (n_tags > 0)
1888  in >> material_id;
1889  else
1890  material_id = 0;
1891 
1892  for (unsigned int i = 1; i < n_tags; ++i)
1893  in >> dummy;
1894 
1895  if (cell_type == 1) // line
1896  nod_num = 2;
1897  else if (cell_type == 2) // tri
1898  nod_num = 3;
1899  else if (cell_type == 3) // quad
1900  nod_num = 4;
1901  else if (cell_type == 4) // tet
1902  nod_num = 4;
1903  else if (cell_type == 5) // hex
1904  nod_num = 8;
1905  }
1906  else
1907  {
1908  // ignore tag
1909  int tag;
1910  in >> tag;
1911 
1912  if (cell_type == 1) // line
1913  nod_num = 2;
1914  else if (cell_type == 2) // tri
1915  nod_num = 3;
1916  else if (cell_type == 3) // quad
1917  nod_num = 4;
1918  else if (cell_type == 4) // tet
1919  nod_num = 4;
1920  else if (cell_type == 5) // hex
1921  nod_num = 8;
1922  }
1923 
1924 
1925  /* `ELM-TYPE'
1926  defines the geometrical type of the N-th element:
1927  `1'
1928  Line (2 nodes, 1 edge).
1929 
1930  `2'
1931  Triangle (3 nodes, 3 edges).
1932 
1933  `3'
1934  Quadrangle (4 nodes, 4 edges).
1935 
1936  `4'
1937  Tetrahedron (4 nodes, 6 edges, 6 faces).
1938 
1939  `5'
1940  Hexahedron (8 nodes, 12 edges, 6 faces).
1941 
1942  `15'
1943  Point (1 node).
1944  */
1945 
1946  if (((cell_type == 1) && (dim == 1)) ||
1947  ((cell_type == 2) && (dim == 2)) ||
1948  ((cell_type == 3) && (dim == 2)) ||
1949  ((cell_type == 4) && (dim == 3)) ||
1950  ((cell_type == 5) && (dim == 3)))
1951  // found a cell
1952  {
1953  unsigned int vertices_per_cell = 0;
1954  if (cell_type == 1) // line
1955  vertices_per_cell = 2;
1956  else if (cell_type == 2) // tri
1957  {
1958  vertices_per_cell = 3;
1959  is_tria_or_tet_mesh = true;
1960  }
1961  else if (cell_type == 3) // quad
1962  {
1963  vertices_per_cell = 4;
1964  is_quad_or_hex_mesh = true;
1965  }
1966  else if (cell_type == 4) // tet
1967  {
1968  vertices_per_cell = 4;
1969  is_tria_or_tet_mesh = true;
1970  }
1971  else if (cell_type == 5) // hex
1972  {
1973  vertices_per_cell = 8;
1974  is_quad_or_hex_mesh = true;
1975  }
1976 
1977  AssertThrow(nod_num == vertices_per_cell,
1978  ExcMessage(
1979  "Number of nodes does not coincide with the "
1980  "number required for this object"));
1981 
1982  // allocate and read indices
1983  cells.emplace_back();
1984  cells.back().vertices.resize(vertices_per_cell);
1985  for (unsigned int i = 0; i < vertices_per_cell; ++i)
1986  in >> cells.back().vertices[i];
1987 
1988  // to make sure that the cast won't fail
1989  Assert(material_id <=
1991  ExcIndexRange(
1992  material_id,
1993  0,
1995  // we use only material_ids in the range from 0 to
1996  // numbers::invalid_material_id-1
1998 
1999  cells.back().material_id = material_id;
2000 
2001  // transform from ucd to
2002  // consecutive numbering
2003  for (unsigned int i = 0; i < vertices_per_cell; ++i)
2004  {
2005  AssertThrow(
2006  vertex_indices.find(cells.back().vertices[i]) !=
2007  vertex_indices.end(),
2008  ExcInvalidVertexIndexGmsh(cell_per_entity,
2009  elm_number,
2010  cells.back().vertices[i]));
2011 
2012  // vertex with this index exists
2013  cells.back().vertices[i] =
2014  vertex_indices[cells.back().vertices[i]];
2015  }
2016  }
2017  else if ((cell_type == 1) && ((dim == 2) || (dim == 3)))
2018  // boundary info
2019  {
2020  subcelldata.boundary_lines.emplace_back();
2021  in >> subcelldata.boundary_lines.back().vertices[0] >>
2022  subcelldata.boundary_lines.back().vertices[1];
2023 
2024  // to make sure that the cast won't fail
2025  Assert(material_id <=
2027  ExcIndexRange(
2028  material_id,
2029  0,
2031  // we use only boundary_ids in the range from 0 to
2032  // numbers::internal_face_boundary_id-1
2033  AssertIndexRange(material_id,
2035 
2036  subcelldata.boundary_lines.back().boundary_id =
2037  static_cast<types::boundary_id>(material_id);
2038 
2039  // transform from ucd to
2040  // consecutive numbering
2041  for (unsigned int &vertex :
2042  subcelldata.boundary_lines.back().vertices)
2043  if (vertex_indices.find(vertex) != vertex_indices.end())
2044  // vertex with this index exists
2045  vertex = vertex_indices[vertex];
2046  else
2047  {
2048  // no such vertex index
2049  AssertThrow(false,
2050  ExcInvalidVertexIndex(cell_per_entity,
2051  vertex));
2053  }
2054  }
2055  else if ((cell_type == 2 || cell_type == 3) && (dim == 3))
2056  // boundary info
2057  {
2058  unsigned int vertices_per_cell = 0;
2059  // check cell type
2060  if (cell_type == 2) // tri
2061  {
2062  vertices_per_cell = 3;
2063  is_tria_or_tet_mesh = true;
2064  }
2065  else if (cell_type == 3) // quad
2066  {
2067  vertices_per_cell = 4;
2068  is_quad_or_hex_mesh = true;
2069  }
2070 
2071  subcelldata.boundary_quads.emplace_back();
2072 
2073  // resize vertices
2074  subcelldata.boundary_quads.back().vertices.resize(
2075  vertices_per_cell);
2076  // for loop
2077  for (unsigned int i = 0; i < vertices_per_cell; ++i)
2078  in >> subcelldata.boundary_quads.back().vertices[i];
2079 
2080  // to make sure that the cast won't fail
2081  Assert(material_id <=
2083  ExcIndexRange(
2084  material_id,
2085  0,
2087  // we use only boundary_ids in the range from 0 to
2088  // numbers::internal_face_boundary_id-1
2089  AssertIndexRange(material_id,
2091 
2092  subcelldata.boundary_quads.back().boundary_id =
2093  static_cast<types::boundary_id>(material_id);
2094 
2095  // transform from gmsh to
2096  // consecutive numbering
2097  for (unsigned int &vertex :
2098  subcelldata.boundary_quads.back().vertices)
2099  if (vertex_indices.find(vertex) != vertex_indices.end())
2100  // vertex with this index exists
2101  vertex = vertex_indices[vertex];
2102  else
2103  {
2104  // no such vertex index
2105  Assert(false,
2106  ExcInvalidVertexIndex(cell_per_entity, vertex));
2108  }
2109  }
2110  else if (cell_type == 15)
2111  {
2112  // read the indices of nodes given
2113  unsigned int node_index = 0;
2114  if (gmsh_file_format < 20)
2115  {
2116  // For points (cell_type==15), we can only ever
2117  // list one node index.
2118  AssertThrow(nod_num == 1, ExcInternalError());
2119  in >> node_index;
2120  }
2121  else
2122  {
2123  in >> node_index;
2124  }
2125 
2126  // we only care about boundary indicators assigned to individual
2127  // vertices in 1d (because otherwise the vertices are not faces)
2128  if (dim == 1)
2129  boundary_ids_1d[vertex_indices[node_index]] = material_id;
2130  }
2131  else
2132  {
2133  AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type));
2134  }
2135  }
2136  }
2137  AssertDimension(global_cell, n_cells);
2138  }
2139  // Assert we reached the end of the block
2140  in >> line;
2141  static const std::string end_elements_marker[] = {"@f$ENDELM", "@f$EndElements"};
2142  AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2143  ExcInvalidGMSHInput(line));
2144 
2145  // check that no forbidden arrays are used
2146  Assert(subcelldata.check_consistency(dim), ExcInternalError());
2147 
2148  AssertThrow(in, ExcIO());
2149 
2150  // check that we actually read some cells.
2151  AssertThrow(cells.size() > 0, ExcGmshNoCellInformation());
2152 
2153  // make sure that only either simplex or hypercube cells are available
2154  //
2155  // TODO: the functions below (GridTools::delete_unused_vertices(),
2156  // GridTools::invert_all_cells_of_negative_grid(),
2157  // GridReordering::reorder_cells(),
2158  // Triangulation::create_triangulation_compatibility()) need to be revisited
2159  // for simplex meshes
2160  AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
2161  ExcNotImplemented());
2162 
2163  if (dim == 1 || is_quad_or_hex_mesh)
2164  {
2165  // do some clean-up on vertices...
2166  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2167  // ... and cells
2168  if (dim == spacedim)
2170  vertices, cells);
2172  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2173  }
2174  else
2175  {
2176  tria->create_triangulation(vertices, cells, subcelldata);
2177  }
2178 
2179  // in 1d, we also have to attach boundary ids to vertices, which does not
2180  // currently work through the call above
2181  if (dim == 1)
2182  assign_1d_boundary_ids(boundary_ids_1d, *tria);
2183 }
2184 
2185 
2186 
2187 template <int dim, int spacedim>
2188 void
2190  std::string & header,
2191  std::vector<unsigned int> &tecplot2deal,
2192  unsigned int & n_vars,
2193  unsigned int & n_vertices,
2194  unsigned int & n_cells,
2195  std::vector<unsigned int> &IJK,
2196  bool & structured,
2197  bool & blocked)
2198 {
2199  Assert(tecplot2deal.size() == dim, ExcInternalError());
2200  Assert(IJK.size() == dim, ExcInternalError());
2201  // initialize the output variables
2202  n_vars = 0;
2203  n_vertices = 0;
2204  n_cells = 0;
2205  switch (dim)
2206  {
2207  case 3:
2208  IJK[2] = 0;
2210  case 2:
2211  IJK[1] = 0;
2213  case 1:
2214  IJK[0] = 0;
2215  }
2216  structured = true;
2217  blocked = false;
2218 
2219  // convert the string to upper case
2220  std::transform(header.begin(), header.end(), header.begin(), ::toupper);
2221 
2222  // replace all tabs, commas, newlines by
2223  // whitespaces
2224  std::replace(header.begin(), header.end(), '\t', ' ');
2225  std::replace(header.begin(), header.end(), ',', ' ');
2226  std::replace(header.begin(), header.end(), '\n', ' ');
2227 
2228  // now remove whitespace in front of and
2229  // after '='
2230  std::string::size_type pos = header.find('=');
2231 
2232  while (pos != static_cast<std::string::size_type>(std::string::npos))
2233  if (header[pos + 1] == ' ')
2234  header.erase(pos + 1, 1);
2235  else if (header[pos - 1] == ' ')
2236  {
2237  header.erase(pos - 1, 1);
2238  --pos;
2239  }
2240  else
2241  pos = header.find('=', ++pos);
2242 
2243  // split the string into individual entries
2244  std::vector<std::string> entries =
2245  Utilities::break_text_into_lines(header, 1, ' ');
2246 
2247  // now go through the list and try to extract
2248  for (unsigned int i = 0; i < entries.size(); ++i)
2249  {
2250  if (Utilities::match_at_string_start(entries[i], "VARIABLES=\""))
2251  {
2252  ++n_vars;
2253  // we assume, that the first variable
2254  // is x or no coordinate at all (not y or z)
2255  if (Utilities::match_at_string_start(entries[i], "VARIABLES=\"X\""))
2256  {
2257  tecplot2deal[0] = 0;
2258  }
2259  ++i;
2260  while (entries[i][0] == '"')
2261  {
2262  if (entries[i] == "\"X\"")
2263  tecplot2deal[0] = n_vars;
2264  else if (entries[i] == "\"Y\"")
2265  {
2266  // we assume, that y contains
2267  // zero data in 1d, so do
2268  // nothing
2269  if (dim > 1)
2270  tecplot2deal[1] = n_vars;
2271  }
2272  else if (entries[i] == "\"Z\"")
2273  {
2274  // we assume, that z contains
2275  // zero data in 1d and 2d, so
2276  // do nothing
2277  if (dim > 2)
2278  tecplot2deal[2] = n_vars;
2279  }
2280  ++n_vars;
2281  ++i;
2282  }
2283  // set i back, so that the next
2284  // string is treated correctly
2285  --i;
2286 
2287  AssertThrow(
2288  n_vars >= dim,
2289  ExcMessage(
2290  "Tecplot file must contain at least one variable for each dimension"));
2291  for (unsigned int d = 1; d < dim; ++d)
2292  AssertThrow(
2293  tecplot2deal[d] > 0,
2294  ExcMessage(
2295  "Tecplot file must contain at least one variable for each dimension."));
2296  }
2297  else if (Utilities::match_at_string_start(entries[i], "ZONETYPE=ORDERED"))
2298  structured = true;
2299  else if (Utilities::match_at_string_start(entries[i],
2300  "ZONETYPE=FELINESEG") &&
2301  dim == 1)
2302  structured = false;
2303  else if (Utilities::match_at_string_start(entries[i],
2304  "ZONETYPE=FEQUADRILATERAL") &&
2305  dim == 2)
2306  structured = false;
2307  else if (Utilities::match_at_string_start(entries[i],
2308  "ZONETYPE=FEBRICK") &&
2309  dim == 3)
2310  structured = false;
2311  else if (Utilities::match_at_string_start(entries[i], "ZONETYPE="))
2312  // unsupported ZONETYPE
2313  {
2314  AssertThrow(false,
2315  ExcMessage(
2316  "The tecplot file contains an unsupported ZONETYPE."));
2317  }
2318  else if (Utilities::match_at_string_start(entries[i],
2319  "DATAPACKING=POINT"))
2320  blocked = false;
2321  else if (Utilities::match_at_string_start(entries[i],
2322  "DATAPACKING=BLOCK"))
2323  blocked = true;
2324  else if (Utilities::match_at_string_start(entries[i], "F=POINT"))
2325  {
2326  structured = true;
2327  blocked = false;
2328  }
2329  else if (Utilities::match_at_string_start(entries[i], "F=BLOCK"))
2330  {
2331  structured = true;
2332  blocked = true;
2333  }
2334  else if (Utilities::match_at_string_start(entries[i], "F=FEPOINT"))
2335  {
2336  structured = false;
2337  blocked = false;
2338  }
2339  else if (Utilities::match_at_string_start(entries[i], "F=FEBLOCK"))
2340  {
2341  structured = false;
2342  blocked = true;
2343  }
2344  else if (Utilities::match_at_string_start(entries[i],
2345  "ET=QUADRILATERAL") &&
2346  dim == 2)
2347  structured = false;
2348  else if (Utilities::match_at_string_start(entries[i], "ET=BRICK") &&
2349  dim == 3)
2350  structured = false;
2351  else if (Utilities::match_at_string_start(entries[i], "ET="))
2352  // unsupported ElementType
2353  {
2354  AssertThrow(
2355  false,
2356  ExcMessage(
2357  "The tecplot file contains an unsupported ElementType."));
2358  }
2359  else if (Utilities::match_at_string_start(entries[i], "I="))
2360  IJK[0] = Utilities::get_integer_at_position(entries[i], 2).first;
2361  else if (Utilities::match_at_string_start(entries[i], "J="))
2362  {
2363  IJK[1] = Utilities::get_integer_at_position(entries[i], 2).first;
2364  AssertThrow(
2365  dim > 1 || IJK[1] == 1,
2366  ExcMessage(
2367  "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
2368  }
2369  else if (Utilities::match_at_string_start(entries[i], "K="))
2370  {
2371  IJK[2] = Utilities::get_integer_at_position(entries[i], 2).first;
2372  AssertThrow(
2373  dim > 2 || IJK[2] == 1,
2374  ExcMessage(
2375  "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
2376  }
2377  else if (Utilities::match_at_string_start(entries[i], "N="))
2378  n_vertices = Utilities::get_integer_at_position(entries[i], 2).first;
2379  else if (Utilities::match_at_string_start(entries[i], "E="))
2380  n_cells = Utilities::get_integer_at_position(entries[i], 2).first;
2381  }
2382 
2383  // now we have read all the fields we are
2384  // interested in. do some checks and
2385  // calculate the variables
2386  if (structured)
2387  {
2388  n_vertices = 1;
2389  n_cells = 1;
2390  for (unsigned int d = 0; d < dim; ++d)
2391  {
2392  AssertThrow(
2393  IJK[d] > 0,
2394  ExcMessage(
2395  "Tecplot file does not contain a complete and consistent set of parameters"));
2396  n_vertices *= IJK[d];
2397  n_cells *= (IJK[d] - 1);
2398  }
2399  }
2400  else
2401  {
2402  AssertThrow(
2403  n_vertices > 0,
2404  ExcMessage(
2405  "Tecplot file does not contain a complete and consistent set of parameters"));
2406  if (n_cells == 0)
2407  // this means an error, although
2408  // tecplot itself accepts entries like
2409  // 'J=20' instead of 'E=20'. therefore,
2410  // take the max of IJK
2411  n_cells = *std::max_element(IJK.begin(), IJK.end());
2412  AssertThrow(
2413  n_cells > 0,
2414  ExcMessage(
2415  "Tecplot file does not contain a complete and consistent set of parameters"));
2416  }
2417 }
2418 
2419 
2420 
2421 template <>
2422 void
2423 GridIn<2>::read_tecplot(std::istream &in)
2424 {
2425  const unsigned int dim = 2;
2426  const unsigned int spacedim = 2;
2427  Assert(tria != nullptr, ExcNoTriangulationSelected());
2428  AssertThrow(in, ExcIO());
2429 
2430  // skip comments at start of file
2431  skip_comment_lines(in, '#');
2432 
2433  // some strings for parsing the header
2434  std::string line, header;
2435 
2436  // first, concatenate all header lines
2437  // create a searchstring with almost all
2438  // letters. exclude e and E from the letters
2439  // to search, as they might appear in
2440  // exponential notation
2441  std::string letters = "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
2442 
2443  getline(in, line);
2444  while (line.find_first_of(letters) != std::string::npos)
2445  {
2446  header += " " + line;
2447  getline(in, line);
2448  }
2449 
2450  // now create some variables holding
2451  // important information on the mesh, get
2452  // this information from the header string
2453  std::vector<unsigned int> tecplot2deal(dim);
2454  std::vector<unsigned int> IJK(dim);
2455  unsigned int n_vars, n_vertices, n_cells;
2456  bool structured, blocked;
2457 
2458  parse_tecplot_header(header,
2459  tecplot2deal,
2460  n_vars,
2461  n_vertices,
2462  n_cells,
2463  IJK,
2464  structured,
2465  blocked);
2466 
2467  // reserve space for vertices. note, that in
2468  // tecplot vertices are ordered beginning
2469  // with 1, whereas in deal all indices start
2470  // with 0. in order not to use -1 for all the
2471  // connectivity information, a 0th vertex
2472  // (unused) is inserted at the origin.
2473  std::vector<Point<spacedim>> vertices(n_vertices + 1);
2474  vertices[0] = Point<spacedim>();
2475  // reserve space for cells
2476  std::vector<CellData<dim>> cells(n_cells);
2477  SubCellData subcelldata;
2478 
2479  if (blocked)
2480  {
2481  // blocked data format. first we get all
2482  // the values of the first variable for
2483  // all points, after that we get all
2484  // values for the second variable and so
2485  // on.
2486 
2487  // dummy variable to read in all the info
2488  // we do not want to use
2489  double dummy;
2490  // which is the first index to read in
2491  // the loop (see below)
2492  unsigned int next_index = 0;
2493 
2494  // note, that we have already read the
2495  // first line containing the first variable
2496  if (tecplot2deal[0] == 0)
2497  {
2498  // we need the information in this
2499  // line, so extract it
2500  std::vector<std::string> first_var =
2502  char *endptr;
2503  for (unsigned int i = 1; i < first_var.size() + 1; ++i)
2504  vertices[i](0) = std::strtod(first_var[i - 1].c_str(), &endptr);
2505 
2506  // if there are many points, the data
2507  // for this var might continue in the
2508  // next line(s)
2509  for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
2510  in >> vertices[j](next_index);
2511  // now we got all values of the first
2512  // variable, so increase the counter
2513  next_index = 1;
2514  }
2515 
2516  // main loop over all variables
2517  for (unsigned int i = 1; i < n_vars; ++i)
2518  {
2519  // if we read all the important
2520  // variables and do not want to
2521  // read further, because we are
2522  // using a structured grid, we can
2523  // stop here (and skip, for
2524  // example, a whole lot of solution
2525  // variables)
2526  if (next_index == dim && structured)
2527  break;
2528 
2529  if ((next_index < dim) && (i == tecplot2deal[next_index]))
2530  {
2531  // we need this line, read it in
2532  for (unsigned int j = 1; j < n_vertices + 1; ++j)
2533  in >> vertices[j](next_index);
2534  ++next_index;
2535  }
2536  else
2537  {
2538  // we do not need this line, read
2539  // it in and discard it
2540  for (unsigned int j = 1; j < n_vertices + 1; ++j)
2541  in >> dummy;
2542  }
2543  }
2544  Assert(next_index == dim, ExcInternalError());
2545  }
2546  else
2547  {
2548  // the data is not blocked, so we get all
2549  // the variables for one point, then the
2550  // next and so on. create a vector to
2551  // hold these components
2552  std::vector<double> vars(n_vars);
2553 
2554  // now fill the first vertex. note, that we
2555  // have already read the first line
2556  // containing the first vertex
2557  std::vector<std::string> first_vertex =
2559  char *endptr;
2560  for (unsigned int d = 0; d < dim; ++d)
2561  vertices[1](d) =
2562  std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
2563 
2564  // read the remaining vertices from the
2565  // list
2566  for (unsigned int v = 2; v < n_vertices + 1; ++v)
2567  {
2568  for (unsigned int i = 0; i < n_vars; ++i)
2569  in >> vars[i];
2570  // fill the vertex
2571  // coordinates. respect the position
2572  // of coordinates in the list of
2573  // variables
2574  for (unsigned int i = 0; i < dim; ++i)
2575  vertices[v](i) = vars[tecplot2deal[i]];
2576  }
2577  }
2578 
2579  if (structured)
2580  {
2581  // this is the part of the code that only
2582  // works in 2d
2583  unsigned int I = IJK[0], J = IJK[1];
2584 
2585  unsigned int cell = 0;
2586  // set up array of cells
2587  for (unsigned int j = 0; j < J - 1; ++j)
2588  for (unsigned int i = 1; i < I; ++i)
2589  {
2590  cells[cell].vertices[0] = i + j * I;
2591  cells[cell].vertices[1] = i + 1 + j * I;
2592  cells[cell].vertices[2] = i + 1 + (j + 1) * I;
2593  cells[cell].vertices[3] = i + (j + 1) * I;
2594  ++cell;
2595  }
2596  Assert(cell == n_cells, ExcInternalError());
2597  std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
2598  unsigned int k = 0;
2599  for (unsigned int i = 1; i < I + 1; ++i)
2600  {
2601  boundary_vertices[k] = i;
2602  ++k;
2603  boundary_vertices[k] = i + (J - 1) * I;
2604  ++k;
2605  }
2606  for (unsigned int j = 1; j < J - 1; ++j)
2607  {
2608  boundary_vertices[k] = 1 + j * I;
2609  ++k;
2610  boundary_vertices[k] = I + j * I;
2611  ++k;
2612  }
2613  Assert(k == boundary_vertices.size(), ExcInternalError());
2614  // delete the duplicated vertices at the
2615  // boundary, which occur, e.g. in c-type
2616  // or o-type grids around a body
2617  // (airfoil). this automatically deletes
2618  // unused vertices as well.
2620  cells,
2621  subcelldata,
2622  boundary_vertices);
2623  }
2624  else
2625  {
2626  // set up array of cells, unstructured
2627  // mode, so the connectivity is
2628  // explicitly given
2629  for (unsigned int i = 0; i < n_cells; ++i)
2630  {
2631  // note that since in the input file
2632  // we found the number of cells at
2633  // the top, there should still be
2634  // input here, so check this:
2635  AssertThrow(in, ExcIO());
2636 
2637  // get the connectivity from the
2638  // input file. the vertices are
2639  // ordered like in the ucd format
2640  for (unsigned int &vertex : cells[i].vertices)
2641  in >> vertex;
2642  }
2643  // do some clean-up on vertices
2644  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2645  }
2646 
2647  // check that no forbidden arrays are
2648  // used. as we do not read in any
2649  // subcelldata, nothing should happen here.
2650  Assert(subcelldata.check_consistency(dim), ExcInternalError());
2651  AssertThrow(in, ExcIO());
2652 
2653  // do some cleanup on cells
2655  cells);
2657  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2658 }
2659 
2660 
2661 
2662 template <int dim, int spacedim>
2663 void
2665 {
2666  Assert(false, ExcNotImplemented());
2667 }
2668 
2669 
2670 
2671 template <int dim, int spacedim>
2672 void
2673 GridIn<dim, spacedim>::read_assimp(const std::string &filename,
2674  const unsigned int mesh_index,
2675  const bool remove_duplicates,
2676  const double tol,
2677  const bool ignore_unsupported_types)
2678 {
2679 #ifdef DEAL_II_WITH_ASSIMP
2680  // Only good for surface grids.
2681  AssertThrow(dim < 3, ExcImpossibleInDim(dim));
2682 
2683  // Create an instance of the Importer class
2684  Assimp::Importer importer;
2685 
2686  // And have it read the given file with some postprocessing
2687  const aiScene *scene =
2688  importer.ReadFile(filename.c_str(),
2689  aiProcess_RemoveComponent |
2690  aiProcess_JoinIdenticalVertices |
2691  aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
2692  aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
2693 
2694  // If the import failed, report it
2695  AssertThrow(scene != nullptr, ExcMessage(importer.GetErrorString()));
2696 
2697  AssertThrow(scene->mNumMeshes != 0,
2698  ExcMessage("Input file contains no meshes."));
2699 
2700  AssertThrow((mesh_index == numbers::invalid_unsigned_int) ||
2701  (mesh_index < scene->mNumMeshes),
2702  ExcMessage("Too few meshes in the file."));
2703 
2704  unsigned int start_mesh =
2705  (mesh_index == numbers::invalid_unsigned_int ? 0 : mesh_index);
2706  unsigned int end_mesh =
2707  (mesh_index == numbers::invalid_unsigned_int ? scene->mNumMeshes :
2708  mesh_index + 1);
2709 
2710  // Deal.II objects are created empty, and then filled with imported file.
2711  std::vector<Point<spacedim>> vertices;
2712  std::vector<CellData<dim>> cells;
2713  SubCellData subcelldata;
2714 
2715  // A series of counters to merge cells.
2716  unsigned int v_offset = 0;
2717  unsigned int c_offset = 0;
2718 
2719  // The index of the mesh will be used as a material index.
2720  for (unsigned int m = start_mesh; m < end_mesh; ++m)
2721  {
2722  const aiMesh *mesh = scene->mMeshes[m];
2723 
2724  // Check that we know what to do with this mesh, otherwise just
2725  // ignore it
2726  if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
2727  {
2728  AssertThrow(ignore_unsupported_types,
2729  ExcMessage("Incompatible mesh " + std::to_string(m) +
2730  "/" + std::to_string(scene->mNumMeshes)));
2731  continue;
2732  }
2733  else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
2734  {
2735  AssertThrow(ignore_unsupported_types,
2736  ExcMessage("Incompatible mesh " + std::to_string(m) +
2737  "/" + std::to_string(scene->mNumMeshes)));
2738  continue;
2739  }
2740  // Vertices
2741  const unsigned int n_vertices = mesh->mNumVertices;
2742  const aiVector3D * mVertices = mesh->mVertices;
2743 
2744  // Faces
2745  const unsigned int n_faces = mesh->mNumFaces;
2746  const aiFace * mFaces = mesh->mFaces;
2747 
2748  vertices.resize(v_offset + n_vertices);
2749  cells.resize(c_offset + n_faces);
2750 
2751  for (unsigned int i = 0; i < n_vertices; ++i)
2752  for (unsigned int d = 0; d < spacedim; ++d)
2753  vertices[i + v_offset][d] = mVertices[i][d];
2754 
2755  unsigned int valid_cell = c_offset;
2756  for (unsigned int i = 0; i < n_faces; ++i)
2757  {
2758  if (mFaces[i].mNumIndices == GeometryInfo<dim>::vertices_per_cell)
2759  {
2760  for (const unsigned int f : GeometryInfo<dim>::vertex_indices())
2761  {
2762  cells[valid_cell].vertices[f] =
2763  mFaces[i].mIndices[f] + v_offset;
2764  }
2765  cells[valid_cell].material_id = m;
2766  ++valid_cell;
2767  }
2768  else
2769  {
2770  AssertThrow(ignore_unsupported_types,
2771  ExcMessage("Face " + std::to_string(i) + " of mesh " +
2772  std::to_string(m) + " has " +
2773  std::to_string(mFaces[i].mNumIndices) +
2774  " vertices. We expected only " +
2777  }
2778  }
2779  cells.resize(valid_cell);
2780 
2781  // The vertices are added all at once. Cells are checked for
2782  // validity, so only valid_cells are now present in the deal.II
2783  // list of cells.
2784  v_offset += n_vertices;
2785  c_offset = valid_cell;
2786  }
2787 
2788  // No cells were read
2789  if (cells.size() == 0)
2790  return;
2791 
2792  if (remove_duplicates)
2793  {
2794  // The function delete_duplicated_vertices() needs to be called more
2795  // than once if a vertex is duplicated more than once. So we keep
2796  // calling it until the number of vertices does not change any more.
2797  unsigned int n_verts = 0;
2798  while (n_verts != vertices.size())
2799  {
2800  n_verts = vertices.size();
2801  std::vector<unsigned int> considered_vertices;
2803  vertices, cells, subcelldata, considered_vertices, tol);
2804  }
2805  }
2806 
2807  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2808  if (dim == spacedim)
2810  cells);
2811 
2813  if (dim == 2)
2814  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2815  else
2816  tria->create_triangulation(vertices, cells, subcelldata);
2817 #else
2818  (void)filename;
2819  (void)mesh_index;
2820  (void)remove_duplicates;
2821  (void)tol;
2822  (void)ignore_unsupported_types;
2823  AssertThrow(false, ExcNeedsAssimp());
2824 #endif
2825 }
2826 
2827 #ifdef DEAL_II_TRILINOS_WITH_SEACAS
2828 // Namespace containing some extra functions for reading ExodusII files
2829 namespace
2830 {
2831  // Convert ExodusII strings to cell types. Use the number of nodes per element
2832  // to disambiguate some cases.
2834  exodusii_name_to_type(const std::string &type_name,
2835  const int n_nodes_per_element)
2836  {
2837  Assert(type_name.size() > 0, ExcInternalError());
2838  // Try to canonify the name by switching to upper case and removing trailing
2839  // numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and PYRAMID13 all
2840  // equal.
2841  std::string type_name_2 = type_name;
2842  std::transform(type_name_2.begin(),
2843  type_name_2.end(),
2844  type_name_2.begin(),
2845  [](unsigned char c) { return std::toupper(c); });
2846  const std::string numbers = "0123456789";
2847  type_name_2.erase(std::find_first_of(type_name_2.begin(),
2848  type_name_2.end(),
2849  numbers.begin(),
2850  numbers.end()),
2851  type_name_2.end());
2852 
2853  if (type_name_2 == "TRI" || type_name_2 == "TRIANGLE")
2854  return ReferenceCell::Type::Tri;
2855  else if (type_name_2 == "QUAD" || type_name_2 == "QUADRILATERAL")
2857  else if (type_name_2 == "SHELL")
2858  {
2859  if (n_nodes_per_element == 3)
2860  return ReferenceCell::Type::Tri;
2861  else
2863  }
2864  else if (type_name_2 == "TET" || type_name_2 == "TETRA" ||
2865  type_name_2 == "TETRAHEDRON")
2866  return ReferenceCell::Type::Tet;
2867  else if (type_name_2 == "PYRA" || type_name_2 == "PYRAMID")
2869  else if (type_name_2 == "WEDGE")
2871  else if (type_name_2 == "HEX" || type_name_2 == "HEXAHEDRON")
2872  return ReferenceCell::Type::Hex;
2873 
2874  Assert(false, ExcNotImplemented());
2876  }
2877 
2878  // Associate deal.II boundary ids with sidesets (a face can be in multiple
2879  // sidesets - to translate we assign each set of side set ids to a
2880  // boundary_id or manifold_id)
2881  template <int dim, int spacedim = dim>
2882  std::pair<SubCellData, std::vector<std::vector<int>>>
2883  read_exodusii_sidesets(const int ex_id,
2884  const int n_side_sets,
2885  const std::vector<CellData<dim>> &cells,
2886  const bool apply_all_indicators_to_manifolds)
2887  {
2888  SubCellData subcelldata;
2889  std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
2890  // boundary id 0 is the default
2891  b_or_m_id_to_sideset_ids.emplace_back();
2892  // deal.II does not support assigning boundary ids with nonzero codimension
2893  // meshes so completely skip this information in that case.
2894  //
2895  // Exodus prints warnings if we try to get empty sets so always check first
2896  if (dim == spacedim && n_side_sets > 0)
2897  {
2898  std::vector<int> side_set_ids(n_side_sets);
2899  int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
2900  AssertThrowExodusII(ierr);
2901 
2902  // First collect all side sets on all boundary faces (indexed here as
2903  // max_faces_per_cell * cell_n + face_n). We then sort and uniquify the
2904  // side sets so that we can convert a set of side set indices into a
2905  // single deal.II boundary or manifold id (and save the correspondence).
2906  constexpr auto max_faces_per_cell = GeometryInfo<dim>::faces_per_cell;
2907  std::map<std::size_t, std::vector<int>> face_side_sets;
2908  for (const int side_set_id : side_set_ids)
2909  {
2910  int n_sides = -1;
2911  int n_distribution_factors = -1;
2912 
2913  ierr = ex_get_set_param(ex_id,
2914  EX_SIDE_SET,
2915  side_set_id,
2916  &n_sides,
2917  &n_distribution_factors);
2918  AssertThrowExodusII(ierr);
2919  if (n_sides > 0)
2920  {
2921  std::vector<int> elements(n_sides);
2922  std::vector<int> faces(n_sides);
2923  ierr = ex_get_set(ex_id,
2924  EX_SIDE_SET,
2925  side_set_id,
2926  elements.data(),
2927  faces.data());
2928  AssertThrowExodusII(ierr);
2929 
2930  // According to the manual (subsection 4.8): "The internal
2931  // number of an element numbering is defined implicitly by the
2932  // order in which it appears in the file. Elements are numbered
2933  // internally (beginning with 1) consecutively across all
2934  // element blocks." Hence element i in Exodus numbering is entry
2935  // i - 1 in the cells array.
2936  for (int side_n = 0; side_n < n_sides; ++side_n)
2937  {
2938  const long element_n = elements[side_n] - 1;
2939  const long face_n = faces[side_n] - 1;
2940  const std::size_t face_id =
2941  element_n * max_faces_per_cell + face_n;
2942  face_side_sets[face_id].push_back(side_set_id);
2943  }
2944  }
2945  }
2946 
2947  // Collect into a sortable data structure:
2948  std::vector<std::pair<std::size_t, std::vector<int>>>
2949  face_id_to_side_sets;
2950  for (auto &pair : face_side_sets)
2951  {
2952  Assert(pair.second.size() > 0, ExcInternalError());
2953  face_id_to_side_sets.push_back(std::move(pair));
2954  }
2955 
2956  // sort by side sets:
2957  std::sort(face_id_to_side_sets.begin(),
2958  face_id_to_side_sets.end(),
2959  [](const auto &a, const auto &b) {
2960  return std::lexicographical_compare(a.second.begin(),
2961  a.second.end(),
2962  b.second.begin(),
2963  b.second.end());
2964  });
2965 
2966  types::boundary_id current_b_or_m_id = 0;
2967  for (const auto &pair : face_id_to_side_sets)
2968  {
2969  const std::size_t face_id = pair.first;
2970  const std::vector<int> &face_sideset_ids = pair.second;
2971  if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
2972  {
2973  // Since we sorted by sideset ids we are guaranteed that if this
2974  // doesn't match the last set then it has not yet been seen
2975  ++current_b_or_m_id;
2976  b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
2977  Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
2978  ExcInternalError());
2979  }
2980  // Record the b_or_m_id of the current face.
2981  const unsigned int local_face_n = face_id % max_faces_per_cell;
2982  const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
2983  const ReferenceCell::Type cell_type =
2984  ReferenceCell::n_vertices_to_type(dim, cell.vertices.size());
2987  const unsigned int deal_face_n =
2988  info.exodusii_face_to_deal_face(local_face_n);
2989  const ReferenceCell::internal::Info::Base &face_info =
2990  ReferenceCell::internal::Info::get_face(cell_type, deal_face_n);
2991 
2992  // The orientation we pick doesn't matter here since when we create
2993  // the Triangulation we will sort the vertices for each CellData
2994  // object created here.
2995  if (dim == 2)
2996  {
2997  CellData<1> boundary_line(face_info.n_vertices());
2998  if (apply_all_indicators_to_manifolds)
2999  boundary_line.manifold_id = current_b_or_m_id;
3000  else
3001  boundary_line.boundary_id = current_b_or_m_id;
3002  for (unsigned int j = 0; j < face_info.n_vertices(); ++j)
3003  boundary_line.vertices[j] =
3004  cell
3005  .vertices[info.face_to_cell_vertices(deal_face_n, j, 0)];
3006 
3007  subcelldata.boundary_lines.push_back(std::move(boundary_line));
3008  }
3009  else if (dim == 3)
3010  {
3011  CellData<2> boundary_quad(face_info.n_vertices());
3012  if (apply_all_indicators_to_manifolds)
3013  boundary_quad.manifold_id = current_b_or_m_id;
3014  else
3015  boundary_quad.boundary_id = current_b_or_m_id;
3016  for (unsigned int j = 0; j < face_info.n_vertices(); ++j)
3017  boundary_quad.vertices[j] =
3018  cell
3019  .vertices[info.face_to_cell_vertices(deal_face_n, j, 0)];
3020 
3021  subcelldata.boundary_quads.push_back(std::move(boundary_quad));
3022  }
3023  }
3024  }
3025 
3026  return std::make_pair(std::move(subcelldata),
3027  std::move(b_or_m_id_to_sideset_ids));
3028  }
3029 } // namespace
3030 #endif
3031 
3032 template <int dim, int spacedim>
3035  const std::string &filename,
3036  const bool apply_all_indicators_to_manifolds)
3037 {
3038 #ifdef DEAL_II_TRILINOS_WITH_SEACAS
3039  // deal.II always uses double precision numbers for geometry
3040  int component_word_size = sizeof(double);
3041  // setting to zero uses the stored word size
3042  int floating_point_word_size = 0;
3043  float ex_version = 0.0;
3044 
3045  const int ex_id = ex_open(filename.c_str(),
3046  EX_READ,
3047  &component_word_size,
3048  &floating_point_word_size,
3049  &ex_version);
3050  AssertThrow(ex_id > 0,
3051  ExcMessage("ExodusII failed to open the specified input file."));
3052 
3053  // Read basic mesh information:
3054  std::vector<char> string_temp(MAX_LINE_LENGTH + 1, '\0');
3055  int mesh_dimension = 0;
3056  int n_nodes = 0;
3057  int n_elements = 0;
3058  int n_element_blocks = 0;
3059  int n_node_sets = 0;
3060  int n_side_sets = 0;
3061 
3062  int ierr = ex_get_init(ex_id,
3063  string_temp.data(),
3064  &mesh_dimension,
3065  &n_nodes,
3066  &n_elements,
3067  &n_element_blocks,
3068  &n_node_sets,
3069  &n_side_sets);
3070  AssertThrowExodusII(ierr);
3071  AssertDimension(mesh_dimension, spacedim);
3072 
3073  // Read nodes:
3074  std::vector<double> xs(n_nodes);
3075  std::vector<double> ys(n_nodes);
3076  std::vector<double> zs(n_nodes);
3077 
3078  ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
3079  AssertThrowExodusII(ierr);
3080 
3081  // Even if there is a node numbering array the values stored inside the
3082  // ExodusII file must use the contiguous, internal ordering (see Section 4.5
3083  // of the manual - "Internal (contiguously numbered) node and element IDs must
3084  // be used for all data structures that contain node or element numbers (IDs),
3085  // including node set node lists, side set element lists, and element
3086  // connectivity.")
3087  std::vector<Point<spacedim>> vertices;
3088  vertices.reserve(n_nodes);
3089  for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
3090  {
3091  switch (spacedim)
3092  {
3093  case 1:
3094  vertices.emplace_back(xs[vertex_n]);
3095  break;
3096  case 2:
3097  vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
3098  break;
3099  case 3:
3100  vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
3101  break;
3102  default:
3103  Assert(spacedim <= 3, ExcNotImplemented());
3104  }
3105  }
3106 
3107  std::vector<int> element_block_ids(n_element_blocks);
3108  ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
3109  AssertThrowExodusII(ierr);
3110 
3111  std::vector<CellData<dim>> cells;
3112  // Elements are grouped together by same reference cell type in element
3113  // blocks. There may be multiple blocks for a single reference cell type,
3114  // but "each element block may contain only one element type".
3115  for (const int element_block_id : element_block_ids)
3116  {
3117  std::fill(string_temp.begin(), string_temp.end(), '\0');
3118  int n_block_elements = 0;
3119  int n_nodes_per_element = 0;
3120  int n_edges_per_element = 0;
3121  int n_faces_per_element = 0;
3122  int n_attributes_per_element = 0;
3123 
3124  // Extract element data.
3125  ierr = ex_get_block(ex_id,
3126  EX_ELEM_BLOCK,
3127  element_block_id,
3128  string_temp.data(),
3129  &n_block_elements,
3130  &n_nodes_per_element,
3131  &n_edges_per_element,
3132  &n_faces_per_element,
3133  &n_attributes_per_element);
3134  AssertThrowExodusII(ierr);
3135  const ReferenceCell::Type type =
3136  exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
3139  // The number of nodes per element may be larger than what we want to
3140  // read - for example, if the Exodus file contains a QUAD9 element, we
3141  // only want to read the first four values and ignore the rest.
3142  Assert(int(info.n_vertices()) <= n_nodes_per_element, ExcInternalError());
3143 
3144  std::vector<int> connection(n_nodes_per_element * n_block_elements);
3145  ierr = ex_get_conn(ex_id,
3146  EX_ELEM_BLOCK,
3147  element_block_id,
3148  connection.data(),
3149  nullptr,
3150  nullptr);
3151  AssertThrowExodusII(ierr);
3152 
3153  for (unsigned int elem_n = 0; elem_n < connection.size();
3154  elem_n += n_nodes_per_element)
3155  {
3156  CellData<dim> cell(info.n_vertices());
3157  for (unsigned int i : info.vertex_indices())
3158  {
3159  cell.vertices[info.exodusii_vertex_to_deal_vertex(i)] =
3160  connection[elem_n + i] - 1;
3161  }
3162  cell.material_id = element_block_id;
3163  cells.push_back(cell);
3164  }
3165  }
3166 
3167  // Extract boundary data.
3168  auto pair = read_exodusii_sidesets<dim, spacedim>(
3169  ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
3170  ierr = ex_close(ex_id);
3171  AssertThrowExodusII(ierr);
3172 
3173  tria->create_triangulation(vertices, cells, pair.first);
3174  ExodusIIData out;
3175  out.id_to_sideset_ids = std::move(pair.second);
3176  return out;
3177 #else
3178  (void)filename;
3179  (void)apply_all_indicators_to_manifolds;
3180  AssertThrow(false, ExcNeedsExodusII());
3181  return {};
3182 #endif
3183 }
3184 
3185 
3186 template <int dim, int spacedim>
3187 void
3189 {
3190  std::string line;
3191  while (in)
3192  {
3193  // get line
3194  getline(in, line);
3195 
3196  // check if this is a line that
3197  // consists only of spaces, and
3198  // if not put the whole thing
3199  // back and return
3200  if (std::find_if(line.begin(), line.end(), [](const char c) {
3201  return c != ' ';
3202  }) != line.end())
3203  {
3204  in.putback('\n');
3205  for (int i = line.length() - 1; i >= 0; --i)
3206  in.putback(line[i]);
3207  return;
3208  }
3209 
3210  // else: go on with next line
3211  }
3212 }
3213 
3214 
3215 
3216 template <int dim, int spacedim>
3217 void
3219  const char comment_start)
3220 {
3221  char c;
3222  // loop over the following comment
3223  // lines
3224  while (in.get(c) && c == comment_start)
3225  // loop over the characters after
3226  // the comment starter
3227  while (in.get() != '\n')
3228  ;
3229 
3230 
3231  // put back first character of
3232  // first non-comment line
3233  if (in)
3234  in.putback(c);
3235 
3236  // at last: skip additional empty lines, if present
3237  skip_empty_lines(in);
3238 }
3239 
3240 
3241 
3242 template <int dim, int spacedim>
3243 void
3245  const std::vector<CellData<dim>> & /*cells*/,
3246  const std::vector<Point<spacedim>> & /*vertices*/,
3247  std::ostream & /*out*/)
3248 {
3249  Assert(false, ExcNotImplemented());
3250 }
3251 
3252 
3253 
3254 template <>
3255 void
3256 GridIn<2>::debug_output_grid(const std::vector<CellData<2>> &cells,
3257  const std::vector<Point<2>> & vertices,
3258  std::ostream & out)
3259 {
3260  double min_x = vertices[cells[0].vertices[0]](0),
3261  max_x = vertices[cells[0].vertices[0]](0),
3262  min_y = vertices[cells[0].vertices[0]](1),
3263  max_y = vertices[cells[0].vertices[0]](1);
3264 
3265  for (unsigned int i = 0; i < cells.size(); ++i)
3266  {
3267  for (const auto vertex : cells[i].vertices)
3268  {
3269  const Point<2> &p = vertices[vertex];
3270 
3271  if (p(0) < min_x)
3272  min_x = p(0);
3273  if (p(0) > max_x)
3274  max_x = p(0);
3275  if (p(1) < min_y)
3276  min_y = p(1);
3277  if (p(1) > max_y)
3278  max_y = p(1);
3279  }
3280 
3281  out << "# cell " << i << std::endl;
3282  Point<2> center;
3283  for (const auto vertex : cells[i].vertices)
3284  center += vertices[vertex];
3285  center /= 4;
3286 
3287  out << "set label \"" << i << "\" at " << center(0) << ',' << center(1)
3288  << " center" << std::endl;
3289 
3290  // first two line right direction
3291  for (unsigned int f = 0; f < 2; ++f)
3292  out << "set arrow from " << vertices[cells[i].vertices[f]](0) << ','
3293  << vertices[cells[i].vertices[f]](1) << " to "
3294  << vertices[cells[i].vertices[(f + 1) % 4]](0) << ','
3295  << vertices[cells[i].vertices[(f + 1) % 4]](1) << std::endl;
3296  // other two lines reverse direction
3297  for (unsigned int f = 2; f < 4; ++f)
3298  out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]](0)
3299  << ',' << vertices[cells[i].vertices[(f + 1) % 4]](1) << " to "
3300  << vertices[cells[i].vertices[f]](0) << ','
3301  << vertices[cells[i].vertices[f]](1) << std::endl;
3302  out << std::endl;
3303  }
3304 
3305 
3306  out << std::endl
3307  << "set nokey" << std::endl
3308  << "pl [" << min_x << ':' << max_x << "][" << min_y << ':' << max_y
3309  << "] " << min_y << std::endl
3310  << "pause -1" << std::endl;
3311 }
3312 
3313 
3314 
3315 template <>
3316 void
3317 GridIn<3>::debug_output_grid(const std::vector<CellData<3>> &cells,
3318  const std::vector<Point<3>> & vertices,
3319  std::ostream & out)
3320 {
3321  for (const auto &cell : cells)
3322  {
3323  // line 0
3324  out << vertices[cell.vertices[0]] << std::endl
3325  << vertices[cell.vertices[1]] << std::endl
3326  << std::endl
3327  << std::endl;
3328  // line 1
3329  out << vertices[cell.vertices[1]] << std::endl
3330  << vertices[cell.vertices[2]] << std::endl
3331  << std::endl
3332  << std::endl;
3333  // line 2
3334  out << vertices[cell.vertices[3]] << std::endl
3335  << vertices[cell.vertices[2]] << std::endl
3336  << std::endl
3337  << std::endl;
3338  // line 3
3339  out << vertices[cell.vertices[0]] << std::endl
3340  << vertices[cell.vertices[3]] << std::endl
3341  << std::endl
3342  << std::endl;
3343  // line 4
3344  out << vertices[cell.vertices[4]] << std::endl
3345  << vertices[cell.vertices[5]] << std::endl
3346  << std::endl
3347  << std::endl;
3348  // line 5
3349  out << vertices[cell.vertices[5]] << std::endl
3350  << vertices[cell.vertices[6]] << std::endl
3351  << std::endl
3352  << std::endl;
3353  // line 6
3354  out << vertices[cell.vertices[7]] << std::endl
3355  << vertices[cell.vertices[6]] << std::endl
3356  << std::endl
3357  << std::endl;
3358  // line 7
3359  out << vertices[cell.vertices[4]] << std::endl
3360  << vertices[cell.vertices[7]] << std::endl
3361  << std::endl
3362  << std::endl;
3363  // line 8
3364  out << vertices[cell.vertices[0]] << std::endl
3365  << vertices[cell.vertices[4]] << std::endl
3366  << std::endl
3367  << std::endl;
3368  // line 9
3369  out << vertices[cell.vertices[1]] << std::endl
3370  << vertices[cell.vertices[5]] << std::endl
3371  << std::endl
3372  << std::endl;
3373  // line 10
3374  out << vertices[cell.vertices[2]] << std::endl
3375  << vertices[cell.vertices[6]] << std::endl
3376  << std::endl
3377  << std::endl;
3378  // line 11
3379  out << vertices[cell.vertices[3]] << std::endl
3380  << vertices[cell.vertices[7]] << std::endl
3381  << std::endl
3382  << std::endl;
3383  }
3384 }
3385 
3386 
3387 
3388 template <int dim, int spacedim>
3389 void
3390 GridIn<dim, spacedim>::read(const std::string &filename, Format format)
3391 {
3392  // Search file class for meshes
3393  PathSearch search("MESH");
3394  std::string name;
3395  // Open the file and remember its name
3396  if (format == Default)
3397  name = search.find(filename);
3398  else
3399  name = search.find(filename, default_suffix(format));
3400 
3401 
3402  if (format == Default)
3403  {
3404  const std::string::size_type slashpos = name.find_last_of('/');
3405  const std::string::size_type dotpos = name.find_last_of('.');
3406  if (dotpos < name.length() &&
3407  (dotpos > slashpos || slashpos == std::string::npos))
3408  {
3409  std::string ext = name.substr(dotpos + 1);
3410  format = parse_format(ext);
3411  }
3412  }
3413 
3414  if (format == assimp)
3415  {
3416  read_assimp(name);
3417  }
3418  else if (format == exodusii)
3419  {
3420  read_exodusii(name);
3421  }
3422  else
3423  {
3424  std::ifstream in(name.c_str());
3425  read(in, format);
3426  }
3427 }
3428 
3429 
3430 template <int dim, int spacedim>
3431 void
3432 GridIn<dim, spacedim>::read(std::istream &in, Format format)
3433 {
3434  if (format == Default)
3435  format = default_format;
3436 
3437  switch (format)
3438  {
3439  case dbmesh:
3440  read_dbmesh(in);
3441  return;
3442 
3443  case msh:
3444  read_msh(in);
3445  return;
3446 
3447  case vtk:
3448  read_vtk(in);
3449  return;
3450 
3451  case vtu:
3452  read_vtu(in);
3453  return;
3454 
3455  case unv:
3456  read_unv(in);
3457  return;
3458 
3459  case ucd:
3460  read_ucd(in);
3461  return;
3462 
3463  case abaqus:
3464  read_abaqus(in);
3465  return;
3466 
3467  case xda:
3468  read_xda(in);
3469  return;
3470 
3471  case tecplot:
3472  read_tecplot(in);
3473  return;
3474 
3475  case assimp:
3476  Assert(false,
3477  ExcMessage("There is no read_assimp(istream &) function. "
3478  "Use the read_assimp(string &filename, ...) "
3479  "functions, instead."));
3480  return;
3481 
3482  case exodusii:
3483  Assert(false,
3484  ExcMessage("There is no read_exodusii(istream &) function. "
3485  "Use the read_exodusii(string &filename, ...) "
3486  "function, instead."));
3487  return;
3488 
3489  case Default:
3490  break;
3491  }
3492  Assert(false, ExcInternalError());
3493 }
3494 
3495 
3496 
3497 template <int dim, int spacedim>
3498 std::string
3500 {
3501  switch (format)
3502  {
3503  case dbmesh:
3504  return ".dbmesh";
3505  case exodusii:
3506  return ".e";
3507  case msh:
3508  return ".msh";
3509  case vtk:
3510  return ".vtk";
3511  case vtu:
3512  return ".vtu";
3513  case unv:
3514  return ".unv";
3515  case ucd:
3516  return ".inp";
3517  case abaqus:
3518  return ".inp"; // Typical suffix for Abaqus mesh files conflicts with
3519  // UCD.
3520  case xda:
3521  return ".xda";
3522  case tecplot:
3523  return ".dat";
3524  default:
3525  Assert(false, ExcNotImplemented());
3526  return ".unknown_format";
3527  }
3528 }
3529 
3530 
3531 
3532 template <int dim, int spacedim>
3534 GridIn<dim, spacedim>::parse_format(const std::string &format_name)
3535 {
3536  if (format_name == "dbmesh")
3537  return dbmesh;
3538 
3539  if (format_name == "exodusii")
3540  return exodusii;
3541 
3542  if (format_name == "msh")
3543  return msh;
3544 
3545  if (format_name == "unv")
3546  return unv;
3547 
3548  if (format_name == "vtk")
3549  return vtk;
3550 
3551  if (format_name == "vtu")
3552  return vtu;
3553 
3554  // This is also the typical extension of Abaqus input files.
3555  if (format_name == "inp")
3556  return ucd;
3557 
3558  if (format_name == "ucd")
3559  return ucd;
3560 
3561  if (format_name == "xda")
3562  return xda;
3563 
3564  if (format_name == "tecplot")
3565  return tecplot;
3566 
3567  if (format_name == "dat")
3568  return tecplot;
3569 
3570  if (format_name == "plt")
3571  // Actually, this is the extension for the
3572  // tecplot binary format, which we do not
3573  // support right now. However, some people
3574  // tend to create tecplot ascii files with
3575  // the extension 'plt' instead of
3576  // 'dat'. Thus, include this extension
3577  // here. If it actually is a binary file,
3578  // the read_tecplot() function will fail
3579  // and throw an exception, anyway.
3580  return tecplot;
3581 
3582  AssertThrow(false, ExcInvalidState());
3583  // return something weird
3584  return Format(Default);
3585 }
3586 
3587 
3588 
3589 template <int dim, int spacedim>
3590 std::string
3592 {
3593  return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
3594 }
3595 
3596 
3597 
3598 namespace
3599 {
3600  template <int dim, int spacedim>
3601  Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
3602  : tolerance(5e-16) // Used to offset Cubit tolerance error when outputting
3603  // value close to zero
3604  {
3605  AssertThrow(spacedim == 2 || spacedim == 3, ExcNotImplemented());
3606  }
3607 
3608 
3609 
3610  // Convert from a string to some other data type
3611  // Reference: http://www.codeguru.com/forum/showthread.php?t=231054
3612  template <class T>
3613  bool
3614  from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &))
3615  {
3616  std::istringstream iss(s);
3617  return !(iss >> f >> t).fail();
3618  }
3619 
3620 
3621 
3622  // Extract an integer from a string
3623  int
3624  extract_int(const std::string &s)
3625  {
3626  std::string tmp;
3627  for (const char c : s)
3628  {
3629  if (isdigit(c))
3630  {
3631  tmp += c;
3632  }
3633  }
3634 
3635  int number = 0;
3636  from_string(number, tmp, std::dec);
3637  return number;
3638  }
3639 
3640 
3641 
3642  template <int dim, int spacedim>
3643  void
3644  Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
3645  {
3646  // References:
3647  // http://www.egr.msu.edu/software/abaqus/Documentation/docs/v6.7/books/usb/default.htm?startat=pt01ch02.html
3648  // http://www.cprogramming.com/tutorial/string.html
3649 
3650  AssertThrow(input_stream, ExcIO());
3651  std::string line;
3652 
3653  while (std::getline(input_stream, line))
3654  {
3655  cont:
3656  std::transform(line.begin(), line.end(), line.begin(), ::toupper);
3657 
3658  if (line.compare("*HEADING") == 0 || line.compare(0, 2, "**") == 0 ||
3659  line.compare(0, 5, "*PART") == 0)
3660  {
3661  // Skip header and comments
3662  while (std::getline(input_stream, line))
3663  {
3664  if (line[0] == '*')
3665  goto cont; // My eyes, they burn!
3666  }
3667  }
3668  else if (line.compare(0, 5, "*NODE") == 0)
3669  {
3670  // Extract list of vertices
3671  // Header line might be:
3672  // *NODE, NSET=ALLNODES
3673  // *NODE
3674 
3675  // Contains lines in the form:
3676  // Index, x, y, z
3677  while (std::getline(input_stream, line))
3678  {
3679  if (line[0] == '*')
3680  goto cont;
3681 
3682  std::vector<double> node(spacedim + 1);
3683 
3684  std::istringstream iss(line);
3685  char comma;
3686  for (unsigned int i = 0; i < spacedim + 1; ++i)
3687  iss >> node[i] >> comma;
3688 
3689  node_list.push_back(node);
3690  }
3691  }
3692  else if (line.compare(0, 8, "*ELEMENT") == 0)
3693  {
3694  // Element construction.
3695  // There are different header formats, the details
3696  // of which we're not particularly interested in except
3697  // whether they represent quads or hexahedrals.
3698  // *ELEMENT, TYPE=S4R, ELSET=EB<material id>
3699  // *ELEMENT, TYPE=C3D8R, ELSET=EB<material id>
3700  // *ELEMENT, TYPE=C3D8
3701  // Elements itself (n=4 or n=8):
3702  // Index, i[0], ..., i[n]
3703 
3704  int material = 0;
3705  // Scan for material id
3706  {
3707  const std::string before_material = "ELSET=EB";
3708  const std::size_t idx = line.find(before_material);
3709  if (idx != std::string::npos)
3710  {
3711  from_string(material,
3712  line.substr(idx + before_material.size()),
3713  std::dec);
3714  }
3715  }
3716 
3717  // Read ELEMENT definition
3718  while (std::getline(input_stream, line))
3719  {
3720  if (line[0] == '*')
3721  goto cont;
3722 
3723  std::istringstream iss(line);
3724  char comma;
3725 
3726  // We will store the material id in the zeroth entry of the
3727  // vector and the rest of the elements represent the global
3728  // node numbers
3729  const unsigned int n_data_per_cell =
3731  std::vector<double> cell(n_data_per_cell);
3732  for (unsigned int i = 0; i < n_data_per_cell; ++i)
3733  iss >> cell[i] >> comma;
3734 
3735  // Overwrite cell index from file by material
3736  cell[0] = static_cast<double>(material);
3737  cell_list.push_back(cell);
3738  }
3739  }
3740  else if (line.compare(0, 8, "*SURFACE") == 0)
3741  {
3742  // Extract the definitions of boundary surfaces
3743  // Old format from Cubit:
3744  // *SURFACE, NAME=SS<boundary indicator>
3745  // <element index>, S<face number>
3746  // Abaqus default format:
3747  // *SURFACE, TYPE=ELEMENT, NAME=SURF-<indicator>
3748 
3749  // Get name of the surface and extract id from it;
3750  // this will be the boundary indicator
3751  const std::string name_key = "NAME=";
3752  const std::size_t name_idx_start =
3753  line.find(name_key) + name_key.size();
3754  std::size_t name_idx_end = line.find(',', name_idx_start);
3755  if (name_idx_end == std::string::npos)
3756  {
3757  name_idx_end = line.size();
3758  }
3759  const int b_indicator = extract_int(
3760  line.substr(name_idx_start, name_idx_end - name_idx_start));
3761 
3762  // Read SURFACE definition
3763  // Note that the orientation of the faces is embedded within the
3764  // definition of each "set" of faces that comprise the surface
3765  // These are either marked by an "S" or "E" in 3d or 2d
3766  // respectively.
3767  while (std::getline(input_stream, line))
3768  {
3769  if (line[0] == '*')
3770  goto cont;
3771 
3772  // Change all characters to upper case
3773  std::transform(line.begin(),
3774  line.end(),
3775  line.begin(),
3776  ::toupper);
3777 
3778  // Surface can be created from ELSET, or directly from cells
3779  // If elsets_list contains a key with specific name - refers to
3780  // that ELSET, otherwise refers to cell
3781  std::istringstream iss(line);
3782  int el_idx;
3783  int face_number;
3784  char temp;
3785 
3786  // Get relevant faces, taking into account the element
3787  // orientation
3788  std::vector<double> quad_node_list;
3789  const std::string elset_name = line.substr(0, line.find(','));
3790  if (elsets_list.count(elset_name) != 0)
3791  {
3792  // Surface refers to ELSET
3793  std::string stmp;
3794  iss >> stmp >> temp >> face_number;
3795 
3796  const std::vector<int> cells = elsets_list[elset_name];
3797  for (const int cell : cells)
3798  {
3799  el_idx = cell;
3800  quad_node_list =
3801  get_global_node_numbers(el_idx, face_number);
3802  quad_node_list.insert(quad_node_list.begin(),
3803  b_indicator);
3804 
3805  face_list.push_back(quad_node_list);
3806  }
3807  }
3808  else
3809  {
3810  // Surface refers directly to elements
3811  char comma;
3812  iss >> el_idx >> comma >> temp >> face_number;
3813  quad_node_list =
3814  get_global_node_numbers(el_idx, face_number);
3815  quad_node_list.insert(quad_node_list.begin(), b_indicator);
3816 
3817  face_list.push_back(quad_node_list);
3818  }
3819  }
3820  }
3821  else if (line.compare(0, 6, "*ELSET") == 0)
3822  {
3823  // Get ELSET name.
3824  // Materials are attached to elsets with specific name
3825  std::string elset_name;
3826  {
3827  const std::string elset_key = "*ELSET, ELSET=";
3828  const std::size_t idx = line.find(elset_key);
3829  if (idx != std::string::npos)
3830  {
3831  const std::string comma = ",";
3832  const std::size_t first_comma = line.find(comma);
3833  const std::size_t second_comma =
3834  line.find(comma, first_comma + 1);
3835  const std::size_t elset_name_start =
3836  line.find(elset_key) + elset_key.size();
3837  elset_name = line.substr(elset_name_start,
3838  second_comma - elset_name_start);
3839  }
3840  }
3841 
3842  // There are two possibilities of storing cells numbers in ELSET:
3843  // 1. If the header contains the 'GENERATE' keyword, then the next
3844  // line describes range of cells as:
3845  // cell_id_start, cell_id_end, cell_step
3846  // 2. If the header does not contain the 'GENERATE' keyword, then
3847  // the next lines contain cells numbers
3848  std::vector<int> elements;
3849  const std::size_t generate_idx = line.find("GENERATE");
3850  if (generate_idx != std::string::npos)
3851  {
3852  // Option (1)
3853  std::getline(input_stream, line);
3854  std::istringstream iss(line);
3855  char comma;
3856  int elid_start;
3857  int elid_end;
3858  int elis_step = 1; // Default if case stride not provided
3859 
3860  // Some files don't have the stride size
3861  // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp to
3862  // ./grids/abaqus/2d/2d_test_abaqus.inp
3863  iss >> elid_start >> comma >> elid_end;
3864  AssertThrow(comma == ',',
3865  ExcMessage(
3866  std::string(
3867  "While reading an ABAQUS file, the reader "
3868  "expected a comma but found a <") +
3869  comma + "> in the line <" + line + ">."));
3870  AssertThrow(
3871  elid_start <= elid_end,
3872  ExcMessage(
3873  std::string(
3874  "While reading an ABAQUS file, the reader encountered "
3875  "a GENERATE statement in which the upper bound <") +
3876  Utilities::int_to_string(elid_end) +
3877  "> for the element numbers is not larger or equal "
3878  "than the lower bound <" +
3879  Utilities::int_to_string(elid_start) + ">."));
3880 
3881  // https://stackoverflow.com/questions/8046357/how-do-i-check-if-a-stringstream-variable-is-empty-null
3882  if (iss.rdbuf()->in_avail() != 0)
3883  iss >> comma >> elis_step;
3884  AssertThrow(comma == ',',
3885  ExcMessage(
3886  std::string(
3887  "While reading an ABAQUS file, the reader "
3888  "expected a comma but found a <") +
3889  comma + "> in the line <" + line + ">."));
3890 
3891  for (int i = elid_start; i <= elid_end; i += elis_step)
3892  elements.push_back(i);
3893  elsets_list[elset_name] = elements;
3894 
3895  std::getline(input_stream, line);
3896  }
3897  else
3898  {
3899  // Option (2)
3900  while (std::getline(input_stream, line))
3901  {
3902  if (line[0] == '*')
3903  break;
3904 
3905  std::istringstream iss(line);
3906  char comma;
3907  int elid;
3908  while (!iss.eof())
3909  {
3910  iss >> elid >> comma;
3911  AssertThrow(
3912  comma == ',',
3913  ExcMessage(
3914  std::string(
3915  "While reading an ABAQUS file, the reader "
3916  "expected a comma but found a <") +
3917  comma + "> in the line <" + line + ">."));
3918 
3919  elements.push_back(elid);
3920  }
3921  }
3922 
3923  elsets_list[elset_name] = elements;
3924  }
3925 
3926  goto cont;
3927  }
3928  else if (line.compare(0, 5, "*NSET") == 0)
3929  {
3930  // Skip nodesets; we have no use for them
3931  while (std::getline(input_stream, line))
3932  {
3933  if (line[0] == '*')
3934  goto cont;
3935  }
3936  }
3937  else if (line.compare(0, 14, "*SOLID SECTION") == 0)
3938  {
3939  // The ELSET name, which describes a section for particular material
3940  const std::string elset_key = "ELSET=";
3941  const std::size_t elset_start =
3942  line.find("ELSET=") + elset_key.size();
3943  const std::size_t elset_end = line.find(',', elset_start + 1);
3944  const std::string elset_name =
3945  line.substr(elset_start, elset_end - elset_start);
3946 
3947  // Solid material definition.
3948  // We assume that material id is taken from material name,
3949  // eg. "Material-1" -> ID=1
3950  const std::string material_key = "MATERIAL=";
3951  const std::size_t last_equal =
3952  line.find("MATERIAL=") + material_key.size();
3953  const std::size_t material_id_start = line.find('-', last_equal);
3954  int material_id = 0;
3955  from_string(material_id,
3956  line.substr(material_id_start + 1),
3957  std::dec);
3958 
3959  // Assign material id to cells
3960  const std::vector<int> &elset_cells = elsets_list[elset_name];
3961  for (const int elset_cell : elset_cells)
3962  {
3963  const int cell_id = elset_cell - 1;
3964  cell_list[cell_id][0] = material_id;
3965  }
3966  }
3967  // Note: All other lines / entries are ignored
3968  }
3969  }
3970 
3971  template <int dim, int spacedim>
3972  std::vector<double>
3973  Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
3974  const int face_cell_no,
3975  const int face_cell_face_no) const
3976  {
3977  std::vector<double> quad_node_list(GeometryInfo<dim>::vertices_per_face);
3978 
3979  // These orderings were reverse engineered by hand and may
3980  // conceivably be erroneous.
3981  // TODO: Currently one test (2d unstructured mesh) in the test
3982  // suite fails, presumably because of an ordering issue.
3983  if (dim == 2)
3984  {
3985  if (face_cell_face_no == 1)
3986  {
3987  quad_node_list[0] = cell_list[face_cell_no - 1][1];
3988  quad_node_list[1] = cell_list[face_cell_no - 1][2];
3989  }
3990  else if (face_cell_face_no == 2)
3991  {
3992  quad_node_list[0] = cell_list[face_cell_no - 1][2];
3993  quad_node_list[1] = cell_list[face_cell_no - 1][3];
3994  }
3995  else if (face_cell_face_no == 3)
3996  {
3997  quad_node_list[0] = cell_list[face_cell_no - 1][3];
3998  quad_node_list[1] = cell_list[face_cell_no - 1][4];
3999  }
4000  else if (face_cell_face_no == 4)
4001  {
4002  quad_node_list[0] = cell_list[face_cell_no - 1][4];
4003  quad_node_list[1] = cell_list[face_cell_no - 1][1];
4004  }
4005  else
4006  {
4007  AssertThrow(face_cell_face_no <= 4,
4008  ExcMessage("Invalid face number in 2d"));
4009  }
4010  }
4011  else if (dim == 3)
4012  {
4013  if (face_cell_face_no == 1)
4014  {
4015  quad_node_list[0] = cell_list[face_cell_no - 1][1];
4016  quad_node_list[1] = cell_list[face_cell_no - 1][4];
4017  quad_node_list[2] = cell_list[face_cell_no - 1][3];
4018  quad_node_list[3] = cell_list[face_cell_no - 1][2];
4019  }
4020  else if (face_cell_face_no == 2)
4021  {
4022  quad_node_list[0] = cell_list[face_cell_no - 1][5];
4023  quad_node_list[1] = cell_list[face_cell_no - 1][8];
4024  quad_node_list[2] = cell_list[face_cell_no - 1][7];
4025  quad_node_list[3] = cell_list[face_cell_no - 1][6];
4026  }
4027  else if (face_cell_face_no == 3)
4028  {
4029  quad_node_list[0] = cell_list[face_cell_no - 1][1];
4030  quad_node_list[1] = cell_list[face_cell_no - 1][2];
4031  quad_node_list[2] = cell_list[face_cell_no - 1][6];
4032  quad_node_list[3] = cell_list[face_cell_no - 1][5];
4033  }
4034  else if (face_cell_face_no == 4)
4035  {
4036  quad_node_list[0] = cell_list[face_cell_no - 1][2];
4037  quad_node_list[1] = cell_list[face_cell_no - 1][3];
4038  quad_node_list[2] = cell_list[face_cell_no - 1][7];
4039  quad_node_list[3] = cell_list[face_cell_no - 1][6];
4040  }
4041  else if (face_cell_face_no == 5)
4042  {
4043  quad_node_list[0] = cell_list[face_cell_no - 1][3];
4044  quad_node_list[1] = cell_list[face_cell_no - 1][4];
4045  quad_node_list[2] = cell_list[face_cell_no - 1][8];
4046  quad_node_list[3] = cell_list[face_cell_no - 1][7];
4047  }
4048  else if (face_cell_face_no == 6)
4049  {
4050  quad_node_list[0] = cell_list[face_cell_no - 1][1];
4051  quad_node_list[1] = cell_list[face_cell_no - 1][5];
4052  quad_node_list[2] = cell_list[face_cell_no - 1][8];
4053  quad_node_list[3] = cell_list[face_cell_no - 1][4];
4054  }
4055  else
4056  {
4057  AssertThrow(face_cell_no <= 6,
4058  ExcMessage("Invalid face number in 3d"));
4059  }
4060  }
4061  else
4062  {
4063  AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
4064  }
4065 
4066  return quad_node_list;
4067  }
4068 
4069  template <int dim, int spacedim>
4070  void
4071  Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output) const
4072  {
4073  // References:
4074  // http://www.dealii.org/developer/doxygen/deal.II/structGeometryInfo.html
4075  // http://people.scs.fsu.edu/~burkardt/data/ucd/ucd.html
4076 
4077  AssertThrow(output, ExcIO());
4078 
4079  // save old formatting options
4080  const boost::io::ios_base_all_saver formatting_saver(output);
4081 
4082  // Write out title - Note: No other commented text can be inserted below the
4083  // title in a UCD file
4084  output << "# Abaqus to UCD mesh conversion" << std::endl;
4085  output << "# Mesh type: AVS UCD" << std::endl;
4086 
4087  // ========================================================
4088  // ASCII UCD File Format
4089  // The input file cannot contain blank lines or lines with leading blanks.
4090  // Comments, if present, must precede all data in the file.
4091  // Comments within the data will cause read errors.
4092  // The general order of the data is as follows:
4093  // 1. Numbers defining the overall structure, including the number of nodes,
4094  // the number of cells, and the length of the vector of data associated
4095  // with the nodes, cells, and the model.
4096  // e.g. 1:
4097  // <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
4098  // e.g. 2:
4099  // n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
4100  // outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
4101  // 2. For each node, its node id and the coordinates of that node in space.
4102  // Node-ids must be integers, but any number including non sequential
4103  // numbers can be used. Mid-edge nodes are treated like any other node.
4104  // 3. For each cell: its cell-id, material, cell type (hexahedral, pyramid,
4105  // etc.), and the list of node-ids that correspond to each of the cell's
4106  // vertices. The below table specifies the different cell types and the
4107  // keyword used to represent them in the file.
4108 
4109  // Write out header
4110  output << node_list.size() << "\t" << (cell_list.size() + face_list.size())
4111  << "\t0\t0\t0" << std::endl;
4112 
4113  output.width(16);
4114  output.precision(8);
4115 
4116  // Write out node numbers
4117  // Loop over all nodes
4118  for (const auto &node : node_list)
4119  {
4120  // Node number
4121  output << node[0] << "\t";
4122 
4123  // Node coordinates
4124  output.setf(std::ios::scientific, std::ios::floatfield);
4125  for (unsigned int jj = 1; jj < spacedim + 1; ++jj)
4126  {
4127  // invoke tolerance -> set points close to zero equal to zero
4128  if (std::abs(node[jj]) > tolerance)
4129  output << static_cast<double>(node[jj]) << "\t";
4130  else
4131  output << 0.0 << "\t";
4132  }
4133  if (spacedim == 2)
4134  output << 0.0 << "\t";
4135 
4136  output << std::endl;
4137  output.unsetf(std::ios::floatfield);
4138  }
4139 
4140  // Write out cell node numbers
4141  for (unsigned int ii = 0; ii < cell_list.size(); ++ii)
4142  {
4143  output << ii + 1 << "\t" << cell_list[ii][0] << "\t"
4144  << (dim == 2 ? "quad" : "hex") << "\t";
4145  for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
4146  ++jj)
4147  output << cell_list[ii][jj] << "\t";
4148 
4149  output << std::endl;
4150  }
4151 
4152  // Write out quad node numbers
4153  for (unsigned int ii = 0; ii < face_list.size(); ++ii)
4154  {
4155  output << ii + 1 << "\t" << face_list[ii][0] << "\t"
4156  << (dim == 2 ? "line" : "quad") << "\t";
4157  for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
4158  ++jj)
4159  output << face_list[ii][jj] << "\t";
4160 
4161  output << std::endl;
4162  }
4163  }
4164 } // namespace
4165 
4166 
4167 // explicit instantiations
4168 #include "grid_in.inst"
4169 
std::vector< CellData< 1 > > boundary_lines
static std::string get_format_names()
Definition: grid_in.cc:3591
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:213
std::vector< std::vector< int > > id_to_sideset_ids
Definition: grid_in.h:586
static ::ExceptionBase & ExcNeedsExodusII()
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
std::vector< unsigned char > decode_base64(const std::string &base64_input)
Definition: utilities.cc:450
Use read_abaqus()
Definition: grid_in.h:322
types::global_dof_index size_type
Definition: cuda_kernels.h:45
static ::ExceptionBase & ExcIO()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
Use read_vtu()
Definition: grid_in.h:334
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11546
void read_msh(std::istream &in)
Definition: grid_in.cc:1494
const ReferenceCell::internal::Info::Base & get_cell(const ReferenceCell::Type &type)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
Format
Definition: grid_in.h:313
std::vector< unsigned int > vertices
std::pair< int, unsigned int > get_integer_at_position(const std::string &name, const unsigned int position)
Definition: utilities.cc:854
Use read_vtk()
Definition: grid_in.h:332
static ::ExceptionBase & ExcUnknownElementType(int arg1)
std::vector< std::string > break_text_into_lines(const std::string &original_text, const unsigned int width, const char delimiter=' ')
Definition: utilities.cc:759
Use read_unv()
Definition: grid_in.h:318
Use read_xda()
Definition: grid_in.h:326
unsigned int material_id
Definition: types.h:152
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
Point< 2 > second
Definition: grid_out.cc:4341
std::string decompress(const std::string &compressed_input)
Definition: utilities.cc:415
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static Format parse_format(const std::string &format_name)
Definition: grid_in.cc:3534
static void skip_empty_lines(std::istream &in)
Definition: grid_in.cc:3188
static ::ExceptionBase & ExcDBMESHWrongDimension(int arg1)
void read_vtk(std::istream &in)
Definition: grid_in.cc:123
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:623
void read_vtu(std::istream &in)
Definition: grid_in.cc:572
const ReferenceCell::internal::Info::Base & get_face(const ReferenceCell::Type &type, const unsigned int face_no)
Use read_ucd()
Definition: grid_in.h:320
void read_dbmesh(std::istream &in)
Definition: grid_in.cc:1169
void read_tecplot(std::istream &in)
Definition: grid_in.cc:2664
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
bool check_consistency(const unsigned int dim) const
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static void parse_tecplot_header(std::string &header, std::vector< unsigned int > &tecplot2deal, unsigned int &n_vars, unsigned int &n_vertices, unsigned int &n_cells, std::vector< unsigned int > &IJK, bool &structured, bool &blocked)
Definition: grid_in.cc:2189
static const char T
static ::ExceptionBase & ExcInvalidVertexIndexGmsh(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1466
SmartPointer< Triangulation< dim, spacedim >, GridIn< dim, spacedim > > tria
Definition: grid_in.h:767
Use read_msh()
Definition: grid_in.h:328
static void skip_comment_lines(std::istream &in, const char comment_start)
Definition: grid_in.cc:3218
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
bool match_at_string_start(const std::string &name, const std::string &pattern)
Definition: utilities.cc:839
void read_ucd(std::istream &in, const bool apply_all_indicators_to_manifolds=false)
Definition: grid_in.cc:848
static ::ExceptionBase & ExcNeedsAssimp()
std::string to_string(const T &t)
Definition: patterns.h:2342
static ::ExceptionBase & ExcInvalidDBMESHInput(std::string arg1)
Point< 3 > vertices[4]
GridIn()
Definition: grid_in.cc:97
static void debug_output_grid(const std::vector< CellData< dim >> &cells, const std::vector< Point< spacedim >> &vertices, std::ostream &out)
Definition: grid_in.cc:3244
Use read_tecplot()
Definition: grid_in.h:330
ExodusIIData read_exodusii(const std::string &filename, const bool apply_all_indicators_to_manifolds=false)
Definition: grid_in.cc:3034
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcUnknownIdentifier(std::string arg1)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:11979
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::string find(const std::string &filename, const char *open_mode="r")
Definition: path_search.cc:171
types::manifold_id manifold_id
void read_xda(std::istream &in)
Definition: grid_in.cc:1336
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
T min(const T &t, const MPI_Comm &mpi_communicator)
std::vector< CellData< 2 > > boundary_quads
Use read_dbmesh()
Definition: grid_in.h:324
virtual unsigned int n_vertices() const
Type n_vertices_to_type(const int dim, const unsigned int n_vertices)
static std::string default_suffix(const Format format)
Definition: grid_in.cc:3499
void read_unv(std::istream &in)
Definition: grid_in.cc:599
void read(std::istream &in, Format format=Default)
Definition: grid_in.cc:3432
Point< 3 > center
static ::ExceptionBase & ExcNoTriangulationSelected()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidDBMeshFormat()
static void invert_all_cells_of_negative_grid(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &original_cells)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
void read_assimp(const std::string &filename, const unsigned int mesh_index=numbers::invalid_unsigned_int, const bool remove_duplicates=true, const double tol=1e-12, const bool ignore_unsupported_element_types=true)
Definition: grid_in.cc:2673
Use GridIn::default_format stored in this object.
Definition: grid_in.h:316
#define DEAL_II_FALLTHROUGH
Definition: config.h:153
void read_abaqus(std::istream &in, const bool apply_all_indicators_to_manifolds=false)
Definition: grid_in.cc:1113
static ::ExceptionBase & ExcInvalidGMSHInput(std::string arg1)
static ::ExceptionBase & ExcGmshUnsupportedGeometry(int arg1)
Use read_exodusii()
Definition: grid_in.h:338
const types::material_id invalid_material_id
Definition: types.h:228
static ::ExceptionBase & ExcUnknownSectionType(int arg1)
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition: grid_in.cc:114
static ::ExceptionBase & ExcInvalidVertexIndex(int arg1, int arg2)
T max(const T &t, const MPI_Comm &mpi_communicator)
Use read_assimp()
Definition: grid_in.h:336
Format default_format
Definition: grid_in.h:832
static ::ExceptionBase & ExcGmshNoCellInformation()
static ::ExceptionBase & ExcInternalError()
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition: grid_tools.cc:728