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