Reference documentation for deal.II version GIT relicensing-468-ge3ce94fd06 2024-04-23 15:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 {
2792 // Patterns::Tools::to_value throws an exception, if it can not
2793 // convert name to a map from string to int.
2794 try
2795 {
2796 std::map<std::string, int> id_names;
2797 Patterns::Tools::to_value(name, id_names);
2798 bool found_unrecognized_tag = false;
2799 bool found_boundary_id = false;
2800 // If the above did not throw, we keep going, and retrieve
2801 // all the information that we know how to translate.
2802 for (const auto &it : id_names)
2803 {
2804 const auto &name = it.first;
2805 const auto &id = it.second;
2806 if (entity_dim == dim && name == "MaterialID")
2807 {
2808 boundary_id = static_cast<types::boundary_id>(id);
2809 found_boundary_id = true;
2810 }
2811 else if (entity_dim < dim && name == "BoundaryID")
2812 {
2813 boundary_id = static_cast<types::boundary_id>(id);
2814 found_boundary_id = true;
2815 }
2816 else if (name == "ManifoldID")
2817 manifold_id = static_cast<types::manifold_id>(id);
2818 else
2819 // We did not recognize one of the keys. We'll fall
2820 // back to setting the boundary id to the physical tag
2821 // after reading all strings.
2822 found_unrecognized_tag = true;
2823 }
2824 // If we didn't find a BoundaryID:XX or MaterialID:XX, and
2825 // something was found but not recognized, then we set the
2826 // material id or boundary using the physical tag directly.
2827 if (found_unrecognized_tag && !found_boundary_id)
2828 boundary_id = physical_tag;
2829 }
2830 catch (...)
2831 {
2832 // When the above didn't work, we revert to the old
2833 // behaviour: the physical tag itself is interpreted either
2834 // as a material_id or a boundary_id, and no manifold id is
2835 // known
2836 boundary_id = physical_tag;
2837 }
2838 }
2839 }
2840
2841 // Get the mesh elements for the entity (dim, tag):
2842 std::vector<int> element_types;
2843 std::vector<std::vector<std::size_t>> element_ids, element_nodes;
2844 gmsh::model::mesh::getElements(
2845 element_types, element_ids, element_nodes, entity_dim, entity_tag);
2846
2847 for (unsigned int i = 0; i < element_types.size(); ++i)
2848 {
2849 const auto &type = gmsh_to_dealii_type.at(element_types[i]);
2850 const auto n_vertices = gmsh_to_dealii[type].size();
2851 const auto &elements = element_ids[i];
2852 const auto &nodes = element_nodes[i];
2853 for (unsigned int j = 0; j < elements.size(); ++j)
2854 {
2855 if (entity_dim == dim)
2856 {
2857 cells.emplace_back(n_vertices);
2858 auto &cell = cells.back();
2859 for (unsigned int v = 0; v < n_vertices; ++v)
2860 {
2861 cell.vertices[v] =
2862 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2863 if (dim == 1)
2864 vertex_counts[cell.vertices[v]] += 1u;
2865 }
2866 cell.manifold_id = manifold_id;
2867 cell.material_id = boundary_id;
2868 }
2869 else if (entity_dim == 2)
2870 {
2871 subcelldata.boundary_quads.emplace_back(n_vertices);
2872 auto &face = subcelldata.boundary_quads.back();
2873 for (unsigned int v = 0; v < n_vertices; ++v)
2874 face.vertices[v] =
2875 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2876
2877 face.manifold_id = manifold_id;
2878 face.boundary_id = boundary_id;
2879 }
2880 else if (entity_dim == 1)
2881 {
2882 subcelldata.boundary_lines.emplace_back(n_vertices);
2883 auto &line = subcelldata.boundary_lines.back();
2884 for (unsigned int v = 0; v < n_vertices; ++v)
2885 line.vertices[v] =
2886 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2887
2888 line.manifold_id = manifold_id;
2889 line.boundary_id = boundary_id;
2890 }
2891 else if (entity_dim == 0)
2892 {
2893 // This should only happen in one dimension.
2894 AssertDimension(dim, 1);
2895 for (unsigned int j = 0; j < elements.size(); ++j)
2896 boundary_ids_1d[nodes[j] - 1] = boundary_id;
2897 }
2898 }
2899 }
2900 }
2901
2902 // apply_grid_fixup_functions() may invalidate the vertex indices (since it
2903 // will delete duplicated or unused vertices). Get around this by storing
2904 // Points directly in that case so that the comparisons are valid.
2905 std::vector<std::pair<Point<spacedim>, types::boundary_id>> boundary_id_pairs;
2906 if (dim == 1)
2907 for (const auto &pair : vertex_counts)
2908 if (pair.second == 1u)
2909 boundary_id_pairs.emplace_back(vertices[pair.first],
2910 boundary_ids_1d[pair.first]);
2911
2912 apply_grid_fixup_functions(vertices, cells, subcelldata);
2913 tria->create_triangulation(vertices, cells, subcelldata);
2914
2915 // in 1d, we also have to attach boundary ids to vertices, which does not
2916 // currently work through the call above.
2917 if (dim == 1)
2918 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2919
2920 gmsh::clear();
2921 gmsh::finalize();
2922}
2923#endif
2924
2925
2926
2927template <int dim, int spacedim>
2928void
2930 std::string &header,
2931 std::vector<unsigned int> &tecplot2deal,
2932 unsigned int &n_vars,
2933 unsigned int &n_vertices,
2934 unsigned int &n_cells,
2935 std::vector<unsigned int> &IJK,
2936 bool &structured,
2937 bool &blocked)
2938{
2939 Assert(tecplot2deal.size() == dim, ExcInternalError());
2940 Assert(IJK.size() == dim, ExcInternalError());
2941 // initialize the output variables
2942 n_vars = 0;
2943 n_vertices = 0;
2944 n_cells = 0;
2945 switch (dim)
2946 {
2947 case 3:
2948 IJK[2] = 0;
2950 case 2:
2951 IJK[1] = 0;
2953 case 1:
2954 IJK[0] = 0;
2955 }
2956 structured = true;
2957 blocked = false;
2958
2959 // convert the string to upper case
2960 std::transform(header.begin(), header.end(), header.begin(), ::toupper);
2961
2962 // replace all tabs, commas, newlines by
2963 // whitespaces
2964 std::replace(header.begin(), header.end(), '\t', ' ');
2965 std::replace(header.begin(), header.end(), ',', ' ');
2966 std::replace(header.begin(), header.end(), '\n', ' ');
2967
2968 // now remove whitespace in front of and
2969 // after '='
2970 std::string::size_type pos = header.find('=');
2971
2972 while (pos != static_cast<std::string::size_type>(std::string::npos))
2973 if (header[pos + 1] == ' ')
2974 header.erase(pos + 1, 1);
2975 else if (header[pos - 1] == ' ')
2976 {
2977 header.erase(pos - 1, 1);
2978 --pos;
2979 }
2980 else
2981 pos = header.find('=', ++pos);
2982
2983 // split the string into individual entries
2984 std::vector<std::string> entries =
2985 Utilities::break_text_into_lines(header, 1, ' ');
2986
2987 // now go through the list and try to extract
2988 for (unsigned int i = 0; i < entries.size(); ++i)
2989 {
2990 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\""))
2991 {
2992 ++n_vars;
2993 // we assume, that the first variable
2994 // is x or no coordinate at all (not y or z)
2995 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\"X\""))
2996 {
2997 tecplot2deal[0] = 0;
2998 }
2999 ++i;
3000 while (entries[i][0] == '"')
3001 {
3002 if (entries[i] == "\"X\"")
3003 tecplot2deal[0] = n_vars;
3004 else if (entries[i] == "\"Y\"")
3005 {
3006 // we assume, that y contains
3007 // zero data in 1d, so do
3008 // nothing
3009 if (dim > 1)
3010 tecplot2deal[1] = n_vars;
3011 }
3012 else if (entries[i] == "\"Z\"")
3013 {
3014 // we assume, that z contains
3015 // zero data in 1d and 2d, so
3016 // do nothing
3017 if (dim > 2)
3018 tecplot2deal[2] = n_vars;
3019 }
3020 ++n_vars;
3021 ++i;
3022 }
3023 // set i back, so that the next
3024 // string is treated correctly
3025 --i;
3026
3028 n_vars >= dim,
3029 ExcMessage(
3030 "Tecplot file must contain at least one variable for each dimension"));
3031 for (unsigned int d = 1; d < dim; ++d)
3033 tecplot2deal[d] > 0,
3034 ExcMessage(
3035 "Tecplot file must contain at least one variable for each dimension."));
3036 }
3037 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE=ORDERED"))
3038 structured = true;
3039 else if (Utilities::match_at_string_start(entries[i],
3040 "ZONETYPE=FELINESEG") &&
3041 dim == 1)
3042 structured = false;
3043 else if (Utilities::match_at_string_start(entries[i],
3044 "ZONETYPE=FEQUADRILATERAL") &&
3045 dim == 2)
3046 structured = false;
3047 else if (Utilities::match_at_string_start(entries[i],
3048 "ZONETYPE=FEBRICK") &&
3049 dim == 3)
3050 structured = false;
3051 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE="))
3052 // unsupported ZONETYPE
3053 {
3054 AssertThrow(false,
3055 ExcMessage(
3056 "The tecplot file contains an unsupported ZONETYPE."));
3057 }
3058 else if (Utilities::match_at_string_start(entries[i],
3059 "DATAPACKING=POINT"))
3060 blocked = false;
3061 else if (Utilities::match_at_string_start(entries[i],
3062 "DATAPACKING=BLOCK"))
3063 blocked = true;
3064 else if (Utilities::match_at_string_start(entries[i], "F=POINT"))
3065 {
3066 structured = true;
3067 blocked = false;
3068 }
3069 else if (Utilities::match_at_string_start(entries[i], "F=BLOCK"))
3070 {
3071 structured = true;
3072 blocked = true;
3073 }
3074 else if (Utilities::match_at_string_start(entries[i], "F=FEPOINT"))
3075 {
3076 structured = false;
3077 blocked = false;
3078 }
3079 else if (Utilities::match_at_string_start(entries[i], "F=FEBLOCK"))
3080 {
3081 structured = false;
3082 blocked = true;
3083 }
3084 else if (Utilities::match_at_string_start(entries[i],
3085 "ET=QUADRILATERAL") &&
3086 dim == 2)
3087 structured = false;
3088 else if (Utilities::match_at_string_start(entries[i], "ET=BRICK") &&
3089 dim == 3)
3090 structured = false;
3091 else if (Utilities::match_at_string_start(entries[i], "ET="))
3092 // unsupported ElementType
3093 {
3095 false,
3096 ExcMessage(
3097 "The tecplot file contains an unsupported ElementType."));
3098 }
3099 else if (Utilities::match_at_string_start(entries[i], "I="))
3100 IJK[0] = Utilities::get_integer_at_position(entries[i], 2).first;
3101 else if (Utilities::match_at_string_start(entries[i], "J="))
3102 {
3103 IJK[1] = Utilities::get_integer_at_position(entries[i], 2).first;
3105 dim > 1 || IJK[1] == 1,
3106 ExcMessage(
3107 "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
3108 }
3109 else if (Utilities::match_at_string_start(entries[i], "K="))
3110 {
3111 IJK[2] = Utilities::get_integer_at_position(entries[i], 2).first;
3113 dim > 2 || IJK[2] == 1,
3114 ExcMessage(
3115 "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
3116 }
3117 else if (Utilities::match_at_string_start(entries[i], "N="))
3118 n_vertices = Utilities::get_integer_at_position(entries[i], 2).first;
3119 else if (Utilities::match_at_string_start(entries[i], "E="))
3120 n_cells = Utilities::get_integer_at_position(entries[i], 2).first;
3121 }
3122
3123 // now we have read all the fields we are
3124 // interested in. do some checks and
3125 // calculate the variables
3126 if (structured)
3127 {
3128 n_vertices = 1;
3129 n_cells = 1;
3130 for (unsigned int d = 0; d < dim; ++d)
3131 {
3133 IJK[d] > 0,
3134 ExcMessage(
3135 "Tecplot file does not contain a complete and consistent set of parameters"));
3136 n_vertices *= IJK[d];
3137 n_cells *= (IJK[d] - 1);
3138 }
3139 }
3140 else
3141 {
3143 n_vertices > 0,
3144 ExcMessage(
3145 "Tecplot file does not contain a complete and consistent set of parameters"));
3146 if (n_cells == 0)
3147 // this means an error, although
3148 // tecplot itself accepts entries like
3149 // 'J=20' instead of 'E=20'. therefore,
3150 // take the max of IJK
3151 n_cells = *std::max_element(IJK.begin(), IJK.end());
3153 n_cells > 0,
3154 ExcMessage(
3155 "Tecplot file does not contain a complete and consistent set of parameters"));
3156 }
3157}
3158
3159
3160
3161template <>
3162void
3164{
3165 const unsigned int dim = 2;
3166 const unsigned int spacedim = 2;
3167 Assert(tria != nullptr, ExcNoTriangulationSelected());
3168 AssertThrow(in.fail() == false, ExcIO());
3169
3170 // skip comments at start of file
3171 skip_comment_lines(in, '#');
3172
3173 // some strings for parsing the header
3174 std::string line, header;
3175
3176 // first, concatenate all header lines
3177 // create a searchstring with almost all
3178 // letters. exclude e and E from the letters
3179 // to search, as they might appear in
3180 // exponential notation
3181 std::string letters = "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
3182
3183 getline(in, line);
3184 while (line.find_first_of(letters) != std::string::npos)
3185 {
3186 header += " " + line;
3187 getline(in, line);
3188 }
3189
3190 // now create some variables holding
3191 // important information on the mesh, get
3192 // this information from the header string
3193 std::vector<unsigned int> tecplot2deal(dim);
3194 std::vector<unsigned int> IJK(dim);
3195 unsigned int n_vars, n_vertices, n_cells;
3196 bool structured, blocked;
3197
3198 parse_tecplot_header(header,
3199 tecplot2deal,
3200 n_vars,
3201 n_vertices,
3202 n_cells,
3203 IJK,
3204 structured,
3205 blocked);
3206
3207 // reserve space for vertices. note, that in
3208 // tecplot vertices are ordered beginning
3209 // with 1, whereas in deal all indices start
3210 // with 0. in order not to use -1 for all the
3211 // connectivity information, a 0th vertex
3212 // (unused) is inserted at the origin.
3213 std::vector<Point<spacedim>> vertices(n_vertices + 1);
3215 // reserve space for cells
3216 std::vector<CellData<dim>> cells(n_cells);
3217 SubCellData subcelldata;
3218
3219 if (blocked)
3220 {
3221 // blocked data format. first we get all
3222 // the values of the first variable for
3223 // all points, after that we get all
3224 // values for the second variable and so
3225 // on.
3226
3227 // dummy variable to read in all the info
3228 // we do not want to use
3229 double dummy;
3230 // which is the first index to read in
3231 // the loop (see below)
3232 unsigned int next_index = 0;
3233
3234 // note, that we have already read the
3235 // first line containing the first variable
3236 if (tecplot2deal[0] == 0)
3237 {
3238 // we need the information in this
3239 // line, so extract it
3240 std::vector<std::string> first_var =
3242 char *endptr;
3243 for (unsigned int i = 1; i < first_var.size() + 1; ++i)
3244 vertices[i][0] = std::strtod(first_var[i - 1].c_str(), &endptr);
3245
3246 // if there are many points, the data
3247 // for this var might continue in the
3248 // next line(s)
3249 for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
3250 in >> vertices[j][next_index];
3251 // now we got all values of the first
3252 // variable, so increase the counter
3253 next_index = 1;
3254 }
3255
3256 // main loop over all variables
3257 for (unsigned int i = 1; i < n_vars; ++i)
3258 {
3259 // if we read all the important
3260 // variables and do not want to
3261 // read further, because we are
3262 // using a structured grid, we can
3263 // stop here (and skip, for
3264 // example, a whole lot of solution
3265 // variables)
3266 if (next_index == dim && structured)
3267 break;
3268
3269 if ((next_index < dim) && (i == tecplot2deal[next_index]))
3270 {
3271 // we need this line, read it in
3272 for (unsigned int j = 1; j < n_vertices + 1; ++j)
3273 in >> vertices[j][next_index];
3274 ++next_index;
3275 }
3276 else
3277 {
3278 // we do not need this line, read
3279 // it in and discard it
3280 for (unsigned int j = 1; j < n_vertices + 1; ++j)
3281 in >> dummy;
3282 }
3283 }
3284 Assert(next_index == dim, ExcInternalError());
3285 }
3286 else
3287 {
3288 // the data is not blocked, so we get all
3289 // the variables for one point, then the
3290 // next and so on. create a vector to
3291 // hold these components
3292 std::vector<double> vars(n_vars);
3293
3294 // now fill the first vertex. note, that we
3295 // have already read the first line
3296 // containing the first vertex
3297 std::vector<std::string> first_vertex =
3299 char *endptr;
3300 for (unsigned int d = 0; d < dim; ++d)
3301 vertices[1][d] =
3302 std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
3303
3304 // read the remaining vertices from the
3305 // list
3306 for (unsigned int v = 2; v < n_vertices + 1; ++v)
3307 {
3308 for (unsigned int i = 0; i < n_vars; ++i)
3309 in >> vars[i];
3310 // fill the vertex
3311 // coordinates. respect the position
3312 // of coordinates in the list of
3313 // variables
3314 for (unsigned int i = 0; i < dim; ++i)
3315 vertices[v][i] = vars[tecplot2deal[i]];
3316 }
3317 }
3318
3319 if (structured)
3320 {
3321 // this is the part of the code that only
3322 // works in 2d
3323 unsigned int I = IJK[0], J = IJK[1];
3324
3325 unsigned int cell = 0;
3326 // set up array of cells
3327 for (unsigned int j = 0; j < J - 1; ++j)
3328 for (unsigned int i = 1; i < I; ++i)
3329 {
3330 cells[cell].vertices[0] = i + j * I;
3331 cells[cell].vertices[1] = i + 1 + j * I;
3332 cells[cell].vertices[2] = i + (j + 1) * I;
3333 cells[cell].vertices[3] = i + 1 + (j + 1) * I;
3334 ++cell;
3335 }
3336 Assert(cell == n_cells, ExcInternalError());
3337 std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
3338 unsigned int k = 0;
3339 for (unsigned int i = 1; i < I + 1; ++i)
3340 {
3341 boundary_vertices[k] = i;
3342 ++k;
3343 boundary_vertices[k] = i + (J - 1) * I;
3344 ++k;
3345 }
3346 for (unsigned int j = 1; j < J - 1; ++j)
3347 {
3348 boundary_vertices[k] = 1 + j * I;
3349 ++k;
3350 boundary_vertices[k] = I + j * I;
3351 ++k;
3352 }
3353 Assert(k == boundary_vertices.size(), ExcInternalError());
3354 // delete the duplicated vertices at the
3355 // boundary, which occur, e.g. in c-type
3356 // or o-type grids around a body
3357 // (airfoil). this automatically deletes
3358 // unused vertices as well.
3360 cells,
3361 subcelldata,
3362 boundary_vertices);
3363 }
3364 else
3365 {
3366 // set up array of cells, unstructured
3367 // mode, so the connectivity is
3368 // explicitly given
3369 for (unsigned int i = 0; i < n_cells; ++i)
3370 {
3371 // note that since in the input file
3372 // we found the number of cells at
3373 // the top, there should still be
3374 // input here, so check this:
3375 AssertThrow(in.fail() == false, ExcIO());
3376
3377 // get the connectivity from the
3378 // input file. the vertices are
3379 // ordered like in the ucd format
3380 for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
3381 in >> cells[i].vertices[GeometryInfo<dim>::ucd_to_deal[j]];
3382 }
3383 }
3384 AssertThrow(in.fail() == false, ExcIO());
3385
3386 apply_grid_fixup_functions(vertices, cells, subcelldata);
3387 tria->create_triangulation(vertices, cells, subcelldata);
3388}
3389
3390
3391
3392template <int dim, int spacedim>
3393void
3398
3399
3400
3401template <int dim, int spacedim>
3402void
3403GridIn<dim, spacedim>::read_assimp(const std::string &filename,
3404 const unsigned int mesh_index,
3405 const bool remove_duplicates,
3406 const double tol,
3407 const bool ignore_unsupported_types)
3408{
3409#ifdef DEAL_II_WITH_ASSIMP
3410 // Only good for surface grids.
3411 AssertThrow(dim < 3, ExcImpossibleInDim(dim));
3412
3413 // Create an instance of the Importer class
3414 Assimp::Importer importer;
3415
3416 // And have it read the given file with some postprocessing
3417 const aiScene *scene =
3418 importer.ReadFile(filename.c_str(),
3419 aiProcess_RemoveComponent |
3420 aiProcess_JoinIdenticalVertices |
3421 aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
3422 aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
3423
3424 // If the import failed, report it
3425 AssertThrow(scene != nullptr, ExcMessage(importer.GetErrorString()));
3426
3427 AssertThrow(scene->mNumMeshes != 0,
3428 ExcMessage("Input file contains no meshes."));
3429
3431 (mesh_index < scene->mNumMeshes),
3432 ExcMessage("Too few meshes in the file."));
3433
3434 unsigned int start_mesh =
3435 (mesh_index == numbers::invalid_unsigned_int ? 0 : mesh_index);
3436 unsigned int end_mesh =
3437 (mesh_index == numbers::invalid_unsigned_int ? scene->mNumMeshes :
3438 mesh_index + 1);
3439
3440 // Deal.II objects are created empty, and then filled with imported file.
3441 std::vector<Point<spacedim>> vertices;
3442 std::vector<CellData<dim>> cells;
3443 SubCellData subcelldata;
3444
3445 // A series of counters to merge cells.
3446 unsigned int v_offset = 0;
3447 unsigned int c_offset = 0;
3448
3449 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
3450 {0, 1, 5, 4, 2, 3, 7, 6}};
3451 // The index of the mesh will be used as a material index.
3452 for (unsigned int m = start_mesh; m < end_mesh; ++m)
3453 {
3454 const aiMesh *mesh = scene->mMeshes[m];
3455
3456 // Check that we know what to do with this mesh, otherwise just
3457 // ignore it
3458 if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
3459 {
3460 AssertThrow(ignore_unsupported_types,
3461 ExcMessage("Incompatible mesh " + std::to_string(m) +
3462 "/" + std::to_string(scene->mNumMeshes)));
3463 continue;
3464 }
3465 else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
3466 {
3467 AssertThrow(ignore_unsupported_types,
3468 ExcMessage("Incompatible mesh " + std::to_string(m) +
3469 "/" + std::to_string(scene->mNumMeshes)));
3470 continue;
3471 }
3472 // Vertices
3473 const unsigned int n_vertices = mesh->mNumVertices;
3474 const aiVector3D *mVertices = mesh->mVertices;
3475
3476 // Faces
3477 const unsigned int n_faces = mesh->mNumFaces;
3478 const aiFace *mFaces = mesh->mFaces;
3479
3480 vertices.resize(v_offset + n_vertices);
3481 cells.resize(c_offset + n_faces);
3482
3483 for (unsigned int i = 0; i < n_vertices; ++i)
3484 for (unsigned int d = 0; d < spacedim; ++d)
3485 vertices[i + v_offset][d] = mVertices[i][d];
3486
3487 unsigned int valid_cell = c_offset;
3488 for (unsigned int i = 0; i < n_faces; ++i)
3489 {
3490 if (mFaces[i].mNumIndices == GeometryInfo<dim>::vertices_per_cell)
3491 {
3492 for (const unsigned int f : GeometryInfo<dim>::vertex_indices())
3493 {
3494 cells[valid_cell]
3495 .vertices[dim == 3 ? local_vertex_numbering[f] :
3497 mFaces[i].mIndices[f] + v_offset;
3498 }
3499 cells[valid_cell].material_id = m;
3500 ++valid_cell;
3501 }
3502 else
3503 {
3504 AssertThrow(ignore_unsupported_types,
3505 ExcMessage("Face " + std::to_string(i) + " of mesh " +
3506 std::to_string(m) + " has " +
3507 std::to_string(mFaces[i].mNumIndices) +
3508 " vertices. We expected only " +
3509 std::to_string(
3511 }
3512 }
3513 cells.resize(valid_cell);
3514
3515 // The vertices are added all at once. Cells are checked for
3516 // validity, so only valid_cells are now present in the deal.II
3517 // list of cells.
3518 v_offset += n_vertices;
3519 c_offset = valid_cell;
3520 }
3521
3522 // No cells were read
3523 if (cells.empty())
3524 return;
3525
3526 if (remove_duplicates)
3527 {
3528 // The function delete_duplicated_vertices() needs to be called more
3529 // than once if a vertex is duplicated more than once. So we keep
3530 // calling it until the number of vertices does not change any more.
3531 unsigned int n_verts = 0;
3532 while (n_verts != vertices.size())
3533 {
3534 n_verts = vertices.size();
3535 std::vector<unsigned int> considered_vertices;
3537 vertices, cells, subcelldata, considered_vertices, tol);
3538 }
3539 }
3540
3541 apply_grid_fixup_functions(vertices, cells, subcelldata);
3542 tria->create_triangulation(vertices, cells, subcelldata);
3543
3544#else
3545 (void)filename;
3546 (void)mesh_index;
3547 (void)remove_duplicates;
3548 (void)tol;
3549 (void)ignore_unsupported_types;
3550 AssertThrow(false, ExcNeedsAssimp());
3551#endif
3552}
3553
3554#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3555// Namespace containing some extra functions for reading ExodusII files
3556namespace
3557{
3558 // Convert ExodusII strings to cell types. Use the number of nodes per
3559 // element to disambiguate some cases.
3561 exodusii_name_to_type(const std::string &type_name,
3562 const int n_nodes_per_element)
3563 {
3564 Assert(type_name.size() > 0, ExcInternalError());
3565 // Try to canonify the name by switching to upper case and removing
3566 // trailing numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and
3567 // PYRAMID13 all equal.
3568 std::string type_name_2 = type_name;
3569 std::transform(type_name_2.begin(),
3570 type_name_2.end(),
3571 type_name_2.begin(),
3572 [](unsigned char c) { return std::toupper(c); });
3573 const std::string numbers = "0123456789";
3574 type_name_2.erase(std::find_first_of(type_name_2.begin(),
3575 type_name_2.end(),
3576 numbers.begin(),
3577 numbers.end()),
3578 type_name_2.end());
3579
3580 // The manual specifies BAR, BEAM, and TRUSS: in practice people use EDGE
3581 if (type_name_2 == "BAR" || type_name_2 == "BEAM" ||
3582 type_name_2 == "EDGE" || type_name_2 == "TRUSS")
3583 return ReferenceCells::Line;
3584 else if (type_name_2 == "TRI" || type_name_2 == "TRIANGLE")
3586 else if (type_name_2 == "QUAD" || type_name_2 == "QUADRILATERAL")
3588 else if (type_name_2 == "SHELL")
3589 {
3590 if (n_nodes_per_element == 3)
3592 else
3594 }
3595 else if (type_name_2 == "TET" || type_name_2 == "TETRA" ||
3596 type_name_2 == "TETRAHEDRON")
3598 else if (type_name_2 == "PYRA" || type_name_2 == "PYRAMID")
3600 else if (type_name_2 == "WEDGE")
3601 return ReferenceCells::Wedge;
3602 else if (type_name_2 == "HEX" || type_name_2 == "HEXAHEDRON")
3604
3607 }
3608
3609 // Associate deal.II boundary ids with sidesets (a face can be in multiple
3610 // sidesets - to translate we assign each set of side set ids to a
3611 // boundary_id or manifold_id)
3612 template <int dim, int spacedim = dim>
3613 std::pair<SubCellData, std::vector<std::vector<int>>>
3614 read_exodusii_sidesets(const int ex_id,
3615 const int n_side_sets,
3616 const std::vector<CellData<dim>> &cells,
3617 const bool apply_all_indicators_to_manifolds)
3618 {
3619 SubCellData subcelldata;
3620 std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
3621 // boundary id 0 is the default
3622 b_or_m_id_to_sideset_ids.emplace_back();
3623 // deal.II does not support assigning boundary ids with nonzero
3624 // codimension meshes so completely skip this information in that case.
3625 //
3626 // Exodus prints warnings if we try to get empty sets so always check
3627 // first
3628 if (dim == spacedim && n_side_sets > 0)
3629 {
3630 std::vector<int> side_set_ids(n_side_sets);
3631 int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
3632 AssertThrowExodusII(ierr);
3633
3634 // First collect all side sets on all boundary faces (indexed here as
3635 // max_faces_per_cell * cell_n + face_n). We then sort and uniquify
3636 // the side sets so that we can convert a set of side set indices into
3637 // a single deal.II boundary or manifold id (and save the
3638 // correspondence).
3639 constexpr auto max_faces_per_cell = GeometryInfo<dim>::faces_per_cell;
3640 std::map<std::size_t, std::vector<int>> face_side_sets;
3641 for (const int side_set_id : side_set_ids)
3642 {
3643 int n_sides = -1;
3644 int n_distribution_factors = -1;
3645
3646 ierr = ex_get_set_param(ex_id,
3647 EX_SIDE_SET,
3648 side_set_id,
3649 &n_sides,
3650 &n_distribution_factors);
3651 AssertThrowExodusII(ierr);
3652 if (n_sides > 0)
3653 {
3654 std::vector<int> elements(n_sides);
3655 std::vector<int> faces(n_sides);
3656 ierr = ex_get_set(ex_id,
3657 EX_SIDE_SET,
3658 side_set_id,
3659 elements.data(),
3660 faces.data());
3661 AssertThrowExodusII(ierr);
3662
3663 // According to the manual (subsection 4.8): "The internal
3664 // number of an element numbering is defined implicitly by the
3665 // order in which it appears in the file. Elements are
3666 // numbered internally (beginning with 1) consecutively across
3667 // all element blocks." Hence element i in Exodus numbering is
3668 // entry i - 1 in the cells array.
3669 for (int side_n = 0; side_n < n_sides; ++side_n)
3670 {
3671 const long element_n = elements[side_n] - 1;
3672 const long face_n = faces[side_n] - 1;
3673 const std::size_t face_id =
3674 element_n * max_faces_per_cell + face_n;
3675 face_side_sets[face_id].push_back(side_set_id);
3676 }
3677 }
3678 }
3679
3680 // Collect into a sortable data structure:
3681 std::vector<std::pair<std::size_t, std::vector<int>>>
3682 face_id_to_side_sets;
3683 face_id_to_side_sets.reserve(face_side_sets.size());
3684 for (auto &pair : face_side_sets)
3685 {
3686 Assert(pair.second.size() > 0, ExcInternalError());
3687 face_id_to_side_sets.push_back(std::move(pair));
3688 }
3689
3690 // sort by side sets:
3691 std::sort(face_id_to_side_sets.begin(),
3692 face_id_to_side_sets.end(),
3693 [](const auto &a, const auto &b) {
3694 return std::lexicographical_compare(a.second.begin(),
3695 a.second.end(),
3696 b.second.begin(),
3697 b.second.end());
3698 });
3699
3700 if (dim == 2)
3701 subcelldata.boundary_lines.reserve(face_id_to_side_sets.size());
3702 else if (dim == 3)
3703 subcelldata.boundary_quads.reserve(face_id_to_side_sets.size());
3704 types::boundary_id current_b_or_m_id = 0;
3705 for (const auto &pair : face_id_to_side_sets)
3706 {
3707 const std::size_t face_id = pair.first;
3708 const std::vector<int> &face_sideset_ids = pair.second;
3709 if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
3710 {
3711 // Since we sorted by sideset ids we are guaranteed that if
3712 // this doesn't match the last set then it has not yet been
3713 // seen
3714 ++current_b_or_m_id;
3715 b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
3716 Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
3718 }
3719 // Record the b_or_m_id of the current face.
3720 const unsigned int local_face_n = face_id % max_faces_per_cell;
3721 const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
3722 const ReferenceCell cell_type =
3724 const unsigned int deal_face_n =
3725 cell_type.exodusii_face_to_deal_face(local_face_n);
3726 const ReferenceCell face_reference_cell =
3727 cell_type.face_reference_cell(deal_face_n);
3728
3729 // The orientation we pick doesn't matter here since when we
3730 // create the Triangulation we will sort the vertices for each
3731 // CellData object created here.
3732 if (dim == 2)
3733 {
3734 CellData<1> boundary_line(face_reference_cell.n_vertices());
3735 if (apply_all_indicators_to_manifolds)
3736 boundary_line.manifold_id = current_b_or_m_id;
3737 else
3738 boundary_line.boundary_id = current_b_or_m_id;
3739 for (unsigned int j = 0; j < face_reference_cell.n_vertices();
3740 ++j)
3741 boundary_line.vertices[j] =
3742 cell.vertices[cell_type.face_to_cell_vertices(
3743 deal_face_n, j, 0)];
3744
3745 subcelldata.boundary_lines.push_back(std::move(boundary_line));
3746 }
3747 else if (dim == 3)
3748 {
3749 CellData<2> boundary_quad(face_reference_cell.n_vertices());
3750 if (apply_all_indicators_to_manifolds)
3751 boundary_quad.manifold_id = current_b_or_m_id;
3752 else
3753 boundary_quad.boundary_id = current_b_or_m_id;
3754 for (unsigned int j = 0; j < face_reference_cell.n_vertices();
3755 ++j)
3756 boundary_quad.vertices[j] =
3757 cell.vertices[cell_type.face_to_cell_vertices(
3758 deal_face_n, j, 0)];
3759
3760 subcelldata.boundary_quads.push_back(std::move(boundary_quad));
3761 }
3762 }
3763 }
3764
3765 return std::make_pair(std::move(subcelldata),
3766 std::move(b_or_m_id_to_sideset_ids));
3767 }
3768} // namespace
3769#endif
3770
3771template <int dim, int spacedim>
3774 const std::string &filename,
3775 const bool apply_all_indicators_to_manifolds)
3776{
3777#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3778 // deal.II always uses double precision numbers for geometry
3779 int component_word_size = sizeof(double);
3780 // setting to zero uses the stored word size
3781 int floating_point_word_size = 0;
3782 float ex_version = 0.0;
3783
3784 const int ex_id = ex_open(filename.c_str(),
3785 EX_READ,
3786 &component_word_size,
3787 &floating_point_word_size,
3788 &ex_version);
3789 AssertThrow(ex_id > 0,
3790 ExcMessage("ExodusII failed to open the specified input file."));
3791
3792 // Read basic mesh information:
3793 std::vector<char> string_temp(MAX_LINE_LENGTH + 1, '\0');
3794 int mesh_dimension = 0;
3795 int n_nodes = 0;
3796 int n_elements = 0;
3797 int n_element_blocks = 0;
3798 int n_node_sets = 0;
3799 int n_side_sets = 0;
3800
3801 int ierr = ex_get_init(ex_id,
3802 string_temp.data(),
3803 &mesh_dimension,
3804 &n_nodes,
3805 &n_elements,
3806 &n_element_blocks,
3807 &n_node_sets,
3808 &n_side_sets);
3809 AssertThrowExodusII(ierr);
3810 AssertDimension(mesh_dimension, spacedim);
3811
3812 // Read nodes:
3813 //
3814 // Even if there is a node numbering array the values stored inside the
3815 // ExodusII file must use the contiguous, internal ordering (see Section 4.5
3816 // of the manual - "Internal (contiguously numbered) node and element IDs
3817 // must be used for all data structures that contain node or element numbers
3818 // (IDs), including node set node lists, side set element lists, and element
3819 // connectivity.")
3820 std::vector<Point<spacedim>> vertices;
3821 vertices.reserve(n_nodes);
3822 {
3823 std::vector<double> xs(n_nodes);
3824 std::vector<double> ys(n_nodes);
3825 std::vector<double> zs(n_nodes);
3826
3827 ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
3828 AssertThrowExodusII(ierr);
3829
3830 for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
3831 {
3832 switch (spacedim)
3833 {
3834 case 1:
3835 vertices.emplace_back(xs[vertex_n]);
3836 break;
3837 case 2:
3838 vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
3839 break;
3840 case 3:
3841 vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
3842 break;
3843 default:
3844 Assert(spacedim <= 3, ExcNotImplemented());
3845 }
3846 }
3847 }
3848
3849 std::vector<int> element_block_ids(n_element_blocks);
3850 ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
3851 AssertThrowExodusII(ierr);
3852
3853 std::vector<CellData<dim>> cells;
3854 cells.reserve(n_elements);
3855 // Elements are grouped together by same reference cell type in element
3856 // blocks. There may be multiple blocks for a single reference cell type,
3857 // but "each element block may contain only one element type".
3858 for (const int element_block_id : element_block_ids)
3859 {
3860 std::fill(string_temp.begin(), string_temp.end(), '\0');
3861 int n_block_elements = 0;
3862 int n_nodes_per_element = 0;
3863 int n_edges_per_element = 0;
3864 int n_faces_per_element = 0;
3865 int n_attributes_per_element = 0;
3866
3867 // Extract element data.
3868 ierr = ex_get_block(ex_id,
3869 EX_ELEM_BLOCK,
3870 element_block_id,
3871 string_temp.data(),
3872 &n_block_elements,
3873 &n_nodes_per_element,
3874 &n_edges_per_element,
3875 &n_faces_per_element,
3876 &n_attributes_per_element);
3877 AssertThrowExodusII(ierr);
3878 const ReferenceCell type =
3879 exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
3880 AssertThrow(type.get_dimension() == dim,
3881 ExcMessage(
3882 "The ExodusII block " + std::to_string(element_block_id) +
3883 " with element type " + std::string(string_temp.data()) +
3884 " has dimension " + std::to_string(type.get_dimension()) +
3885 ", which does not match the topological mesh dimension " +
3886 std::to_string(dim) + "."));
3887
3888 // The number of nodes per element may be larger than what we want to
3889 // read - for example, if the Exodus file contains a QUAD9 element, we
3890 // only want to read the first four values and ignore the rest.
3891 Assert(int(type.n_vertices()) <= n_nodes_per_element, ExcInternalError());
3892
3893 std::vector<int> connection(n_nodes_per_element * n_block_elements);
3894 ierr = ex_get_conn(ex_id,
3895 EX_ELEM_BLOCK,
3896 element_block_id,
3897 connection.data(),
3898 nullptr,
3899 nullptr);
3900 AssertThrowExodusII(ierr);
3901
3902 for (unsigned int elem_n = 0; elem_n < connection.size();
3903 elem_n += n_nodes_per_element)
3904 {
3905 CellData<dim> cell(type.n_vertices());
3906 for (const unsigned int i : type.vertex_indices())
3907 {
3909 connection[elem_n + i] - 1;
3910 }
3911 cell.material_id = element_block_id;
3912 cells.push_back(std::move(cell));
3913 }
3914 }
3915
3916 // Extract boundary data.
3917 auto pair = read_exodusii_sidesets<dim, spacedim>(
3918 ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
3919 ierr = ex_close(ex_id);
3920 AssertThrowExodusII(ierr);
3921
3922 apply_grid_fixup_functions(vertices, cells, pair.first);
3923 tria->create_triangulation(vertices, cells, pair.first);
3924 ExodusIIData out;
3925 out.id_to_sideset_ids = std::move(pair.second);
3926 return out;
3927#else
3928 (void)filename;
3929 (void)apply_all_indicators_to_manifolds;
3930 AssertThrow(false, ExcNeedsExodusII());
3931 return {};
3932#endif
3933}
3934
3935
3936template <int dim, int spacedim>
3937void
3939{
3940 std::string line;
3941 while (in)
3942 {
3943 // get line
3944 getline(in, line);
3945
3946 // check if this is a line that
3947 // consists only of spaces, and
3948 // if not put the whole thing
3949 // back and return
3950 if (std::find_if(line.begin(), line.end(), [](const char c) {
3951 return c != ' ';
3952 }) != line.end())
3953 {
3954 in.putback('\n');
3955 for (int i = line.size() - 1; i >= 0; --i)
3956 in.putback(line[i]);
3957 return;
3958 }
3959
3960 // else: go on with next line
3961 }
3962}
3963
3964
3965
3966template <int dim, int spacedim>
3967void
3969 const char comment_start)
3970{
3971 char c;
3972 // loop over the following comment
3973 // lines
3974 while (in.get(c) && c == comment_start)
3975 // loop over the characters after
3976 // the comment starter
3977 while (in.get() != '\n')
3978 ;
3979
3980
3981 // put back first character of
3982 // first non-comment line
3983 if (in)
3984 in.putback(c);
3985
3986 // at last: skip additional empty lines, if present
3987 skip_empty_lines(in);
3988}
3989
3990
3991
3992template <int dim, int spacedim>
3993void
3995 const std::vector<CellData<dim>> & /*cells*/,
3996 const std::vector<Point<spacedim>> & /*vertices*/,
3997 std::ostream & /*out*/)
3998{
4000}
4001
4002
4003
4004template <>
4005void
4007 const std::vector<Point<2>> &vertices,
4008 std::ostream &out)
4009{
4010 double min_x = vertices[cells[0].vertices[0]][0],
4011 max_x = vertices[cells[0].vertices[0]][0],
4012 min_y = vertices[cells[0].vertices[0]][1],
4013 max_y = vertices[cells[0].vertices[0]][1];
4014
4015 for (unsigned int i = 0; i < cells.size(); ++i)
4016 {
4017 for (const auto vertex : cells[i].vertices)
4018 {
4019 const Point<2> &p = vertices[vertex];
4020
4021 if (p[0] < min_x)
4022 min_x = p[0];
4023 if (p[0] > max_x)
4024 max_x = p[0];
4025 if (p[1] < min_y)
4026 min_y = p[1];
4027 if (p[1] > max_y)
4028 max_y = p[1];
4029 }
4030
4031 out << "# cell " << i << std::endl;
4033 for (const auto vertex : cells[i].vertices)
4034 center += vertices[vertex];
4035 center /= 4;
4036
4037 out << "set label \"" << i << "\" at " << center[0] << ',' << center[1]
4038 << " center" << std::endl;
4039
4040 // first two line right direction
4041 for (unsigned int f = 0; f < 2; ++f)
4042 out << "set arrow from " << vertices[cells[i].vertices[f]][0] << ','
4043 << vertices[cells[i].vertices[f]][1] << " to "
4044 << vertices[cells[i].vertices[(f + 1) % 4]][0] << ','
4045 << vertices[cells[i].vertices[(f + 1) % 4]][1] << std::endl;
4046 // other two lines reverse direction
4047 for (unsigned int f = 2; f < 4; ++f)
4048 out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]][0]
4049 << ',' << vertices[cells[i].vertices[(f + 1) % 4]][1] << " to "
4050 << vertices[cells[i].vertices[f]][0] << ','
4051 << vertices[cells[i].vertices[f]][1] << std::endl;
4052 out << std::endl;
4053 }
4054
4055
4056 out << std::endl
4057 << "set nokey" << std::endl
4058 << "pl [" << min_x << ':' << max_x << "][" << min_y << ':' << max_y
4059 << "] " << min_y << std::endl
4060 << "pause -1" << std::endl;
4061}
4062
4063
4064
4065template <>
4066void
4068 const std::vector<Point<3>> &vertices,
4069 std::ostream &out)
4070{
4071 for (const auto &cell : cells)
4072 {
4073 // line 0
4074 out << vertices[cell.vertices[0]] << std::endl
4075 << vertices[cell.vertices[1]] << std::endl
4076 << std::endl
4077 << std::endl;
4078 // line 1
4079 out << vertices[cell.vertices[1]] << std::endl
4080 << vertices[cell.vertices[2]] << std::endl
4081 << std::endl
4082 << std::endl;
4083 // line 2
4084 out << vertices[cell.vertices[3]] << std::endl
4085 << vertices[cell.vertices[2]] << std::endl
4086 << std::endl
4087 << std::endl;
4088 // line 3
4089 out << vertices[cell.vertices[0]] << std::endl
4090 << vertices[cell.vertices[3]] << std::endl
4091 << std::endl
4092 << std::endl;
4093 // line 4
4094 out << vertices[cell.vertices[4]] << std::endl
4095 << vertices[cell.vertices[5]] << std::endl
4096 << std::endl
4097 << std::endl;
4098 // line 5
4099 out << vertices[cell.vertices[5]] << std::endl
4100 << vertices[cell.vertices[6]] << std::endl
4101 << std::endl
4102 << std::endl;
4103 // line 6
4104 out << vertices[cell.vertices[7]] << std::endl
4105 << vertices[cell.vertices[6]] << std::endl
4106 << std::endl
4107 << std::endl;
4108 // line 7
4109 out << vertices[cell.vertices[4]] << std::endl
4110 << vertices[cell.vertices[7]] << std::endl
4111 << std::endl
4112 << std::endl;
4113 // line 8
4114 out << vertices[cell.vertices[0]] << std::endl
4115 << vertices[cell.vertices[4]] << std::endl
4116 << std::endl
4117 << std::endl;
4118 // line 9
4119 out << vertices[cell.vertices[1]] << std::endl
4120 << vertices[cell.vertices[5]] << std::endl
4121 << std::endl
4122 << std::endl;
4123 // line 10
4124 out << vertices[cell.vertices[2]] << std::endl
4125 << vertices[cell.vertices[6]] << std::endl
4126 << std::endl
4127 << std::endl;
4128 // line 11
4129 out << vertices[cell.vertices[3]] << std::endl
4130 << vertices[cell.vertices[7]] << std::endl
4131 << std::endl
4132 << std::endl;
4133 }
4134}
4135
4136
4137
4138template <int dim, int spacedim>
4139void
4140GridIn<dim, spacedim>::read(const std::string &filename, Format format)
4141{
4142 // Search file class for meshes
4143 PathSearch search("MESH");
4144 std::string name;
4145 // Open the file and remember its name
4146 if (format == Default)
4147 name = search.find(filename);
4148 else
4149 name = search.find(filename, default_suffix(format));
4150
4151
4152 if (format == Default)
4153 {
4154 const std::string::size_type slashpos = name.find_last_of('/');
4155 const std::string::size_type dotpos = name.find_last_of('.');
4156 if (dotpos < name.size() &&
4157 (dotpos > slashpos || slashpos == std::string::npos))
4158 {
4159 std::string ext = name.substr(dotpos + 1);
4160 format = parse_format(ext);
4161 }
4162 }
4163
4164 if (format == assimp)
4165 {
4166 read_assimp(name);
4167 }
4168 else if (format == exodusii)
4169 {
4170 read_exodusii(name);
4171 }
4172 else
4173 {
4174 std::ifstream in(name);
4175 read(in, format);
4176 }
4177}
4178
4179
4180template <int dim, int spacedim>
4181void
4182GridIn<dim, spacedim>::read(std::istream &in, Format format)
4183{
4184 if (format == Default)
4185 format = default_format;
4186
4187 switch (format)
4188 {
4189 case dbmesh:
4190 read_dbmesh(in);
4191 return;
4192
4193 case msh:
4194 read_msh(in);
4195 return;
4196
4197 case vtk:
4198 read_vtk(in);
4199 return;
4200
4201 case vtu:
4202 read_vtu(in);
4203 return;
4204
4205 case unv:
4206 read_unv(in);
4207 return;
4208
4209 case ucd:
4210 read_ucd(in);
4211 return;
4212
4213 case abaqus:
4214 read_abaqus(in);
4215 return;
4216
4217 case xda:
4218 read_xda(in);
4219 return;
4220
4221 case tecplot:
4222 read_tecplot(in);
4223 return;
4224
4225 case assimp:
4226 Assert(false,
4227 ExcMessage("There is no read_assimp(istream &) function. "
4228 "Use the read_assimp(string &filename, ...) "
4229 "functions, instead."));
4230 return;
4231
4232 case exodusii:
4233 Assert(false,
4234 ExcMessage("There is no read_exodusii(istream &) function. "
4235 "Use the read_exodusii(string &filename, ...) "
4236 "function, instead."));
4237 return;
4238
4239 case Default:
4240 break;
4241 }
4243}
4244
4245
4246
4247template <int dim, int spacedim>
4248std::string
4250{
4251 switch (format)
4252 {
4253 case dbmesh:
4254 return ".dbmesh";
4255 case exodusii:
4256 return ".e";
4257 case msh:
4258 return ".msh";
4259 case vtk:
4260 return ".vtk";
4261 case vtu:
4262 return ".vtu";
4263 case unv:
4264 return ".unv";
4265 case ucd:
4266 return ".inp";
4267 case abaqus:
4268 return ".inp"; // Typical suffix for Abaqus mesh files conflicts with
4269 // UCD.
4270 case xda:
4271 return ".xda";
4272 case tecplot:
4273 return ".dat";
4274 default:
4276 return ".unknown_format";
4277 }
4278}
4279
4280
4281
4282template <int dim, int spacedim>
4284GridIn<dim, spacedim>::parse_format(const std::string &format_name)
4285{
4286 if (format_name == "dbmesh")
4287 return dbmesh;
4288
4289 if (format_name == "exodusii")
4290 return exodusii;
4291
4292 if (format_name == "msh")
4293 return msh;
4294
4295 if (format_name == "unv")
4296 return unv;
4297
4298 if (format_name == "vtk")
4299 return vtk;
4300
4301 if (format_name == "vtu")
4302 return vtu;
4303
4304 // This is also the typical extension of Abaqus input files.
4305 if (format_name == "inp")
4306 return ucd;
4307
4308 if (format_name == "ucd")
4309 return ucd;
4310
4311 if (format_name == "xda")
4312 return xda;
4313
4314 if (format_name == "tecplot")
4315 return tecplot;
4316
4317 if (format_name == "dat")
4318 return tecplot;
4319
4320 if (format_name == "plt")
4321 // Actually, this is the extension for the
4322 // tecplot binary format, which we do not
4323 // support right now. However, some people
4324 // tend to create tecplot ascii files with
4325 // the extension 'plt' instead of
4326 // 'dat'. Thus, include this extension
4327 // here. If it actually is a binary file,
4328 // the read_tecplot() function will fail
4329 // and throw an exception, anyway.
4330 return tecplot;
4331
4332 AssertThrow(false, ExcInvalidState());
4333 // return something weird
4334 return Format(Default);
4335}
4336
4337
4338
4339template <int dim, int spacedim>
4340std::string
4342{
4343 return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
4344}
4345
4346
4347
4348namespace
4349{
4350 template <int dim, int spacedim>
4351 Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
4352 : tolerance(5e-16) // Used to offset Cubit tolerance error when outputting
4353 // value close to zero
4354 {
4355 AssertThrow(spacedim == 2 || spacedim == 3, ExcNotImplemented());
4356 }
4357
4358
4359
4360 // Convert from a string to some other data type
4361 // Reference: http://www.codeguru.com/forum/showthread.php?t=231054
4362 template <class T>
4363 bool
4364 from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &))
4365 {
4366 std::istringstream iss(s);
4367 return !(iss >> f >> t).fail();
4368 }
4369
4370
4371
4372 // Extract an integer from a string
4373 int
4374 extract_int(const std::string &s)
4375 {
4376 std::string tmp;
4377 for (const char c : s)
4378 {
4379 if (isdigit(c) != 0)
4380 {
4381 tmp += c;
4382 }
4383 }
4384
4385 int number = 0;
4386 from_string(number, tmp, std::dec);
4387 return number;
4388 }
4389
4390
4391
4392 template <int dim, int spacedim>
4393 void
4394 Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
4395 {
4396 // References:
4397 // http://www.egr.msu.edu/software/abaqus/Documentation/docs/v6.7/books/usb/default.htm?startat=pt01ch02.html
4398 // http://www.cprogramming.com/tutorial/string.html
4399
4400 AssertThrow(input_stream.fail() == false, ExcIO());
4401 std::string line;
4402
4403 while (std::getline(input_stream, line))
4404 {
4405 cont:
4406 std::transform(line.begin(), line.end(), line.begin(), ::toupper);
4407
4408 if (line.compare("*HEADING") == 0 || line.compare(0, 2, "**") == 0 ||
4409 line.compare(0, 5, "*PART") == 0)
4410 {
4411 // Skip header and comments
4412 while (std::getline(input_stream, line))
4413 {
4414 if (line[0] == '*')
4415 goto cont; // My eyes, they burn!
4416 }
4417 }
4418 else if (line.compare(0, 5, "*NODE") == 0)
4419 {
4420 // Extract list of vertices
4421 // Header line might be:
4422 // *NODE, NSET=ALLNODES
4423 // *NODE
4424
4425 // Contains lines in the form:
4426 // Index, x, y, z
4427 while (std::getline(input_stream, line))
4428 {
4429 if (line[0] == '*')
4430 goto cont;
4431
4432 std::vector<double> node(spacedim + 1);
4433
4434 std::istringstream iss(line);
4435 char comma;
4436 for (unsigned int i = 0; i < spacedim + 1; ++i)
4437 iss >> node[i] >> comma;
4438
4439 node_list.push_back(node);
4440 }
4441 }
4442 else if (line.compare(0, 8, "*ELEMENT") == 0)
4443 {
4444 // Element construction.
4445 // There are different header formats, the details
4446 // of which we're not particularly interested in except
4447 // whether they represent quads or hexahedrals.
4448 // *ELEMENT, TYPE=S4R, ELSET=EB<material id>
4449 // *ELEMENT, TYPE=C3d8R, ELSET=EB<material id>
4450 // *ELEMENT, TYPE=C3d8
4451 // Elements itself (n=4 or n=8):
4452 // Index, i[0], ..., i[n]
4453
4454 int material = 0;
4455 // Scan for material id
4456 {
4457 const std::string before_material = "ELSET=EB";
4458 const std::size_t idx = line.find(before_material);
4459 if (idx != std::string::npos)
4460 {
4461 from_string(material,
4462 line.substr(idx + before_material.size()),
4463 std::dec);
4464 }
4465 }
4466
4467 // Read ELEMENT definition
4468 while (std::getline(input_stream, line))
4469 {
4470 if (line[0] == '*')
4471 goto cont;
4472
4473 std::istringstream iss(line);
4474 char comma;
4475
4476 // We will store the material id in the zeroth entry of the
4477 // vector and the rest of the elements represent the global
4478 // node numbers
4479 const unsigned int n_data_per_cell =
4481 std::vector<double> cell(n_data_per_cell);
4482 for (unsigned int i = 0; i < n_data_per_cell; ++i)
4483 iss >> cell[i] >> comma;
4484
4485 // Overwrite cell index from file by material
4486 cell[0] = static_cast<double>(material);
4487 cell_list.push_back(cell);
4488 }
4489 }
4490 else if (line.compare(0, 8, "*SURFACE") == 0)
4491 {
4492 // Extract the definitions of boundary surfaces
4493 // Old format from Cubit:
4494 // *SURFACE, NAME=SS<boundary indicator>
4495 // <element index>, S<face number>
4496 // Abaqus default format:
4497 // *SURFACE, TYPE=ELEMENT, NAME=SURF-<indicator>
4498
4499 // Get name of the surface and extract id from it;
4500 // this will be the boundary indicator
4501 const std::string name_key = "NAME=";
4502 const std::size_t name_idx_start =
4503 line.find(name_key) + name_key.size();
4504 std::size_t name_idx_end = line.find(',', name_idx_start);
4505 if (name_idx_end == std::string::npos)
4506 {
4507 name_idx_end = line.size();
4508 }
4509 const int b_indicator = extract_int(
4510 line.substr(name_idx_start, name_idx_end - name_idx_start));
4511
4512 // Read SURFACE definition
4513 // Note that the orientation of the faces is embedded within the
4514 // definition of each "set" of faces that comprise the surface
4515 // These are either marked by an "S" or "E" in 3d or 2d
4516 // respectively.
4517 while (std::getline(input_stream, line))
4518 {
4519 if (line[0] == '*')
4520 goto cont;
4521
4522 // Change all characters to upper case
4523 std::transform(line.begin(),
4524 line.end(),
4525 line.begin(),
4526 ::toupper);
4527
4528 // Surface can be created from ELSET, or directly from cells
4529 // If elsets_list contains a key with specific name - refers
4530 // to that ELSET, otherwise refers to cell
4531 std::istringstream iss(line);
4532 int el_idx;
4533 int face_number;
4534 char temp;
4535
4536 // Get relevant faces, taking into account the element
4537 // orientation
4538 std::vector<double> quad_node_list;
4539 const std::string elset_name = line.substr(0, line.find(','));
4540 if (elsets_list.count(elset_name) != 0)
4541 {
4542 // Surface refers to ELSET
4543 std::string stmp;
4544 iss >> stmp >> temp >> face_number;
4545
4546 const std::vector<int> cells = elsets_list[elset_name];
4547 for (const int cell : cells)
4548 {
4549 el_idx = cell;
4550 quad_node_list =
4551 get_global_node_numbers(el_idx, face_number);
4552 quad_node_list.insert(quad_node_list.begin(),
4553 b_indicator);
4554
4555 face_list.push_back(quad_node_list);
4556 }
4557 }
4558 else
4559 {
4560 // Surface refers directly to elements
4561 char comma;
4562 iss >> el_idx >> comma >> temp >> face_number;
4563 quad_node_list =
4564 get_global_node_numbers(el_idx, face_number);
4565 quad_node_list.insert(quad_node_list.begin(), b_indicator);
4566
4567 face_list.push_back(quad_node_list);
4568 }
4569 }
4570 }
4571 else if (line.compare(0, 6, "*ELSET") == 0)
4572 {
4573 // Get ELSET name.
4574 // Materials are attached to elsets with specific name
4575 std::string elset_name;
4576 {
4577 const std::string elset_key = "*ELSET, ELSET=";
4578 const std::size_t idx = line.find(elset_key);
4579 if (idx != std::string::npos)
4580 {
4581 const std::string comma = ",";
4582 const std::size_t first_comma = line.find(comma);
4583 const std::size_t second_comma =
4584 line.find(comma, first_comma + 1);
4585 const std::size_t elset_name_start =
4586 line.find(elset_key) + elset_key.size();
4587 elset_name = line.substr(elset_name_start,
4588 second_comma - elset_name_start);
4589 }
4590 }
4591
4592 // There are two possibilities of storing cells numbers in ELSET:
4593 // 1. If the header contains the 'GENERATE' keyword, then the next
4594 // line describes range of cells as:
4595 // cell_id_start, cell_id_end, cell_step
4596 // 2. If the header does not contain the 'GENERATE' keyword, then
4597 // the next lines contain cells numbers
4598 std::vector<int> elements;
4599 const std::size_t generate_idx = line.find("GENERATE");
4600 if (generate_idx != std::string::npos)
4601 {
4602 // Option (1)
4603 std::getline(input_stream, line);
4604 std::istringstream iss(line);
4605 char comma;
4606 int elid_start;
4607 int elid_end;
4608 int elis_step = 1; // Default if case stride not provided
4609
4610 // Some files don't have the stride size
4611 // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp
4612 // to
4613 // ./grids/abaqus/2d/2d_test_abaqus.inp
4614 iss >> elid_start >> comma >> elid_end;
4615 AssertThrow(comma == ',',
4616 ExcMessage(
4617 std::string(
4618 "While reading an ABAQUS file, the reader "
4619 "expected a comma but found a <") +
4620 comma + "> in the line <" + line + ">."));
4622 elid_start <= elid_end,
4623 ExcMessage(
4624 std::string(
4625 "While reading an ABAQUS file, the reader encountered "
4626 "a GENERATE statement in which the upper bound <") +
4627 Utilities::int_to_string(elid_end) +
4628 "> for the element numbers is not larger or equal "
4629 "than the lower bound <" +
4630 Utilities::int_to_string(elid_start) + ">."));
4631
4632 // https://stackoverflow.com/questions/8046357/how-do-i-check-if-a-stringstream-variable-is-empty-null
4633 if (iss.rdbuf()->in_avail() != 0)
4634 iss >> comma >> elis_step;
4635 AssertThrow(comma == ',',
4636 ExcMessage(
4637 std::string(
4638 "While reading an ABAQUS file, the reader "
4639 "expected a comma but found a <") +
4640 comma + "> in the line <" + line + ">."));
4641
4642 for (int i = elid_start; i <= elid_end; i += elis_step)
4643 elements.push_back(i);
4644 elsets_list[elset_name] = elements;
4645
4646 std::getline(input_stream, line);
4647 }
4648 else
4649 {
4650 // Option (2)
4651 while (std::getline(input_stream, line))
4652 {
4653 if (line[0] == '*')
4654 break;
4655
4656 std::istringstream iss(line);
4657 char comma;
4658 int elid;
4659 while (!iss.eof())
4660 {
4661 iss >> elid >> comma;
4663 comma == ',',
4664 ExcMessage(
4665 std::string(
4666 "While reading an ABAQUS file, the reader "
4667 "expected a comma but found a <") +
4668 comma + "> in the line <" + line + ">."));
4669
4670 elements.push_back(elid);
4671 }
4672 }
4673
4674 elsets_list[elset_name] = elements;
4675 }
4676
4677 goto cont;
4678 }
4679 else if (line.compare(0, 5, "*NSET") == 0)
4680 {
4681 // Skip nodesets; we have no use for them
4682 while (std::getline(input_stream, line))
4683 {
4684 if (line[0] == '*')
4685 goto cont;
4686 }
4687 }
4688 else if (line.compare(0, 14, "*SOLID SECTION") == 0)
4689 {
4690 // The ELSET name, which describes a section for particular
4691 // material
4692 const std::string elset_key = "ELSET=";
4693 const std::size_t elset_start =
4694 line.find("ELSET=") + elset_key.size();
4695 const std::size_t elset_end = line.find(',', elset_start + 1);
4696 const std::string elset_name =
4697 line.substr(elset_start, elset_end - elset_start);
4698
4699 // Solid material definition.
4700 // We assume that material id is taken from material name,
4701 // eg. "Material-1" -> ID=1
4702 const std::string material_key = "MATERIAL=";
4703 const std::size_t last_equal =
4704 line.find("MATERIAL=") + material_key.size();
4705 const std::size_t material_id_start = line.find('-', last_equal);
4706 int material_id = 0;
4707 from_string(material_id,
4708 line.substr(material_id_start + 1),
4709 std::dec);
4710
4711 // Assign material id to cells
4712 const std::vector<int> &elset_cells = elsets_list[elset_name];
4713 for (const int elset_cell : elset_cells)
4714 {
4715 const int cell_id = elset_cell - 1;
4716 cell_list[cell_id][0] = material_id;
4717 }
4718 }
4719 // Note: All other lines / entries are ignored
4720 }
4721 }
4722
4723 template <int dim, int spacedim>
4724 std::vector<double>
4725 Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
4726 const int face_cell_no,
4727 const int face_cell_face_no) const
4728 {
4729 std::vector<double> quad_node_list(GeometryInfo<dim>::vertices_per_face);
4730
4731 // These orderings were reverse engineered by hand and may
4732 // conceivably be erroneous.
4733 // TODO: Currently one test (2d unstructured mesh) in the test
4734 // suite fails, presumably because of an ordering issue.
4735 if (dim == 2)
4736 {
4737 if (face_cell_face_no == 1)
4738 {
4739 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4740 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4741 }
4742 else if (face_cell_face_no == 2)
4743 {
4744 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4745 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4746 }
4747 else if (face_cell_face_no == 3)
4748 {
4749 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4750 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4751 }
4752 else if (face_cell_face_no == 4)
4753 {
4754 quad_node_list[0] = cell_list[face_cell_no - 1][4];
4755 quad_node_list[1] = cell_list[face_cell_no - 1][1];
4756 }
4757 else
4758 {
4759 AssertThrow(face_cell_face_no <= 4,
4760 ExcMessage("Invalid face number in 2d"));
4761 }
4762 }
4763 else if (dim == 3)
4764 {
4765 if (face_cell_face_no == 1)
4766 {
4767 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4768 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4769 quad_node_list[2] = cell_list[face_cell_no - 1][3];
4770 quad_node_list[3] = cell_list[face_cell_no - 1][2];
4771 }
4772 else if (face_cell_face_no == 2)
4773 {
4774 quad_node_list[0] = cell_list[face_cell_no - 1][5];
4775 quad_node_list[1] = cell_list[face_cell_no - 1][8];
4776 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4777 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4778 }
4779 else if (face_cell_face_no == 3)
4780 {
4781 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4782 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4783 quad_node_list[2] = cell_list[face_cell_no - 1][6];
4784 quad_node_list[3] = cell_list[face_cell_no - 1][5];
4785 }
4786 else if (face_cell_face_no == 4)
4787 {
4788 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4789 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4790 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4791 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4792 }
4793 else if (face_cell_face_no == 5)
4794 {
4795 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4796 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4797 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4798 quad_node_list[3] = cell_list[face_cell_no - 1][7];
4799 }
4800 else if (face_cell_face_no == 6)
4801 {
4802 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4803 quad_node_list[1] = cell_list[face_cell_no - 1][5];
4804 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4805 quad_node_list[3] = cell_list[face_cell_no - 1][4];
4806 }
4807 else
4808 {
4809 AssertThrow(face_cell_no <= 6,
4810 ExcMessage("Invalid face number in 3d"));
4811 }
4812 }
4813 else
4814 {
4815 AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
4816 }
4817
4818 return quad_node_list;
4819 }
4820
4821 template <int dim, int spacedim>
4822 void
4823 Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output) const
4824 {
4825 // References:
4826 // http://www.dealii.org/developer/doxygen/deal.II/structGeometryInfo.html
4827 // http://people.scs.fsu.edu/~burkardt/data/ucd/ucd.html
4828
4829 AssertThrow(output.fail() == false, ExcIO());
4830
4831 // save old formatting options
4832 const boost::io::ios_base_all_saver formatting_saver(output);
4833
4834 // Write out title - Note: No other commented text can be inserted below
4835 // the title in a UCD file
4836 output << "# Abaqus to UCD mesh conversion" << std::endl;
4837 output << "# Mesh type: AVS UCD" << std::endl;
4838
4839 // ========================================================
4840 // ASCII UCD File Format
4841 // The input file cannot contain blank lines or lines with leading blanks.
4842 // Comments, if present, must precede all data in the file.
4843 // Comments within the data will cause read errors.
4844 // The general order of the data is as follows:
4845 // 1. Numbers defining the overall structure, including the number of
4846 // nodes,
4847 // the number of cells, and the length of the vector of data associated
4848 // with the nodes, cells, and the model.
4849 // e.g. 1:
4850 // <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
4851 // e.g. 2:
4852 // n_elements = n_hex_cells + n_bc_quads + n_quad_cells +
4853 // n_bc_edges outfile.write(str(n_nodes) + " " + str(n_elements) +
4854 // " 0 0 0\n")
4855 // 2. For each node, its node id and the coordinates of that node in
4856 // space.
4857 // Node-ids must be integers, but any number including non sequential
4858 // numbers can be used. Mid-edge nodes are treated like any other node.
4859 // 3. For each cell: its cell-id, material, cell type (hexahedral,
4860 // pyramid,
4861 // etc.), and the list of node-ids that correspond to each of the
4862 // cell's vertices. The below table specifies the different cell types
4863 // and the keyword used to represent them in the file.
4864
4865 // Write out header
4866 output << node_list.size() << "\t" << (cell_list.size() + face_list.size())
4867 << "\t0\t0\t0" << std::endl;
4868
4869 output.width(16);
4870 output.precision(8);
4871
4872 // Write out node numbers
4873 // Loop over all nodes
4874 for (const auto &node : node_list)
4875 {
4876 // Node number
4877 output << node[0] << "\t";
4878
4879 // Node coordinates
4880 output.setf(std::ios::scientific, std::ios::floatfield);
4881 for (unsigned int jj = 1; jj < spacedim + 1; ++jj)
4882 {
4883 // invoke tolerance -> set points close to zero equal to zero
4884 if (std::abs(node[jj]) > tolerance)
4885 output << static_cast<double>(node[jj]) << "\t";
4886 else
4887 output << 0.0 << "\t";
4888 }
4889 if (spacedim == 2)
4890 output << 0.0 << "\t";
4891
4892 output << std::endl;
4893 output.unsetf(std::ios::floatfield);
4894 }
4895
4896 // Write out cell node numbers
4897 for (unsigned int ii = 0; ii < cell_list.size(); ++ii)
4898 {
4899 output << ii + 1 << "\t" << cell_list[ii][0] << "\t"
4900 << (dim == 2 ? "quad" : "hex") << "\t";
4901 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
4902 ++jj)
4903 output << cell_list[ii][jj] << "\t";
4904
4905 output << std::endl;
4906 }
4907
4908 // Write out quad node numbers
4909 for (unsigned int ii = 0; ii < face_list.size(); ++ii)
4910 {
4911 output << ii + 1 << "\t" << face_list[ii][0] << "\t"
4912 << (dim == 2 ? "line" : "quad") << "\t";
4913 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
4914 ++jj)
4915 output << face_list[ii][jj] << "\t";
4916
4917 output << std::endl;
4918 }
4919 }
4920} // namespace
4921
4922
4923// explicit instantiations
4924#include "grid_in.inst"
4925
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:3938
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:3403
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:4249
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:3968
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:4284
void read_vtu(std::istream &in)
Definition grid_in.cc:588
void read_tecplot(std::istream &in)
Definition grid_in.cc:3394
ExodusIIData read_exodusii(const std::string &filename, const bool apply_all_indicators_to_manifolds=false)
Definition grid_in.cc:3773
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:4182
static void debug_output_grid(const std::vector< CellData< dim > > &cells, const std::vector< Point< spacedim > > &vertices, std::ostream &out)
Definition grid_in.cc:3994
static std::string get_format_names()
Definition grid_in.cc:4341
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:2929
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:4623
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