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