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