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