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