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