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