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