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