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