deal.II version GIT relicensing-1936-ge55124254b 2024-10-04 21:10:00+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
2005GridIn<dim, spacedim>::read_msh(std::istream &input_stream)
2006{
2007 Assert(tria != nullptr, ExcNoTriangulationSelected());
2008 AssertThrow(input_stream.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 // contain the content of the file stripped of the comments
2020 std::string stripped_file;
2021
2022 // Comments can be included by mesh generating software and must be deleted,
2023 // a string is filed with the content of the file stripped of the comments
2024 while (std::getline(input_stream, line))
2025 {
2026 if (line == "@f$Comments")
2027 {
2028 while (std::getline(input_stream, line))
2029 {
2030 if (line == "@f$EndComments")
2031 {
2032 break;
2033 }
2034 }
2035 continue;
2036 }
2037 stripped_file += line + '\n';
2038 }
2039
2040 // Restart reading the file normally since it has been stripped of comments
2041 std::istringstream in(stripped_file);
2042
2043 in >> line;
2044
2045 // first determine file format
2046 unsigned int gmsh_file_format = 0;
2047 if (line == "@f$NOD")
2048 gmsh_file_format = 10;
2049 else if (line == "@f$MeshFormat")
2050 gmsh_file_format = 20;
2051 else
2052 AssertThrow(false, ExcInvalidGMSHInput(line));
2053
2054 // if file format is 2.0 or greater then we also have to read the rest of
2055 // the header
2056 if (gmsh_file_format == 20)
2057 {
2058 double version;
2059 unsigned int file_type, data_size;
2060
2061 in >> version >> file_type >> data_size;
2062
2063 Assert((version >= 2.0) && (version <= 4.1), ExcNotImplemented());
2064 gmsh_file_format = static_cast<unsigned int>(version * 10);
2065
2066 Assert(file_type == 0, ExcNotImplemented());
2067 Assert(data_size == sizeof(double), ExcNotImplemented());
2068
2069 // read the end of the header and the first line of the nodes
2070 // description to synch ourselves with the format 1 handling above
2071 in >> line;
2072 AssertThrow(line == "@f$EndMeshFormat", ExcInvalidGMSHInput(line));
2073
2074 in >> line;
2075 // if the next block is of kind @f$PhysicalNames, ignore it
2076 if (line == "@f$PhysicalNames")
2077 {
2078 do
2079 {
2080 in >> line;
2081 }
2082 while (line != "@f$EndPhysicalNames");
2083 in >> line;
2084 }
2085
2086 // if the next block is of kind @f$Entities, parse it
2087 if (line == "@f$Entities")
2088 {
2089 unsigned long n_points, n_curves, n_surfaces, n_volumes;
2090
2091 in >> n_points >> n_curves >> n_surfaces >> n_volumes;
2092 for (unsigned int i = 0; i < n_points; ++i)
2093 {
2094 // parse point ids
2095 int tag;
2096 unsigned int n_physicals;
2097 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2098 box_max_z;
2099
2100 // we only care for 'tag' as key for tag_maps[0]
2101 if (gmsh_file_format > 40)
2102 {
2103 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2104 n_physicals;
2105 box_max_x = box_min_x;
2106 box_max_y = box_min_y;
2107 box_max_z = box_min_z;
2108 }
2109 else
2110 {
2111 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2112 box_max_x >> box_max_y >> box_max_z >> n_physicals;
2113 }
2114 // if there is a physical tag, we will use it as boundary id
2115 // below
2116 AssertThrow(n_physicals < 2,
2117 ExcMessage("More than one tag is not supported!"));
2118 // if there is no physical tag, use 0 as default
2119 int physical_tag = 0;
2120 for (unsigned int j = 0; j < n_physicals; ++j)
2121 in >> physical_tag;
2122 tag_maps[0][tag] = physical_tag;
2123 }
2124 for (unsigned int i = 0; i < n_curves; ++i)
2125 {
2126 // parse curve ids
2127 int tag;
2128 unsigned int n_physicals;
2129 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2130 box_max_z;
2131
2132 // we only care for 'tag' as key for tag_maps[1]
2133 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2134 box_max_y >> box_max_z >> n_physicals;
2135 // if there is a physical tag, we will use it as boundary id
2136 // below
2137 AssertThrow(n_physicals < 2,
2138 ExcMessage("More than one tag is not supported!"));
2139 // if there is no physical tag, use 0 as default
2140 int physical_tag = 0;
2141 for (unsigned int j = 0; j < n_physicals; ++j)
2142 in >> physical_tag;
2143 tag_maps[1][tag] = physical_tag;
2144 // we don't care about the points associated to a curve, but
2145 // have to parse them anyway because their format is
2146 // unstructured
2147 in >> n_points;
2148 for (unsigned int j = 0; j < n_points; ++j)
2149 in >> tag;
2150 }
2151
2152 for (unsigned int i = 0; i < n_surfaces; ++i)
2153 {
2154 // parse surface ids
2155 int tag;
2156 unsigned int n_physicals;
2157 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2158 box_max_z;
2159
2160 // we only care for 'tag' as key for tag_maps[2]
2161 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2162 box_max_y >> box_max_z >> n_physicals;
2163 // if there is a physical tag, we will use it as boundary id
2164 // below
2165 AssertThrow(n_physicals < 2,
2166 ExcMessage("More than one tag is not supported!"));
2167 // if there is no physical tag, use 0 as default
2168 int physical_tag = 0;
2169 for (unsigned int j = 0; j < n_physicals; ++j)
2170 in >> physical_tag;
2171 tag_maps[2][tag] = physical_tag;
2172 // we don't care about the curves associated to a surface, but
2173 // have to parse them anyway because their format is
2174 // unstructured
2175 in >> n_curves;
2176 for (unsigned int j = 0; j < n_curves; ++j)
2177 in >> tag;
2178 }
2179 for (unsigned int i = 0; i < n_volumes; ++i)
2180 {
2181 // parse volume ids
2182 int tag;
2183 unsigned int n_physicals;
2184 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2185 box_max_z;
2186
2187 // we only care for 'tag' as key for tag_maps[3]
2188 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2189 box_max_y >> box_max_z >> n_physicals;
2190 // if there is a physical tag, we will use it as boundary id
2191 // below
2192 AssertThrow(n_physicals < 2,
2193 ExcMessage("More than one tag is not supported!"));
2194 // if there is no physical tag, use 0 as default
2195 int physical_tag = 0;
2196 for (unsigned int j = 0; j < n_physicals; ++j)
2197 in >> physical_tag;
2198 tag_maps[3][tag] = physical_tag;
2199 // we don't care about the surfaces associated to a volume, but
2200 // have to parse them anyway because their format is
2201 // unstructured
2202 in >> n_surfaces;
2203 for (unsigned int j = 0; j < n_surfaces; ++j)
2204 in >> tag;
2205 }
2206 in >> line;
2207 AssertThrow(line == "@f$EndEntities", ExcInvalidGMSHInput(line));
2208 in >> line;
2209 }
2210
2211 // if the next block is of kind @f$PartitionedEntities, ignore it
2212 if (line == "@f$PartitionedEntities")
2213 {
2214 do
2215 {
2216 in >> line;
2217 }
2218 while (line != "@f$EndPartitionedEntities");
2219 in >> line;
2220 }
2221
2222 // but the next thing should,
2223 // in any case, be the list of
2224 // nodes:
2225 AssertThrow(line == "@f$Nodes", ExcInvalidGMSHInput(line));
2226 }
2227
2228 // now read the nodes list
2229 int n_entity_blocks = 1;
2230 if (gmsh_file_format > 40)
2231 {
2232 int min_node_tag;
2233 int max_node_tag;
2234 in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag;
2235 }
2236 else if (gmsh_file_format == 40)
2237 {
2238 in >> n_entity_blocks >> n_vertices;
2239 }
2240 else
2241 in >> n_vertices;
2242 std::vector<Point<spacedim>> vertices(n_vertices);
2243 // set up mapping between numbering
2244 // in msh-file (nod) and in the
2245 // vertices vector
2246 std::map<int, int> vertex_indices;
2247
2248 {
2249 unsigned int global_vertex = 0;
2250 for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2251 {
2252 int parametric;
2253 unsigned long numNodes;
2254
2255 if (gmsh_file_format < 40)
2256 {
2257 numNodes = n_vertices;
2258 parametric = 0;
2259 }
2260 else
2261 {
2262 // for gmsh_file_format 4.1 the order of tag and dim is reversed,
2263 // but we are ignoring both anyway.
2264 int tagEntity, dimEntity;
2265 in >> tagEntity >> dimEntity >> parametric >> numNodes;
2266 }
2267
2268 std::vector<int> vertex_numbers;
2269 int vertex_number;
2270 if (gmsh_file_format > 40)
2271 for (unsigned long vertex_per_entity = 0;
2272 vertex_per_entity < numNodes;
2273 ++vertex_per_entity)
2274 {
2275 in >> vertex_number;
2276 vertex_numbers.push_back(vertex_number);
2277 }
2278
2279 for (unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes;
2280 ++vertex_per_entity, ++global_vertex)
2281 {
2282 int vertex_number;
2283 double x[3];
2284
2285 // read vertex
2286 if (gmsh_file_format > 40)
2287 {
2288 vertex_number = vertex_numbers[vertex_per_entity];
2289 in >> x[0] >> x[1] >> x[2];
2290 }
2291 else
2292 in >> vertex_number >> x[0] >> x[1] >> x[2];
2293
2294 for (unsigned int d = 0; d < spacedim; ++d)
2295 vertices[global_vertex][d] = x[d];
2296 // store mapping
2297 vertex_indices[vertex_number] = global_vertex;
2298
2299 // ignore parametric coordinates
2300 if (parametric != 0)
2301 {
2302 double u = 0.;
2303 double v = 0.;
2304 in >> u >> v;
2305 (void)u;
2306 (void)v;
2307 }
2308 }
2309 }
2310 AssertDimension(global_vertex, n_vertices);
2311 }
2312
2313 // Assert we reached the end of the block
2314 in >> line;
2315 static const std::string end_nodes_marker[] = {"@f$ENDNOD", "@f$EndNodes"};
2316 AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1],
2317 ExcInvalidGMSHInput(line));
2318
2319 // Now read in next bit
2320 in >> line;
2321 static const std::string begin_elements_marker[] = {"@f$ELM", "@f$Elements"};
2322 AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2323 ExcInvalidGMSHInput(line));
2324
2325 // now read the cell list
2326 if (gmsh_file_format > 40)
2327 {
2328 int min_node_tag;
2329 int max_node_tag;
2330 in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag;
2331 }
2332 else if (gmsh_file_format == 40)
2333 {
2334 in >> n_entity_blocks >> n_cells;
2335 }
2336 else
2337 {
2338 n_entity_blocks = 1;
2339 in >> n_cells;
2340 }
2341
2342 // set up array of cells and subcells (faces). In 1d, there is currently no
2343 // standard way in deal.II to pass boundary indicators attached to
2344 // individual vertices, so do this by hand via the boundary_ids_1d array
2345 std::vector<CellData<dim>> cells;
2346 SubCellData subcelldata;
2347 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2348
2349 // Track the number of times each vertex is used in 1D. This determines
2350 // whether or not we can assign a boundary id to a vertex. This is necessary
2351 // because sometimes gmsh saves internal vertices in the @f$ELEM list in codim
2352 // 1 or codim 2.
2353 std::map<unsigned int, unsigned int> vertex_counts;
2354
2355 {
2356 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
2357 {0, 1, 5, 4, 2, 3, 7, 6}};
2358 unsigned int global_cell = 0;
2359 for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2360 {
2361 unsigned int material_id;
2362 unsigned long numElements;
2363 int cell_type;
2364
2365 if (gmsh_file_format < 40)
2366 {
2367 material_id = 0;
2368 cell_type = 0;
2369 numElements = n_cells;
2370 }
2371 else if (gmsh_file_format == 40)
2372 {
2373 int tagEntity, dimEntity;
2374 in >> tagEntity >> dimEntity >> cell_type >> numElements;
2375 material_id = tag_maps[dimEntity][tagEntity];
2376 }
2377 else
2378 {
2379 // for gmsh_file_format 4.1 the order of tag and dim is reversed,
2380 int tagEntity, dimEntity;
2381 in >> dimEntity >> tagEntity >> cell_type >> numElements;
2382 material_id = tag_maps[dimEntity][tagEntity];
2383 }
2384
2385 for (unsigned int cell_per_entity = 0; cell_per_entity < numElements;
2386 ++cell_per_entity, ++global_cell)
2387 {
2388 // note that since in the input
2389 // file we found the number of
2390 // cells at the top, there
2391 // should still be input here,
2392 // so check this:
2393 AssertThrow(in.fail() == false, ExcIO());
2394
2395 unsigned int nod_num;
2396
2397 /*
2398 For file format version 1, the format of each cell is as
2399 follows: elm-number elm-type reg-phys reg-elem number-of-nodes
2400 node-number-list
2401
2402 However, for version 2, the format reads like this:
2403 elm-number elm-type number-of-tags < tag > ...
2404 node-number-list
2405
2406 For version 4, we have:
2407 tag(int) numVert(int) ...
2408
2409 In the following, we will ignore the element number (we simply
2410 enumerate them in the order in which we read them, and we will
2411 take reg-phys (version 1) or the first tag (version 2, if any
2412 tag is given at all) as material id. For version 4, we already
2413 read the material and the cell type in above.
2414 */
2415
2416 unsigned int elm_number = 0;
2417 if (gmsh_file_format < 40)
2418 {
2419 in >> elm_number // ELM-NUMBER
2420 >> cell_type; // ELM-TYPE
2421 }
2422
2423 if (gmsh_file_format < 20)
2424 {
2425 in >> material_id // REG-PHYS
2426 >> dummy // reg_elm
2427 >> nod_num;
2428 }
2429 else if (gmsh_file_format < 40)
2430 {
2431 // read the tags; ignore all but the first one which we will
2432 // interpret as the material_id (for cells) or boundary_id
2433 // (for faces)
2434 unsigned int n_tags;
2435 in >> n_tags;
2436 if (n_tags > 0)
2437 in >> material_id;
2438 else
2439 material_id = 0;
2440
2441 for (unsigned int i = 1; i < n_tags; ++i)
2442 in >> dummy;
2443
2444 if (cell_type == 1) // line
2445 nod_num = 2;
2446 else if (cell_type == 2) // tri
2447 nod_num = 3;
2448 else if (cell_type == 3) // quad
2449 nod_num = 4;
2450 else if (cell_type == 4) // tet
2451 nod_num = 4;
2452 else if (cell_type == 5) // hex
2453 nod_num = 8;
2454 }
2455 else // file format version 4.0 and later
2456 {
2457 // ignore tag
2458 int tag;
2459 in >> tag;
2460
2461 if (cell_type == 1) // line
2462 nod_num = 2;
2463 else if (cell_type == 2) // tri
2464 nod_num = 3;
2465 else if (cell_type == 3) // quad
2466 nod_num = 4;
2467 else if (cell_type == 4) // tet
2468 nod_num = 4;
2469 else if (cell_type == 5) // hex
2470 nod_num = 8;
2471 }
2472
2473
2474 /* `ELM-TYPE'
2475 defines the geometrical type of the N-th element:
2476 `1'
2477 Line (2 nodes, 1 edge).
2478
2479 `2'
2480 Triangle (3 nodes, 3 edges).
2481
2482 `3'
2483 Quadrangle (4 nodes, 4 edges).
2484
2485 `4'
2486 Tetrahedron (4 nodes, 6 edges, 6 faces).
2487
2488 `5'
2489 Hexahedron (8 nodes, 12 edges, 6 faces).
2490
2491 `15'
2492 Point (1 node).
2493 */
2494
2495 if (((cell_type == 1) && (dim == 1)) || // a line in 1d
2496 ((cell_type == 2) && (dim == 2)) || // a triangle in 2d
2497 ((cell_type == 3) && (dim == 2)) || // a quadrilateral in 2d
2498 ((cell_type == 4) && (dim == 3)) || // a tet in 3d
2499 ((cell_type == 5) && (dim == 3))) // a hex in 3d
2500 // found a cell
2501 {
2502 unsigned int vertices_per_cell = 0;
2503 if (cell_type == 1) // line
2504 vertices_per_cell = 2;
2505 else if (cell_type == 2) // tri
2506 vertices_per_cell = 3;
2507 else if (cell_type == 3) // quad
2508 vertices_per_cell = 4;
2509 else if (cell_type == 4) // tet
2510 vertices_per_cell = 4;
2511 else if (cell_type == 5) // hex
2512 vertices_per_cell = 8;
2513
2514 AssertThrow(nod_num == vertices_per_cell,
2515 ExcMessage(
2516 "Number of nodes does not coincide with the "
2517 "number required for this object"));
2518
2519 // allocate and read indices
2520 cells.emplace_back();
2521 CellData<dim> &cell = cells.back();
2522 cell.vertices.resize(vertices_per_cell);
2523 for (unsigned int i = 0; i < vertices_per_cell; ++i)
2524 {
2525 // hypercube cells need to be reordered
2526 if (vertices_per_cell ==
2528 in >> cell.vertices[dim == 3 ?
2529 local_vertex_numbering[i] :
2531 else
2532 in >> cell.vertices[i];
2533 }
2534
2535 // to make sure that the cast won't fail
2536 Assert(material_id <=
2537 std::numeric_limits<types::material_id>::max(),
2539 material_id,
2540 0,
2541 std::numeric_limits<types::material_id>::max()));
2542 // we use only material_ids in the range from 0 to
2543 // numbers::invalid_material_id-1
2545
2546 cell.material_id = material_id;
2547
2548 // transform from gmsh to consecutive numbering
2549 for (unsigned int i = 0; i < vertices_per_cell; ++i)
2550 {
2551 AssertThrow(vertex_indices.find(cell.vertices[i]) !=
2552 vertex_indices.end(),
2553 ExcInvalidVertexIndexGmsh(cell_per_entity,
2554 elm_number,
2555 cell.vertices[i]));
2556
2557 const auto vertex = vertex_indices[cell.vertices[i]];
2558 if (dim == 1)
2559 vertex_counts[vertex] += 1u;
2560 cell.vertices[i] = vertex;
2561 }
2562 }
2563 else if ((cell_type == 1) &&
2564 ((dim == 2) || (dim == 3))) // a line in 2d or 3d
2565 // boundary info
2566 {
2567 subcelldata.boundary_lines.emplace_back();
2568 in >> subcelldata.boundary_lines.back().vertices[0] >>
2569 subcelldata.boundary_lines.back().vertices[1];
2570
2571 // to make sure that the cast won't fail
2572 Assert(material_id <=
2573 std::numeric_limits<types::boundary_id>::max(),
2575 material_id,
2576 0,
2577 std::numeric_limits<types::boundary_id>::max()));
2578 // we use only boundary_ids in the range from 0 to
2579 // numbers::internal_face_boundary_id-1
2580 AssertIndexRange(material_id,
2582
2583 subcelldata.boundary_lines.back().boundary_id =
2584 static_cast<types::boundary_id>(material_id);
2585
2586 // transform from ucd to
2587 // consecutive numbering
2588 for (unsigned int &vertex :
2589 subcelldata.boundary_lines.back().vertices)
2590 if (vertex_indices.find(vertex) != vertex_indices.end())
2591 // vertex with this index exists
2592 vertex = vertex_indices[vertex];
2593 else
2594 {
2595 // no such vertex index
2596 AssertThrow(false,
2597 ExcInvalidVertexIndex(cell_per_entity,
2598 vertex));
2600 }
2601 }
2602 else if ((cell_type == 2 || cell_type == 3) &&
2603 (dim == 3)) // a triangle or a quad in 3d
2604 // boundary info
2605 {
2606 unsigned int vertices_per_cell = 0;
2607 // check cell type
2608 if (cell_type == 2) // tri
2609 vertices_per_cell = 3;
2610 else if (cell_type == 3) // quad
2611 vertices_per_cell = 4;
2612
2613 subcelldata.boundary_quads.emplace_back();
2614
2615 // resize vertices
2616 subcelldata.boundary_quads.back().vertices.resize(
2617 vertices_per_cell);
2618 // for loop
2619 for (unsigned int i = 0; i < vertices_per_cell; ++i)
2620 in >> subcelldata.boundary_quads.back().vertices[i];
2621
2622 // to make sure that the cast won't fail
2623 Assert(material_id <=
2624 std::numeric_limits<types::boundary_id>::max(),
2626 material_id,
2627 0,
2628 std::numeric_limits<types::boundary_id>::max()));
2629 // we use only boundary_ids in the range from 0 to
2630 // numbers::internal_face_boundary_id-1
2631 AssertIndexRange(material_id,
2633
2634 subcelldata.boundary_quads.back().boundary_id =
2635 static_cast<types::boundary_id>(material_id);
2636
2637 // transform from gmsh to
2638 // consecutive numbering
2639 for (unsigned int &vertex :
2640 subcelldata.boundary_quads.back().vertices)
2641 if (vertex_indices.find(vertex) != vertex_indices.end())
2642 // vertex with this index exists
2643 vertex = vertex_indices[vertex];
2644 else
2645 {
2646 // no such vertex index
2647 Assert(false,
2648 ExcInvalidVertexIndex(cell_per_entity, vertex));
2650 }
2651 }
2652 else if (cell_type == 15)
2653 {
2654 // read the indices of nodes given
2655 unsigned int node_index = 0;
2656 if (gmsh_file_format < 20)
2657 {
2658 // For points (cell_type==15), we can only ever
2659 // list one node index.
2660 AssertThrow(nod_num == 1, ExcInternalError());
2661 in >> node_index;
2662 }
2663 else
2664 {
2665 in >> node_index;
2666 }
2667
2668 // we only care about boundary indicators assigned to
2669 // individual vertices in 1d (because otherwise the vertices
2670 // are not faces)
2671 if (dim == 1)
2672 boundary_ids_1d[vertex_indices[node_index]] = material_id;
2673 }
2674 else
2675 {
2676 AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type));
2677 }
2678 }
2679 }
2680 AssertDimension(global_cell, n_cells);
2681 }
2682 // Assert that we reached the end of the block
2683 in >> line;
2684 static const std::string end_elements_marker[] = {"@f$ENDELM", "@f$EndElements"};
2685 AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2686 ExcInvalidGMSHInput(line));
2687 AssertThrow(in.fail() == false, ExcIO());
2688
2689 // check that we actually read some cells.
2690 AssertThrow(cells.size() > 0,
2691 ExcGmshNoCellInformation(subcelldata.boundary_lines.size(),
2692 subcelldata.boundary_quads.size()));
2693
2694 // apply_grid_fixup_functions() may invalidate the vertex indices (since it
2695 // will delete duplicated or unused vertices). Get around this by storing
2696 // Points directly in that case so that the comparisons are valid.
2697 std::vector<std::pair<Point<spacedim>, types::boundary_id>> boundary_id_pairs;
2698 if (dim == 1)
2699 for (const auto &pair : vertex_counts)
2700 if (pair.second == 1u)
2701 boundary_id_pairs.emplace_back(vertices[pair.first],
2702 boundary_ids_1d[pair.first]);
2703
2704 apply_grid_fixup_functions(vertices, cells, subcelldata);
2705 tria->create_triangulation(vertices, cells, subcelldata);
2706
2707 // in 1d, we also have to attach boundary ids to vertices, which does not
2708 // currently work through the call above.
2709 if (dim == 1)
2710 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2711}
2712
2713
2714
2715#ifdef DEAL_II_GMSH_WITH_API
2716template <int dim, int spacedim>
2717void
2718GridIn<dim, spacedim>::read_msh(const std::string &fname)
2719{
2720 Assert(tria != nullptr, ExcNoTriangulationSelected());
2721 // gmsh -> deal.II types
2722 const std::map<int, std::uint8_t> gmsh_to_dealii_type = {
2723 {{15, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {7, 5}, {6, 6}, {5, 7}}};
2724
2725 // Vertex renumbering, by dealii type
2726 const std::array<std::vector<unsigned int>, 8> gmsh_to_dealii = {
2727 {{0},
2728 {{0, 1}},
2729 {{0, 1, 2}},
2730 {{0, 1, 3, 2}},
2731 {{0, 1, 2, 3}},
2732 {{0, 1, 3, 2, 4}},
2733 {{0, 1, 2, 3, 4, 5}},
2734 {{0, 1, 3, 2, 4, 5, 7, 6}}}};
2735
2736 std::vector<Point<spacedim>> vertices;
2737 std::vector<CellData<dim>> cells;
2738 SubCellData subcelldata;
2739 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2740
2741 // Track the number of times each vertex is used in 1D. This determines
2742 // whether or not we can assign a boundary id to a vertex. This is necessary
2743 // because sometimes gmsh saves internal vertices in the @f$ELEM list in codim
2744 // 1 or codim 2.
2745 std::map<unsigned int, unsigned int> vertex_counts;
2746
2747 gmsh::initialize();
2748 gmsh::option::setNumber("General.Verbosity", 0);
2749 gmsh::open(fname);
2750
2751 AssertThrow(gmsh::model::getDimension() == dim,
2752 ExcMessage("You are trying to read a gmsh file with dimension " +
2753 std::to_string(gmsh::model::getDimension()) +
2754 " into a grid of dimension " + std::to_string(dim)));
2755
2756 // Read all nodes, and store them in our vector of vertices. Before we do
2757 // that, make sure all tags are consecutive
2758 {
2759 gmsh::model::mesh::removeDuplicateNodes();
2760 gmsh::model::mesh::renumberNodes();
2761 std::vector<std::size_t> node_tags;
2762 std::vector<double> coord;
2763 std::vector<double> parametricCoord;
2764 gmsh::model::mesh::getNodes(
2765 node_tags, coord, parametricCoord, -1, -1, false, false);
2766 vertices.resize(node_tags.size());
2767 for (unsigned int i = 0; i < node_tags.size(); ++i)
2768 {
2769 // Check that renumbering worked!
2770 AssertDimension(node_tags[i], i + 1);
2771 for (unsigned int d = 0; d < spacedim; ++d)
2772 vertices[i][d] = coord[i * 3 + d];
2773# ifdef DEBUG
2774 // Make sure the embedded dimension is right
2775 for (unsigned int d = spacedim; d < 3; ++d)
2776 Assert(std::abs(coord[i * 3 + d]) < 1e-10,
2777 ExcMessage("The grid you are reading contains nodes that are "
2778 "nonzero in the coordinate with index " +
2779 std::to_string(d) +
2780 ", but you are trying to save "
2781 "it on a grid embedded in a " +
2782 std::to_string(spacedim) + " dimensional space."));
2783# endif
2784 }
2785 }
2786
2787 // Get all the elementary entities in the model, as a vector of (dimension,
2788 // tag) pairs:
2789 std::vector<std::pair<int, int>> entities;
2790 gmsh::model::getEntities(entities);
2791
2792 for (const auto &e : entities)
2793 {
2794 // Dimension and tag of the entity:
2795 const int &entity_dim = e.first;
2796 const int &entity_tag = e.second;
2797
2799 types::boundary_id boundary_id = 0;
2800
2801 // Get the physical tags, to deduce boundary, material, and manifold_id
2802 std::vector<int> physical_tags;
2803 gmsh::model::getPhysicalGroupsForEntity(entity_dim,
2804 entity_tag,
2805 physical_tags);
2806
2807 // Now fill manifold id and boundary or material id
2808 if (physical_tags.size())
2809 for (auto physical_tag : physical_tags)
2810 {
2811 std::string name;
2812 gmsh::model::getPhysicalName(entity_dim, physical_tag, name);
2813 if (!name.empty())
2814 {
2815 // Patterns::Tools::to_value throws an exception, if it can not
2816 // convert name to a map from string to int.
2817 try
2818 {
2819 std::map<std::string, int> id_names;
2820 Patterns::Tools::to_value(name, id_names);
2821 bool found_unrecognized_tag = false;
2822 bool found_boundary_id = false;
2823 // If the above did not throw, we keep going, and retrieve
2824 // all the information that we know how to translate.
2825 for (const auto &it : id_names)
2826 {
2827 const auto &name = it.first;
2828 const auto &id = it.second;
2829 if (entity_dim == dim && name == "MaterialID")
2830 {
2831 boundary_id = static_cast<types::boundary_id>(id);
2832 found_boundary_id = true;
2833 }
2834 else if (entity_dim < dim && name == "BoundaryID")
2835 {
2836 boundary_id = static_cast<types::boundary_id>(id);
2837 found_boundary_id = true;
2838 }
2839 else if (name == "ManifoldID")
2840 manifold_id = static_cast<types::manifold_id>(id);
2841 else
2842 // We did not recognize one of the keys. We'll fall
2843 // back to setting the boundary id to the physical tag
2844 // after reading all strings.
2845 found_unrecognized_tag = true;
2846 }
2847 // If we didn't find a BoundaryID:XX or MaterialID:XX, and
2848 // something was found but not recognized, then we set the
2849 // material id or boundary using the physical tag directly.
2850 if (found_unrecognized_tag && !found_boundary_id)
2851 boundary_id = physical_tag;
2852 }
2853 catch (...)
2854 {
2855 // When the above didn't work, we revert to the old
2856 // behaviour: the physical tag itself is interpreted either
2857 // as a material_id or a boundary_id, and no manifold id is
2858 // known
2859 boundary_id = physical_tag;
2860 }
2861 }
2862 }
2863
2864 // Get the mesh elements for the entity (dim, tag):
2865 std::vector<int> element_types;
2866 std::vector<std::vector<std::size_t>> element_ids, element_nodes;
2867 gmsh::model::mesh::getElements(
2868 element_types, element_ids, element_nodes, entity_dim, entity_tag);
2869
2870 for (unsigned int i = 0; i < element_types.size(); ++i)
2871 {
2872 const auto &type = gmsh_to_dealii_type.at(element_types[i]);
2873 const auto n_vertices = gmsh_to_dealii[type].size();
2874 const auto &elements = element_ids[i];
2875 const auto &nodes = element_nodes[i];
2876 for (unsigned int j = 0; j < elements.size(); ++j)
2877 {
2878 if (entity_dim == dim)
2879 {
2880 cells.emplace_back(n_vertices);
2881 auto &cell = cells.back();
2882 for (unsigned int v = 0; v < n_vertices; ++v)
2883 {
2884 cell.vertices[v] =
2885 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2886 if (dim == 1)
2887 vertex_counts[cell.vertices[v]] += 1u;
2888 }
2889 cell.manifold_id = manifold_id;
2890 cell.material_id = boundary_id;
2891 }
2892 else if (entity_dim == 2)
2893 {
2894 subcelldata.boundary_quads.emplace_back(n_vertices);
2895 auto &face = subcelldata.boundary_quads.back();
2896 for (unsigned int v = 0; v < n_vertices; ++v)
2897 face.vertices[v] =
2898 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2899
2900 face.manifold_id = manifold_id;
2901 face.boundary_id = boundary_id;
2902 }
2903 else if (entity_dim == 1)
2904 {
2905 subcelldata.boundary_lines.emplace_back(n_vertices);
2906 auto &line = subcelldata.boundary_lines.back();
2907 for (unsigned int v = 0; v < n_vertices; ++v)
2908 line.vertices[v] =
2909 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2910
2911 line.manifold_id = manifold_id;
2912 line.boundary_id = boundary_id;
2913 }
2914 else if (entity_dim == 0)
2915 {
2916 // This should only happen in one dimension.
2917 AssertDimension(dim, 1);
2918 for (unsigned int j = 0; j < elements.size(); ++j)
2919 boundary_ids_1d[nodes[j] - 1] = boundary_id;
2920 }
2921 }
2922 }
2923 }
2924
2925 // apply_grid_fixup_functions() may invalidate the vertex indices (since it
2926 // will delete duplicated or unused vertices). Get around this by storing
2927 // Points directly in that case so that the comparisons are valid.
2928 std::vector<std::pair<Point<spacedim>, types::boundary_id>> boundary_id_pairs;
2929 if (dim == 1)
2930 for (const auto &pair : vertex_counts)
2931 if (pair.second == 1u)
2932 boundary_id_pairs.emplace_back(vertices[pair.first],
2933 boundary_ids_1d[pair.first]);
2934
2935 apply_grid_fixup_functions(vertices, cells, subcelldata);
2936 tria->create_triangulation(vertices, cells, subcelldata);
2937
2938 // in 1d, we also have to attach boundary ids to vertices, which does not
2939 // currently work through the call above.
2940 if (dim == 1)
2941 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2942
2943 gmsh::clear();
2944 gmsh::finalize();
2945}
2946#endif
2947
2948
2949
2950template <int dim, int spacedim>
2951void
2953 std::string &header,
2954 std::vector<unsigned int> &tecplot2deal,
2955 unsigned int &n_vars,
2956 unsigned int &n_vertices,
2957 unsigned int &n_cells,
2958 std::vector<unsigned int> &IJK,
2959 bool &structured,
2960 bool &blocked)
2961{
2962 Assert(tecplot2deal.size() == dim, ExcInternalError());
2963 Assert(IJK.size() == dim, ExcInternalError());
2964 // initialize the output variables
2965 n_vars = 0;
2966 n_vertices = 0;
2967 n_cells = 0;
2968 switch (dim)
2969 {
2970 case 3:
2971 IJK[2] = 0;
2973 case 2:
2974 IJK[1] = 0;
2976 case 1:
2977 IJK[0] = 0;
2978 }
2979 structured = true;
2980 blocked = false;
2981
2982 // convert the string to upper case
2983 std::transform(header.begin(), header.end(), header.begin(), ::toupper);
2984
2985 // replace all tabs, commas, newlines by
2986 // whitespaces
2987 std::replace(header.begin(), header.end(), '\t', ' ');
2988 std::replace(header.begin(), header.end(), ',', ' ');
2989 std::replace(header.begin(), header.end(), '\n', ' ');
2990
2991 // now remove whitespace in front of and
2992 // after '='
2993 std::string::size_type pos = header.find('=');
2994
2995 while (pos != static_cast<std::string::size_type>(std::string::npos))
2996 if (header[pos + 1] == ' ')
2997 header.erase(pos + 1, 1);
2998 else if (header[pos - 1] == ' ')
2999 {
3000 header.erase(pos - 1, 1);
3001 --pos;
3002 }
3003 else
3004 pos = header.find('=', ++pos);
3005
3006 // split the string into individual entries
3007 std::vector<std::string> entries =
3008 Utilities::break_text_into_lines(header, 1, ' ');
3009
3010 // now go through the list and try to extract
3011 for (unsigned int i = 0; i < entries.size(); ++i)
3012 {
3013 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\""))
3014 {
3015 ++n_vars;
3016 // we assume, that the first variable
3017 // is x or no coordinate at all (not y or z)
3018 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\"X\""))
3019 {
3020 tecplot2deal[0] = 0;
3021 }
3022 ++i;
3023 while (entries[i][0] == '"')
3024 {
3025 if (entries[i] == "\"X\"")
3026 tecplot2deal[0] = n_vars;
3027 else if (entries[i] == "\"Y\"")
3028 {
3029 // we assume, that y contains
3030 // zero data in 1d, so do
3031 // nothing
3032 if (dim > 1)
3033 tecplot2deal[1] = n_vars;
3034 }
3035 else if (entries[i] == "\"Z\"")
3036 {
3037 // we assume, that z contains
3038 // zero data in 1d and 2d, so
3039 // do nothing
3040 if (dim > 2)
3041 tecplot2deal[2] = n_vars;
3042 }
3043 ++n_vars;
3044 ++i;
3045 }
3046 // set i back, so that the next
3047 // string is treated correctly
3048 --i;
3049
3051 n_vars >= dim,
3052 ExcMessage(
3053 "Tecplot file must contain at least one variable for each dimension"));
3054 for (unsigned int d = 1; d < dim; ++d)
3056 tecplot2deal[d] > 0,
3057 ExcMessage(
3058 "Tecplot file must contain at least one variable for each dimension."));
3059 }
3060 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE=ORDERED"))
3061 structured = true;
3062 else if (Utilities::match_at_string_start(entries[i],
3063 "ZONETYPE=FELINESEG") &&
3064 dim == 1)
3065 structured = false;
3066 else if (Utilities::match_at_string_start(entries[i],
3067 "ZONETYPE=FEQUADRILATERAL") &&
3068 dim == 2)
3069 structured = false;
3070 else if (Utilities::match_at_string_start(entries[i],
3071 "ZONETYPE=FEBRICK") &&
3072 dim == 3)
3073 structured = false;
3074 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE="))
3075 // unsupported ZONETYPE
3076 {
3077 AssertThrow(false,
3078 ExcMessage(
3079 "The tecplot file contains an unsupported ZONETYPE."));
3080 }
3081 else if (Utilities::match_at_string_start(entries[i],
3082 "DATAPACKING=POINT"))
3083 blocked = false;
3084 else if (Utilities::match_at_string_start(entries[i],
3085 "DATAPACKING=BLOCK"))
3086 blocked = true;
3087 else if (Utilities::match_at_string_start(entries[i], "F=POINT"))
3088 {
3089 structured = true;
3090 blocked = false;
3091 }
3092 else if (Utilities::match_at_string_start(entries[i], "F=BLOCK"))
3093 {
3094 structured = true;
3095 blocked = true;
3096 }
3097 else if (Utilities::match_at_string_start(entries[i], "F=FEPOINT"))
3098 {
3099 structured = false;
3100 blocked = false;
3101 }
3102 else if (Utilities::match_at_string_start(entries[i], "F=FEBLOCK"))
3103 {
3104 structured = false;
3105 blocked = true;
3106 }
3107 else if (Utilities::match_at_string_start(entries[i],
3108 "ET=QUADRILATERAL") &&
3109 dim == 2)
3110 structured = false;
3111 else if (Utilities::match_at_string_start(entries[i], "ET=BRICK") &&
3112 dim == 3)
3113 structured = false;
3114 else if (Utilities::match_at_string_start(entries[i], "ET="))
3115 // unsupported ElementType
3116 {
3118 false,
3119 ExcMessage(
3120 "The tecplot file contains an unsupported ElementType."));
3121 }
3122 else if (Utilities::match_at_string_start(entries[i], "I="))
3123 IJK[0] = Utilities::get_integer_at_position(entries[i], 2).first;
3124 else if (Utilities::match_at_string_start(entries[i], "J="))
3125 {
3126 IJK[1] = Utilities::get_integer_at_position(entries[i], 2).first;
3128 dim > 1 || IJK[1] == 1,
3129 ExcMessage(
3130 "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
3131 }
3132 else if (Utilities::match_at_string_start(entries[i], "K="))
3133 {
3134 IJK[2] = Utilities::get_integer_at_position(entries[i], 2).first;
3136 dim > 2 || IJK[2] == 1,
3137 ExcMessage(
3138 "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
3139 }
3140 else if (Utilities::match_at_string_start(entries[i], "N="))
3141 n_vertices = Utilities::get_integer_at_position(entries[i], 2).first;
3142 else if (Utilities::match_at_string_start(entries[i], "E="))
3143 n_cells = Utilities::get_integer_at_position(entries[i], 2).first;
3144 }
3145
3146 // now we have read all the fields we are
3147 // interested in. do some checks and
3148 // calculate the variables
3149 if (structured)
3150 {
3151 n_vertices = 1;
3152 n_cells = 1;
3153 for (unsigned int d = 0; d < dim; ++d)
3154 {
3156 IJK[d] > 0,
3157 ExcMessage(
3158 "Tecplot file does not contain a complete and consistent set of parameters"));
3159 n_vertices *= IJK[d];
3160 n_cells *= (IJK[d] - 1);
3161 }
3162 }
3163 else
3164 {
3166 n_vertices > 0,
3167 ExcMessage(
3168 "Tecplot file does not contain a complete and consistent set of parameters"));
3169 if (n_cells == 0)
3170 // this means an error, although
3171 // tecplot itself accepts entries like
3172 // 'J=20' instead of 'E=20'. therefore,
3173 // take the max of IJK
3174 n_cells = *std::max_element(IJK.begin(), IJK.end());
3176 n_cells > 0,
3177 ExcMessage(
3178 "Tecplot file does not contain a complete and consistent set of parameters"));
3179 }
3180}
3181
3182
3183
3184template <>
3185void
3187{
3188 const unsigned int dim = 2;
3189 const unsigned int spacedim = 2;
3190 Assert(tria != nullptr, ExcNoTriangulationSelected());
3191 AssertThrow(in.fail() == false, ExcIO());
3192
3193 // skip comments at start of file
3194 skip_comment_lines(in, '#');
3195
3196 // some strings for parsing the header
3197 std::string line, header;
3198
3199 // first, concatenate all header lines
3200 // create a searchstring with almost all
3201 // letters. exclude e and E from the letters
3202 // to search, as they might appear in
3203 // exponential notation
3204 std::string letters = "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
3205
3206 getline(in, line);
3207 while (line.find_first_of(letters) != std::string::npos)
3208 {
3209 header += " " + line;
3210 getline(in, line);
3211 }
3212
3213 // now create some variables holding
3214 // important information on the mesh, get
3215 // this information from the header string
3216 std::vector<unsigned int> tecplot2deal(dim);
3217 std::vector<unsigned int> IJK(dim);
3218 unsigned int n_vars, n_vertices, n_cells;
3219 bool structured, blocked;
3220
3221 parse_tecplot_header(header,
3222 tecplot2deal,
3223 n_vars,
3224 n_vertices,
3225 n_cells,
3226 IJK,
3227 structured,
3228 blocked);
3229
3230 // reserve space for vertices. note, that in
3231 // tecplot vertices are ordered beginning
3232 // with 1, whereas in deal all indices start
3233 // with 0. in order not to use -1 for all the
3234 // connectivity information, a 0th vertex
3235 // (unused) is inserted at the origin.
3236 std::vector<Point<spacedim>> vertices(n_vertices + 1);
3238 // reserve space for cells
3239 std::vector<CellData<dim>> cells(n_cells);
3240 SubCellData subcelldata;
3241
3242 if (blocked)
3243 {
3244 // blocked data format. first we get all
3245 // the values of the first variable for
3246 // all points, after that we get all
3247 // values for the second variable and so
3248 // on.
3249
3250 // dummy variable to read in all the info
3251 // we do not want to use
3252 double dummy;
3253 // which is the first index to read in
3254 // the loop (see below)
3255 unsigned int next_index = 0;
3256
3257 // note, that we have already read the
3258 // first line containing the first variable
3259 if (tecplot2deal[0] == 0)
3260 {
3261 // we need the information in this
3262 // line, so extract it
3263 std::vector<std::string> first_var =
3265 char *endptr;
3266 for (unsigned int i = 1; i < first_var.size() + 1; ++i)
3267 vertices[i][0] = std::strtod(first_var[i - 1].c_str(), &endptr);
3268
3269 // if there are many points, the data
3270 // for this var might continue in the
3271 // next line(s)
3272 for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
3273 in >> vertices[j][next_index];
3274 // now we got all values of the first
3275 // variable, so increase the counter
3276 next_index = 1;
3277 }
3278
3279 // main loop over all variables
3280 for (unsigned int i = 1; i < n_vars; ++i)
3281 {
3282 // if we read all the important
3283 // variables and do not want to
3284 // read further, because we are
3285 // using a structured grid, we can
3286 // stop here (and skip, for
3287 // example, a whole lot of solution
3288 // variables)
3289 if (next_index == dim && structured)
3290 break;
3291
3292 if ((next_index < dim) && (i == tecplot2deal[next_index]))
3293 {
3294 // we need this line, read it in
3295 for (unsigned int j = 1; j < n_vertices + 1; ++j)
3296 in >> vertices[j][next_index];
3297 ++next_index;
3298 }
3299 else
3300 {
3301 // we do not need this line, read
3302 // it in and discard it
3303 for (unsigned int j = 1; j < n_vertices + 1; ++j)
3304 in >> dummy;
3305 }
3306 }
3307 Assert(next_index == dim, ExcInternalError());
3308 }
3309 else
3310 {
3311 // the data is not blocked, so we get all
3312 // the variables for one point, then the
3313 // next and so on. create a vector to
3314 // hold these components
3315 std::vector<double> vars(n_vars);
3316
3317 // now fill the first vertex. note, that we
3318 // have already read the first line
3319 // containing the first vertex
3320 std::vector<std::string> first_vertex =
3322 char *endptr;
3323 for (unsigned int d = 0; d < dim; ++d)
3324 vertices[1][d] =
3325 std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
3326
3327 // read the remaining vertices from the
3328 // list
3329 for (unsigned int v = 2; v < n_vertices + 1; ++v)
3330 {
3331 for (unsigned int i = 0; i < n_vars; ++i)
3332 in >> vars[i];
3333 // fill the vertex
3334 // coordinates. respect the position
3335 // of coordinates in the list of
3336 // variables
3337 for (unsigned int i = 0; i < dim; ++i)
3338 vertices[v][i] = vars[tecplot2deal[i]];
3339 }
3340 }
3341
3342 if (structured)
3343 {
3344 // this is the part of the code that only
3345 // works in 2d
3346 unsigned int I = IJK[0], J = IJK[1];
3347
3348 unsigned int cell = 0;
3349 // set up array of cells
3350 for (unsigned int j = 0; j < J - 1; ++j)
3351 for (unsigned int i = 1; i < I; ++i)
3352 {
3353 cells[cell].vertices[0] = i + j * I;
3354 cells[cell].vertices[1] = i + 1 + j * I;
3355 cells[cell].vertices[2] = i + (j + 1) * I;
3356 cells[cell].vertices[3] = i + 1 + (j + 1) * I;
3357 ++cell;
3358 }
3359 Assert(cell == n_cells, ExcInternalError());
3360 std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
3361 unsigned int k = 0;
3362 for (unsigned int i = 1; i < I + 1; ++i)
3363 {
3364 boundary_vertices[k] = i;
3365 ++k;
3366 boundary_vertices[k] = i + (J - 1) * I;
3367 ++k;
3368 }
3369 for (unsigned int j = 1; j < J - 1; ++j)
3370 {
3371 boundary_vertices[k] = 1 + j * I;
3372 ++k;
3373 boundary_vertices[k] = I + j * I;
3374 ++k;
3375 }
3376 Assert(k == boundary_vertices.size(), ExcInternalError());
3377 // delete the duplicated vertices at the
3378 // boundary, which occur, e.g. in c-type
3379 // or o-type grids around a body
3380 // (airfoil). this automatically deletes
3381 // unused vertices as well.
3383 cells,
3384 subcelldata,
3385 boundary_vertices);
3386 }
3387 else
3388 {
3389 // set up array of cells, unstructured
3390 // mode, so the connectivity is
3391 // explicitly given
3392 for (unsigned int i = 0; i < n_cells; ++i)
3393 {
3394 // note that since in the input file
3395 // we found the number of cells at
3396 // the top, there should still be
3397 // input here, so check this:
3398 AssertThrow(in.fail() == false, ExcIO());
3399
3400 // get the connectivity from the
3401 // input file. the vertices are
3402 // ordered like in the ucd format
3403 for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
3404 in >> cells[i].vertices[GeometryInfo<dim>::ucd_to_deal[j]];
3405 }
3406 }
3407 AssertThrow(in.fail() == false, ExcIO());
3408
3409 apply_grid_fixup_functions(vertices, cells, subcelldata);
3410 tria->create_triangulation(vertices, cells, subcelldata);
3411}
3412
3413
3414
3415template <int dim, int spacedim>
3416void
3421
3422
3423
3424template <int dim, int spacedim>
3425void
3426GridIn<dim, spacedim>::read_assimp(const std::string &filename,
3427 const unsigned int mesh_index,
3428 const bool remove_duplicates,
3429 const double tol,
3430 const bool ignore_unsupported_types)
3431{
3432#ifdef DEAL_II_WITH_ASSIMP
3433 // Only good for surface grids.
3434 AssertThrow(dim < 3, ExcImpossibleInDim(dim));
3435
3436 // Create an instance of the Importer class
3437 Assimp::Importer importer;
3438
3439 // And have it read the given file with some postprocessing
3440 const aiScene *scene =
3441 importer.ReadFile(filename.c_str(),
3442 aiProcess_RemoveComponent |
3443 aiProcess_JoinIdenticalVertices |
3444 aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
3445 aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
3446
3447 // If the import failed, report it
3448 AssertThrow(scene != nullptr, ExcMessage(importer.GetErrorString()));
3449
3450 AssertThrow(scene->mNumMeshes != 0,
3451 ExcMessage("Input file contains no meshes."));
3452
3454 (mesh_index < scene->mNumMeshes),
3455 ExcMessage("Too few meshes in the file."));
3456
3457 unsigned int start_mesh =
3458 (mesh_index == numbers::invalid_unsigned_int ? 0 : mesh_index);
3459 unsigned int end_mesh =
3460 (mesh_index == numbers::invalid_unsigned_int ? scene->mNumMeshes :
3461 mesh_index + 1);
3462
3463 // Deal.II objects are created empty, and then filled with imported file.
3464 std::vector<Point<spacedim>> vertices;
3465 std::vector<CellData<dim>> cells;
3466 SubCellData subcelldata;
3467
3468 // A series of counters to merge cells.
3469 unsigned int v_offset = 0;
3470 unsigned int c_offset = 0;
3471
3472 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
3473 {0, 1, 5, 4, 2, 3, 7, 6}};
3474 // The index of the mesh will be used as a material index.
3475 for (unsigned int m = start_mesh; m < end_mesh; ++m)
3476 {
3477 const aiMesh *mesh = scene->mMeshes[m];
3478
3479 // Check that we know what to do with this mesh, otherwise just
3480 // ignore it
3481 if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
3482 {
3483 AssertThrow(ignore_unsupported_types,
3484 ExcMessage("Incompatible mesh " + std::to_string(m) +
3485 "/" + std::to_string(scene->mNumMeshes)));
3486 continue;
3487 }
3488 else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
3489 {
3490 AssertThrow(ignore_unsupported_types,
3491 ExcMessage("Incompatible mesh " + std::to_string(m) +
3492 "/" + std::to_string(scene->mNumMeshes)));
3493 continue;
3494 }
3495 // Vertices
3496 const unsigned int n_vertices = mesh->mNumVertices;
3497 const aiVector3D *mVertices = mesh->mVertices;
3498
3499 // Faces
3500 const unsigned int n_faces = mesh->mNumFaces;
3501 const aiFace *mFaces = mesh->mFaces;
3502
3503 vertices.resize(v_offset + n_vertices);
3504 cells.resize(c_offset + n_faces);
3505
3506 for (unsigned int i = 0; i < n_vertices; ++i)
3507 for (unsigned int d = 0; d < spacedim; ++d)
3508 vertices[i + v_offset][d] = mVertices[i][d];
3509
3510 unsigned int valid_cell = c_offset;
3511 for (unsigned int i = 0; i < n_faces; ++i)
3512 {
3513 if (mFaces[i].mNumIndices == GeometryInfo<dim>::vertices_per_cell)
3514 {
3515 for (const unsigned int f : GeometryInfo<dim>::vertex_indices())
3516 {
3517 cells[valid_cell]
3518 .vertices[dim == 3 ? local_vertex_numbering[f] :
3520 mFaces[i].mIndices[f] + v_offset;
3521 }
3522 cells[valid_cell].material_id = m;
3523 ++valid_cell;
3524 }
3525 else
3526 {
3527 AssertThrow(ignore_unsupported_types,
3528 ExcMessage("Face " + std::to_string(i) + " of mesh " +
3529 std::to_string(m) + " has " +
3530 std::to_string(mFaces[i].mNumIndices) +
3531 " vertices. We expected only " +
3532 std::to_string(
3534 }
3535 }
3536 cells.resize(valid_cell);
3537
3538 // The vertices are added all at once. Cells are checked for
3539 // validity, so only valid_cells are now present in the deal.II
3540 // list of cells.
3541 v_offset += n_vertices;
3542 c_offset = valid_cell;
3543 }
3544
3545 // No cells were read
3546 if (cells.empty())
3547 return;
3548
3549 if (remove_duplicates)
3550 {
3551 // The function delete_duplicated_vertices() needs to be called more
3552 // than once if a vertex is duplicated more than once. So we keep
3553 // calling it until the number of vertices does not change any more.
3554 unsigned int n_verts = 0;
3555 while (n_verts != vertices.size())
3556 {
3557 n_verts = vertices.size();
3558 std::vector<unsigned int> considered_vertices;
3560 vertices, cells, subcelldata, considered_vertices, tol);
3561 }
3562 }
3563
3564 apply_grid_fixup_functions(vertices, cells, subcelldata);
3565 tria->create_triangulation(vertices, cells, subcelldata);
3566
3567#else
3568 (void)filename;
3569 (void)mesh_index;
3570 (void)remove_duplicates;
3571 (void)tol;
3572 (void)ignore_unsupported_types;
3573 AssertThrow(false, ExcNeedsAssimp());
3574#endif
3575}
3576
3577#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3578// Namespace containing some extra functions for reading ExodusII files
3579namespace
3580{
3581 // Convert ExodusII strings to cell types. Use the number of nodes per
3582 // element to disambiguate some cases.
3584 exodusii_name_to_type(const std::string &type_name,
3585 const int n_nodes_per_element)
3586 {
3587 Assert(type_name.size() > 0, ExcInternalError());
3588 // Try to canonify the name by switching to upper case and removing
3589 // trailing numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and
3590 // PYRAMID13 all equal.
3591 std::string type_name_2 = type_name;
3592 std::transform(type_name_2.begin(),
3593 type_name_2.end(),
3594 type_name_2.begin(),
3595 [](unsigned char c) { return std::toupper(c); });
3596 const std::string numbers = "0123456789";
3597 type_name_2.erase(std::find_first_of(type_name_2.begin(),
3598 type_name_2.end(),
3599 numbers.begin(),
3600 numbers.end()),
3601 type_name_2.end());
3602
3603 // The manual specifies BAR, BEAM, and TRUSS: in practice people use EDGE
3604 if (type_name_2 == "BAR" || type_name_2 == "BEAM" ||
3605 type_name_2 == "EDGE" || type_name_2 == "TRUSS")
3606 return ReferenceCells::Line;
3607 else if (type_name_2 == "TRI" || type_name_2 == "TRIANGLE")
3609 else if (type_name_2 == "QUAD" || type_name_2 == "QUADRILATERAL")
3611 else if (type_name_2 == "SHELL")
3612 {
3613 if (n_nodes_per_element == 3)
3615 else
3617 }
3618 else if (type_name_2 == "TET" || type_name_2 == "TETRA" ||
3619 type_name_2 == "TETRAHEDRON")
3621 else if (type_name_2 == "PYRA" || type_name_2 == "PYRAMID")
3623 else if (type_name_2 == "WEDGE")
3624 return ReferenceCells::Wedge;
3625 else if (type_name_2 == "HEX" || type_name_2 == "HEXAHEDRON")
3627
3630 }
3631
3632 // Associate deal.II boundary ids with sidesets (a face can be in multiple
3633 // sidesets - to translate we assign each set of side set ids to a
3634 // boundary_id or manifold_id)
3635 template <int dim, int spacedim = dim>
3636 std::pair<SubCellData, std::vector<std::vector<int>>>
3637 read_exodusii_sidesets(const int ex_id,
3638 const int n_side_sets,
3639 const std::vector<CellData<dim>> &cells,
3640 const bool apply_all_indicators_to_manifolds)
3641 {
3642 SubCellData subcelldata;
3643 std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
3644 // boundary id 0 is the default
3645 b_or_m_id_to_sideset_ids.emplace_back();
3646 // deal.II does not support assigning boundary ids with nonzero
3647 // codimension meshes so completely skip this information in that case.
3648 //
3649 // Exodus prints warnings if we try to get empty sets so always check
3650 // first
3651 if (dim == spacedim && n_side_sets > 0)
3652 {
3653 std::vector<int> side_set_ids(n_side_sets);
3654 int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
3655 AssertThrowExodusII(ierr);
3656
3657 // First collect all side sets on all boundary faces (indexed here as
3658 // max_faces_per_cell * cell_n + face_n). We then sort and uniquify
3659 // the side sets so that we can convert a set of side set indices into
3660 // a single deal.II boundary or manifold id (and save the
3661 // correspondence).
3662 constexpr auto max_faces_per_cell = GeometryInfo<dim>::faces_per_cell;
3663 std::map<std::size_t, std::vector<int>> face_side_sets;
3664 for (const int side_set_id : side_set_ids)
3665 {
3666 int n_sides = -1;
3667 int n_distribution_factors = -1;
3668
3669 ierr = ex_get_set_param(ex_id,
3670 EX_SIDE_SET,
3671 side_set_id,
3672 &n_sides,
3673 &n_distribution_factors);
3674 AssertThrowExodusII(ierr);
3675 if (n_sides > 0)
3676 {
3677 std::vector<int> elements(n_sides);
3678 std::vector<int> faces(n_sides);
3679 ierr = ex_get_set(ex_id,
3680 EX_SIDE_SET,
3681 side_set_id,
3682 elements.data(),
3683 faces.data());
3684 AssertThrowExodusII(ierr);
3685
3686 // According to the manual (subsection 4.8): "The internal
3687 // number of an element numbering is defined implicitly by the
3688 // order in which it appears in the file. Elements are
3689 // numbered internally (beginning with 1) consecutively across
3690 // all element blocks." Hence element i in Exodus numbering is
3691 // entry i - 1 in the cells array.
3692 for (int side_n = 0; side_n < n_sides; ++side_n)
3693 {
3694 const long element_n = elements[side_n] - 1;
3695 const long face_n = faces[side_n] - 1;
3696 const std::size_t face_id =
3697 element_n * max_faces_per_cell + face_n;
3698 face_side_sets[face_id].push_back(side_set_id);
3699 }
3700 }
3701 }
3702
3703 // Collect into a sortable data structure:
3704 std::vector<std::pair<std::size_t, std::vector<int>>>
3705 face_id_to_side_sets;
3706 face_id_to_side_sets.reserve(face_side_sets.size());
3707 for (auto &pair : face_side_sets)
3708 {
3709 Assert(pair.second.size() > 0, ExcInternalError());
3710 face_id_to_side_sets.emplace_back(std::move(pair));
3711 }
3712
3713 // sort by side sets:
3714 std::sort(face_id_to_side_sets.begin(),
3715 face_id_to_side_sets.end(),
3716 [](const auto &a, const auto &b) {
3717 return std::lexicographical_compare(a.second.begin(),
3718 a.second.end(),
3719 b.second.begin(),
3720 b.second.end());
3721 });
3722
3723 if (dim == 2)
3724 subcelldata.boundary_lines.reserve(face_id_to_side_sets.size());
3725 else if (dim == 3)
3726 subcelldata.boundary_quads.reserve(face_id_to_side_sets.size());
3727 types::boundary_id current_b_or_m_id = 0;
3728 for (const auto &pair : face_id_to_side_sets)
3729 {
3730 const std::size_t face_id = pair.first;
3731 const std::vector<int> &face_sideset_ids = pair.second;
3732 if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
3733 {
3734 // Since we sorted by sideset ids we are guaranteed that if
3735 // this doesn't match the last set then it has not yet been
3736 // seen
3737 ++current_b_or_m_id;
3738 b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
3739 Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
3741 }
3742 // Record the b_or_m_id of the current face.
3743 const unsigned int local_face_n = face_id % max_faces_per_cell;
3744 const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
3745 const ReferenceCell cell_type =
3747 const unsigned int deal_face_n =
3748 cell_type.exodusii_face_to_deal_face(local_face_n);
3749 const ReferenceCell face_reference_cell =
3750 cell_type.face_reference_cell(deal_face_n);
3751
3752 // The orientation we pick doesn't matter here since when we
3753 // create the Triangulation we will sort the vertices for each
3754 // CellData object created here.
3755 if (dim == 2)
3756 {
3757 CellData<1> boundary_line(face_reference_cell.n_vertices());
3758 if (apply_all_indicators_to_manifolds)
3759 boundary_line.manifold_id = current_b_or_m_id;
3760 else
3761 boundary_line.boundary_id = current_b_or_m_id;
3762 for (unsigned int j = 0; j < face_reference_cell.n_vertices();
3763 ++j)
3764 boundary_line.vertices[j] =
3765 cell.vertices[cell_type.face_to_cell_vertices(
3766 deal_face_n, j, 0)];
3767
3768 subcelldata.boundary_lines.push_back(std::move(boundary_line));
3769 }
3770 else if (dim == 3)
3771 {
3772 CellData<2> boundary_quad(face_reference_cell.n_vertices());
3773 if (apply_all_indicators_to_manifolds)
3774 boundary_quad.manifold_id = current_b_or_m_id;
3775 else
3776 boundary_quad.boundary_id = current_b_or_m_id;
3777 for (unsigned int j = 0; j < face_reference_cell.n_vertices();
3778 ++j)
3779 boundary_quad.vertices[j] =
3780 cell.vertices[cell_type.face_to_cell_vertices(
3781 deal_face_n, j, 0)];
3782
3783 subcelldata.boundary_quads.push_back(std::move(boundary_quad));
3784 }
3785 }
3786 }
3787
3788 return std::make_pair(std::move(subcelldata),
3789 std::move(b_or_m_id_to_sideset_ids));
3790 }
3791} // namespace
3792#endif
3793
3794template <int dim, int spacedim>
3797 const std::string &filename,
3798 const bool apply_all_indicators_to_manifolds)
3799{
3800#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3801 // deal.II always uses double precision numbers for geometry
3802 int component_word_size = sizeof(double);
3803 // setting to zero uses the stored word size
3804 int floating_point_word_size = 0;
3805 float ex_version = 0.0;
3806
3807 const int ex_id = ex_open(filename.c_str(),
3808 EX_READ,
3809 &component_word_size,
3810 &floating_point_word_size,
3811 &ex_version);
3812 AssertThrow(ex_id > 0,
3813 ExcMessage("ExodusII failed to open the specified input file."));
3814
3815 // Read basic mesh information:
3816 std::vector<char> string_temp(MAX_LINE_LENGTH + 1, '\0');
3817 int mesh_dimension = 0;
3818 int n_nodes = 0;
3819 int n_elements = 0;
3820 int n_element_blocks = 0;
3821 int n_node_sets = 0;
3822 int n_side_sets = 0;
3823
3824 int ierr = ex_get_init(ex_id,
3825 string_temp.data(),
3826 &mesh_dimension,
3827 &n_nodes,
3828 &n_elements,
3829 &n_element_blocks,
3830 &n_node_sets,
3831 &n_side_sets);
3832 AssertThrowExodusII(ierr);
3833 AssertDimension(mesh_dimension, spacedim);
3834
3835 // Read nodes:
3836 //
3837 // Even if there is a node numbering array the values stored inside the
3838 // ExodusII file must use the contiguous, internal ordering (see Section 4.5
3839 // of the manual - "Internal (contiguously numbered) node and element IDs
3840 // must be used for all data structures that contain node or element numbers
3841 // (IDs), including node set node lists, side set element lists, and element
3842 // connectivity.")
3843 std::vector<Point<spacedim>> vertices;
3844 vertices.reserve(n_nodes);
3845 {
3846 std::vector<double> xs(n_nodes);
3847 std::vector<double> ys(n_nodes);
3848 std::vector<double> zs(n_nodes);
3849
3850 ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
3851 AssertThrowExodusII(ierr);
3852
3853 for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
3854 {
3855 switch (spacedim)
3856 {
3857 case 1:
3858 vertices.emplace_back(xs[vertex_n]);
3859 break;
3860 case 2:
3861 vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
3862 break;
3863 case 3:
3864 vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
3865 break;
3866 default:
3867 Assert(spacedim <= 3, ExcNotImplemented());
3868 }
3869 }
3870 }
3871
3872 std::vector<int> element_block_ids(n_element_blocks);
3873 ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
3874 AssertThrowExodusII(ierr);
3875
3876 std::vector<CellData<dim>> cells;
3877 cells.reserve(n_elements);
3878 // Elements are grouped together by same reference cell type in element
3879 // blocks. There may be multiple blocks for a single reference cell type,
3880 // but "each element block may contain only one element type".
3881 for (const int element_block_id : element_block_ids)
3882 {
3883 std::fill(string_temp.begin(), string_temp.end(), '\0');
3884 int n_block_elements = 0;
3885 int n_nodes_per_element = 0;
3886 int n_edges_per_element = 0;
3887 int n_faces_per_element = 0;
3888 int n_attributes_per_element = 0;
3889
3890 // Extract element data.
3891 ierr = ex_get_block(ex_id,
3892 EX_ELEM_BLOCK,
3893 element_block_id,
3894 string_temp.data(),
3895 &n_block_elements,
3896 &n_nodes_per_element,
3897 &n_edges_per_element,
3898 &n_faces_per_element,
3899 &n_attributes_per_element);
3900 AssertThrowExodusII(ierr);
3901 const ReferenceCell type =
3902 exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
3903 AssertThrow(type.get_dimension() == dim,
3904 ExcMessage(
3905 "The ExodusII block " + std::to_string(element_block_id) +
3906 " with element type " + std::string(string_temp.data()) +
3907 " has dimension " + std::to_string(type.get_dimension()) +
3908 ", which does not match the topological mesh dimension " +
3909 std::to_string(dim) + "."));
3910
3911 // The number of nodes per element may be larger than what we want to
3912 // read - for example, if the Exodus file contains a QUAD9 element, we
3913 // only want to read the first four values and ignore the rest.
3914 Assert(int(type.n_vertices()) <= n_nodes_per_element, ExcInternalError());
3915
3916 std::vector<int> connection(n_nodes_per_element * n_block_elements);
3917 ierr = ex_get_conn(ex_id,
3918 EX_ELEM_BLOCK,
3919 element_block_id,
3920 connection.data(),
3921 nullptr,
3922 nullptr);
3923 AssertThrowExodusII(ierr);
3924
3925 for (unsigned int elem_n = 0; elem_n < connection.size();
3926 elem_n += n_nodes_per_element)
3927 {
3928 CellData<dim> cell(type.n_vertices());
3929 for (const unsigned int i : type.vertex_indices())
3930 {
3932 connection[elem_n + i] - 1;
3933 }
3934 cell.material_id = element_block_id;
3935 cells.push_back(std::move(cell));
3936 }
3937 }
3938
3939 // Extract boundary data.
3940 auto pair = read_exodusii_sidesets<dim, spacedim>(
3941 ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
3942 ierr = ex_close(ex_id);
3943 AssertThrowExodusII(ierr);
3944
3945 apply_grid_fixup_functions(vertices, cells, pair.first);
3946 tria->create_triangulation(vertices, cells, pair.first);
3947 ExodusIIData out;
3948 out.id_to_sideset_ids = std::move(pair.second);
3949 return out;
3950#else
3951 (void)filename;
3952 (void)apply_all_indicators_to_manifolds;
3953 AssertThrow(false, ExcNeedsExodusII());
3954 return {};
3955#endif
3956}
3957
3958
3959template <int dim, int spacedim>
3960void
3962{
3963 std::string line;
3964 while (in)
3965 {
3966 // get line
3967 getline(in, line);
3968
3969 // check if this is a line that
3970 // consists only of spaces, and
3971 // if not put the whole thing
3972 // back and return
3973 if (std::find_if(line.begin(), line.end(), [](const char c) {
3974 return c != ' ';
3975 }) != line.end())
3976 {
3977 in.putback('\n');
3978 for (int i = line.size() - 1; i >= 0; --i)
3979 in.putback(line[i]);
3980 return;
3981 }
3982
3983 // else: go on with next line
3984 }
3985}
3986
3987
3988
3989template <int dim, int spacedim>
3990void
3992 const char comment_start)
3993{
3994 char c;
3995 // loop over the following comment
3996 // lines
3997 while (in.get(c) && c == comment_start)
3998 // loop over the characters after
3999 // the comment starter
4000 while (in.get() != '\n')
4001 ;
4002
4003
4004 // put back first character of
4005 // first non-comment line
4006 if (in)
4007 in.putback(c);
4008
4009 // at last: skip additional empty lines, if present
4010 skip_empty_lines(in);
4011}
4012
4013
4014
4015template <int dim, int spacedim>
4016void
4018 const std::vector<CellData<dim>> & /*cells*/,
4019 const std::vector<Point<spacedim>> & /*vertices*/,
4020 std::ostream & /*out*/)
4021{
4023}
4024
4025
4026
4027template <>
4028void
4030 const std::vector<Point<2>> &vertices,
4031 std::ostream &out)
4032{
4033 double min_x = vertices[cells[0].vertices[0]][0],
4034 max_x = vertices[cells[0].vertices[0]][0],
4035 min_y = vertices[cells[0].vertices[0]][1],
4036 max_y = vertices[cells[0].vertices[0]][1];
4037
4038 for (unsigned int i = 0; i < cells.size(); ++i)
4039 {
4040 for (const auto vertex : cells[i].vertices)
4041 {
4042 const Point<2> &p = vertices[vertex];
4043
4044 if (p[0] < min_x)
4045 min_x = p[0];
4046 if (p[0] > max_x)
4047 max_x = p[0];
4048 if (p[1] < min_y)
4049 min_y = p[1];
4050 if (p[1] > max_y)
4051 max_y = p[1];
4052 }
4053
4054 out << "# cell " << i << std::endl;
4056 for (const auto vertex : cells[i].vertices)
4057 center += vertices[vertex];
4058 center /= 4;
4059
4060 out << "set label \"" << i << "\" at " << center[0] << ',' << center[1]
4061 << " center" << std::endl;
4062
4063 // first two line right direction
4064 for (unsigned int f = 0; f < 2; ++f)
4065 out << "set arrow from " << vertices[cells[i].vertices[f]][0] << ','
4066 << vertices[cells[i].vertices[f]][1] << " to "
4067 << vertices[cells[i].vertices[(f + 1) % 4]][0] << ','
4068 << vertices[cells[i].vertices[(f + 1) % 4]][1] << std::endl;
4069 // other two lines reverse direction
4070 for (unsigned int f = 2; f < 4; ++f)
4071 out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]][0]
4072 << ',' << vertices[cells[i].vertices[(f + 1) % 4]][1] << " to "
4073 << vertices[cells[i].vertices[f]][0] << ','
4074 << vertices[cells[i].vertices[f]][1] << std::endl;
4075 out << std::endl;
4076 }
4077
4078
4079 out << std::endl
4080 << "set nokey" << std::endl
4081 << "pl [" << min_x << ':' << max_x << "][" << min_y << ':' << max_y
4082 << "] " << min_y << std::endl
4083 << "pause -1" << std::endl;
4084}
4085
4086
4087
4088template <>
4089void
4091 const std::vector<Point<3>> &vertices,
4092 std::ostream &out)
4093{
4094 for (const auto &cell : cells)
4095 {
4096 // line 0
4097 out << vertices[cell.vertices[0]] << std::endl
4098 << vertices[cell.vertices[1]] << std::endl
4099 << std::endl
4100 << std::endl;
4101 // line 1
4102 out << vertices[cell.vertices[1]] << std::endl
4103 << vertices[cell.vertices[2]] << std::endl
4104 << std::endl
4105 << std::endl;
4106 // line 2
4107 out << vertices[cell.vertices[3]] << std::endl
4108 << vertices[cell.vertices[2]] << std::endl
4109 << std::endl
4110 << std::endl;
4111 // line 3
4112 out << vertices[cell.vertices[0]] << std::endl
4113 << vertices[cell.vertices[3]] << std::endl
4114 << std::endl
4115 << std::endl;
4116 // line 4
4117 out << vertices[cell.vertices[4]] << std::endl
4118 << vertices[cell.vertices[5]] << std::endl
4119 << std::endl
4120 << std::endl;
4121 // line 5
4122 out << vertices[cell.vertices[5]] << std::endl
4123 << vertices[cell.vertices[6]] << std::endl
4124 << std::endl
4125 << std::endl;
4126 // line 6
4127 out << vertices[cell.vertices[7]] << std::endl
4128 << vertices[cell.vertices[6]] << std::endl
4129 << std::endl
4130 << std::endl;
4131 // line 7
4132 out << vertices[cell.vertices[4]] << std::endl
4133 << vertices[cell.vertices[7]] << std::endl
4134 << std::endl
4135 << std::endl;
4136 // line 8
4137 out << vertices[cell.vertices[0]] << std::endl
4138 << vertices[cell.vertices[4]] << std::endl
4139 << std::endl
4140 << std::endl;
4141 // line 9
4142 out << vertices[cell.vertices[1]] << std::endl
4143 << vertices[cell.vertices[5]] << std::endl
4144 << std::endl
4145 << std::endl;
4146 // line 10
4147 out << vertices[cell.vertices[2]] << std::endl
4148 << vertices[cell.vertices[6]] << std::endl
4149 << std::endl
4150 << std::endl;
4151 // line 11
4152 out << vertices[cell.vertices[3]] << std::endl
4153 << vertices[cell.vertices[7]] << std::endl
4154 << std::endl
4155 << std::endl;
4156 }
4157}
4158
4159
4160
4161template <int dim, int spacedim>
4162void
4163GridIn<dim, spacedim>::read(const std::string &filename, Format format)
4164{
4165 // Check early that the file actually exists and if not throw ExcFileNotOpen.
4166 AssertThrow(std::filesystem::exists(filename), ExcFileNotOpen(filename));
4167
4168 if (format == Default)
4169 {
4170 const std::string::size_type slashpos = filename.find_last_of('/');
4171 const std::string::size_type dotpos = filename.find_last_of('.');
4172 if (dotpos < filename.size() &&
4173 (dotpos > slashpos || slashpos == std::string::npos))
4174 {
4175 std::string ext = filename.substr(dotpos + 1);
4176 format = parse_format(ext);
4177 }
4178 }
4179
4180 if (format == assimp)
4181 {
4182 read_assimp(filename);
4183 }
4184 else if (format == exodusii)
4185 {
4186 read_exodusii(filename);
4187 }
4188 else
4189 {
4190 std::ifstream in(filename);
4191 read(in, format);
4192 }
4193}
4194
4195
4196template <int dim, int spacedim>
4197void
4198GridIn<dim, spacedim>::read(std::istream &in, Format format)
4199{
4200 if (format == Default)
4201 format = default_format;
4202
4203 switch (format)
4204 {
4205 case dbmesh:
4206 read_dbmesh(in);
4207 return;
4208
4209 case msh:
4210 read_msh(in);
4211 return;
4212
4213 case vtk:
4214 read_vtk(in);
4215 return;
4216
4217 case vtu:
4218 read_vtu(in);
4219 return;
4220
4221 case unv:
4222 read_unv(in);
4223 return;
4224
4225 case ucd:
4226 read_ucd(in);
4227 return;
4228
4229 case abaqus:
4230 read_abaqus(in);
4231 return;
4232
4233 case xda:
4234 read_xda(in);
4235 return;
4236
4237 case tecplot:
4238 read_tecplot(in);
4239 return;
4240
4241 case assimp:
4242 Assert(false,
4243 ExcMessage("There is no read_assimp(istream &) function. "
4244 "Use the read_assimp(string &filename, ...) "
4245 "functions, instead."));
4246 return;
4247
4248 case exodusii:
4249 Assert(false,
4250 ExcMessage("There is no read_exodusii(istream &) function. "
4251 "Use the read_exodusii(string &filename, ...) "
4252 "function, instead."));
4253 return;
4254
4255 case Default:
4256 break;
4257 }
4259}
4260
4261
4262
4263template <int dim, int spacedim>
4264std::string
4266{
4267 switch (format)
4268 {
4269 case dbmesh:
4270 return ".dbmesh";
4271 case exodusii:
4272 return ".e";
4273 case msh:
4274 return ".msh";
4275 case vtk:
4276 return ".vtk";
4277 case vtu:
4278 return ".vtu";
4279 case unv:
4280 return ".unv";
4281 case ucd:
4282 return ".inp";
4283 case abaqus:
4284 return ".inp"; // Typical suffix for Abaqus mesh files conflicts with
4285 // UCD.
4286 case xda:
4287 return ".xda";
4288 case tecplot:
4289 return ".dat";
4290 default:
4292 return ".unknown_format";
4293 }
4294}
4295
4296
4297
4298template <int dim, int spacedim>
4300GridIn<dim, spacedim>::parse_format(const std::string &format_name)
4301{
4302 if (format_name == "dbmesh")
4303 return dbmesh;
4304
4305 if (format_name == "exodusii")
4306 return exodusii;
4307
4308 if (format_name == "msh")
4309 return msh;
4310
4311 if (format_name == "unv")
4312 return unv;
4313
4314 if (format_name == "vtk")
4315 return vtk;
4316
4317 if (format_name == "vtu")
4318 return vtu;
4319
4320 // This is also the typical extension of Abaqus input files.
4321 if (format_name == "inp")
4322 return ucd;
4323
4324 if (format_name == "ucd")
4325 return ucd;
4326
4327 if (format_name == "xda")
4328 return xda;
4329
4330 if (format_name == "tecplot")
4331 return tecplot;
4332
4333 if (format_name == "dat")
4334 return tecplot;
4335
4336 if (format_name == "plt")
4337 // Actually, this is the extension for the
4338 // tecplot binary format, which we do not
4339 // support right now. However, some people
4340 // tend to create tecplot ascii files with
4341 // the extension 'plt' instead of
4342 // 'dat'. Thus, include this extension
4343 // here. If it actually is a binary file,
4344 // the read_tecplot() function will fail
4345 // and throw an exception, anyway.
4346 return tecplot;
4347
4348 AssertThrow(false, ExcInvalidState());
4349 // return something weird
4350 return Format(Default);
4351}
4352
4353
4354
4355template <int dim, int spacedim>
4356std::string
4358{
4359 return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
4360}
4361
4362
4363
4364namespace
4365{
4366 template <int dim, int spacedim>
4367 Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
4368 : tolerance(5e-16) // Used to offset Cubit tolerance error when outputting
4369 // value close to zero
4370 {
4371 AssertThrow(spacedim == 2 || spacedim == 3, ExcNotImplemented());
4372 }
4373
4374
4375
4376 // Convert from a string to some other data type
4377 // Reference: http://www.codeguru.com/forum/showthread.php?t=231054
4378 template <class T>
4379 bool
4380 from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &))
4381 {
4382 std::istringstream iss(s);
4383 return !(iss >> f >> t).fail();
4384 }
4385
4386
4387
4388 // Extract an integer from a string
4389 int
4390 extract_int(const std::string &s)
4391 {
4392 std::string tmp;
4393 for (const char c : s)
4394 {
4395 if (isdigit(c) != 0)
4396 {
4397 tmp += c;
4398 }
4399 }
4400
4401 int number = 0;
4402 from_string(number, tmp, std::dec);
4403 return number;
4404 }
4405
4406
4407
4408 template <int dim, int spacedim>
4409 void
4410 Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
4411 {
4412 // References:
4413 // http://www.egr.msu.edu/software/abaqus/Documentation/docs/v6.7/books/usb/default.htm?startat=pt01ch02.html
4414 // http://www.cprogramming.com/tutorial/string.html
4415
4416 AssertThrow(input_stream.fail() == false, ExcIO());
4417 std::string line;
4418
4419 while (std::getline(input_stream, line))
4420 {
4421 cont:
4422 std::transform(line.begin(), line.end(), line.begin(), ::toupper);
4423
4424 if (line.compare("*HEADING") == 0 || line.compare(0, 2, "**") == 0 ||
4425 line.compare(0, 5, "*PART") == 0)
4426 {
4427 // Skip header and comments
4428 while (std::getline(input_stream, line))
4429 {
4430 if (line[0] == '*')
4431 goto cont; // My eyes, they burn!
4432 }
4433 }
4434 else if (line.compare(0, 5, "*NODE") == 0)
4435 {
4436 // Extract list of vertices
4437 // Header line might be:
4438 // *NODE, NSET=ALLNODES
4439 // *NODE
4440
4441 // Contains lines in the form:
4442 // Index, x, y, z
4443 while (std::getline(input_stream, line))
4444 {
4445 if (line[0] == '*')
4446 goto cont;
4447
4448 std::vector<double> node(spacedim + 1);
4449
4450 std::istringstream iss(line);
4451 char comma;
4452 for (unsigned int i = 0; i < spacedim + 1; ++i)
4453 iss >> node[i] >> comma;
4454
4455 node_list.push_back(node);
4456 }
4457 }
4458 else if (line.compare(0, 8, "*ELEMENT") == 0)
4459 {
4460 // Element construction.
4461 // There are different header formats, the details
4462 // of which we're not particularly interested in except
4463 // whether they represent quads or hexahedrals.
4464 // *ELEMENT, TYPE=S4R, ELSET=EB<material id>
4465 // *ELEMENT, TYPE=C3d8R, ELSET=EB<material id>
4466 // *ELEMENT, TYPE=C3d8
4467 // Elements itself (n=4 or n=8):
4468 // Index, i[0], ..., i[n]
4469
4470 int material = 0;
4471 // Scan for material id
4472 {
4473 const std::string before_material = "ELSET=EB";
4474 const std::size_t idx = line.find(before_material);
4475 if (idx != std::string::npos)
4476 {
4477 from_string(material,
4478 line.substr(idx + before_material.size()),
4479 std::dec);
4480 }
4481 }
4482
4483 // Read ELEMENT definition
4484 while (std::getline(input_stream, line))
4485 {
4486 if (line[0] == '*')
4487 goto cont;
4488
4489 std::istringstream iss(line);
4490 char comma;
4491
4492 // We will store the material id in the zeroth entry of the
4493 // vector and the rest of the elements represent the global
4494 // node numbers
4495 const unsigned int n_data_per_cell =
4497 std::vector<double> cell(n_data_per_cell);
4498 for (unsigned int i = 0; i < n_data_per_cell; ++i)
4499 iss >> cell[i] >> comma;
4500
4501 // Overwrite cell index from file by material
4502 cell[0] = static_cast<double>(material);
4503 cell_list.push_back(cell);
4504 }
4505 }
4506 else if (line.compare(0, 8, "*SURFACE") == 0)
4507 {
4508 // Extract the definitions of boundary surfaces
4509 // Old format from Cubit:
4510 // *SURFACE, NAME=SS<boundary indicator>
4511 // <element index>, S<face number>
4512 // Abaqus default format:
4513 // *SURFACE, TYPE=ELEMENT, NAME=SURF-<indicator>
4514
4515 // Get name of the surface and extract id from it;
4516 // this will be the boundary indicator
4517 const std::string name_key = "NAME=";
4518 const std::size_t name_idx_start =
4519 line.find(name_key) + name_key.size();
4520 std::size_t name_idx_end = line.find(',', name_idx_start);
4521 if (name_idx_end == std::string::npos)
4522 {
4523 name_idx_end = line.size();
4524 }
4525 const int b_indicator = extract_int(
4526 line.substr(name_idx_start, name_idx_end - name_idx_start));
4527
4528 // Read SURFACE definition
4529 // Note that the orientation of the faces is embedded within the
4530 // definition of each "set" of faces that comprise the surface
4531 // These are either marked by an "S" or "E" in 3d or 2d
4532 // respectively.
4533 while (std::getline(input_stream, line))
4534 {
4535 if (line[0] == '*')
4536 goto cont;
4537
4538 // Change all characters to upper case
4539 std::transform(line.begin(),
4540 line.end(),
4541 line.begin(),
4542 ::toupper);
4543
4544 // Surface can be created from ELSET, or directly from cells
4545 // If elsets_list contains a key with specific name - refers
4546 // to that ELSET, otherwise refers to cell
4547 std::istringstream iss(line);
4548 int el_idx;
4549 int face_number;
4550 char temp;
4551
4552 // Get relevant faces, taking into account the element
4553 // orientation
4554 std::vector<double> quad_node_list;
4555 const std::string elset_name = line.substr(0, line.find(','));
4556 if (elsets_list.count(elset_name) != 0)
4557 {
4558 // Surface refers to ELSET
4559 std::string stmp;
4560 iss >> stmp >> temp >> face_number;
4561
4562 const std::vector<int> cells = elsets_list[elset_name];
4563 for (const int cell : cells)
4564 {
4565 el_idx = cell;
4566 quad_node_list =
4567 get_global_node_numbers(el_idx, face_number);
4568 quad_node_list.insert(quad_node_list.begin(),
4569 b_indicator);
4570
4571 face_list.push_back(quad_node_list);
4572 }
4573 }
4574 else
4575 {
4576 // Surface refers directly to elements
4577 char comma;
4578 iss >> el_idx >> comma >> temp >> face_number;
4579 quad_node_list =
4580 get_global_node_numbers(el_idx, face_number);
4581 quad_node_list.insert(quad_node_list.begin(), b_indicator);
4582
4583 face_list.push_back(quad_node_list);
4584 }
4585 }
4586 }
4587 else if (line.compare(0, 6, "*ELSET") == 0)
4588 {
4589 // Get ELSET name.
4590 // Materials are attached to elsets with specific name
4591 std::string elset_name;
4592 {
4593 const std::string elset_key = "*ELSET, ELSET=";
4594 const std::size_t idx = line.find(elset_key);
4595 if (idx != std::string::npos)
4596 {
4597 const std::string comma = ",";
4598 const std::size_t first_comma = line.find(comma);
4599 const std::size_t second_comma =
4600 line.find(comma, first_comma + 1);
4601 const std::size_t elset_name_start =
4602 line.find(elset_key) + elset_key.size();
4603 elset_name = line.substr(elset_name_start,
4604 second_comma - elset_name_start);
4605 }
4606 }
4607
4608 // There are two possibilities of storing cells numbers in ELSET:
4609 // 1. If the header contains the 'GENERATE' keyword, then the next
4610 // line describes range of cells as:
4611 // cell_id_start, cell_id_end, cell_step
4612 // 2. If the header does not contain the 'GENERATE' keyword, then
4613 // the next lines contain cells numbers
4614 std::vector<int> elements;
4615 const std::size_t generate_idx = line.find("GENERATE");
4616 if (generate_idx != std::string::npos)
4617 {
4618 // Option (1)
4619 std::getline(input_stream, line);
4620 std::istringstream iss(line);
4621 char comma;
4622 int elid_start;
4623 int elid_end;
4624 int elis_step = 1; // Default if case stride not provided
4625
4626 // Some files don't have the stride size
4627 // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp
4628 // to
4629 // ./grids/abaqus/2d/2d_test_abaqus.inp
4630 iss >> elid_start >> comma >> elid_end;
4631 AssertThrow(comma == ',',
4632 ExcMessage(
4633 std::string(
4634 "While reading an ABAQUS file, the reader "
4635 "expected a comma but found a <") +
4636 comma + "> in the line <" + line + ">."));
4638 elid_start <= elid_end,
4639 ExcMessage(
4640 std::string(
4641 "While reading an ABAQUS file, the reader encountered "
4642 "a GENERATE statement in which the upper bound <") +
4643 Utilities::int_to_string(elid_end) +
4644 "> for the element numbers is not larger or equal "
4645 "than the lower bound <" +
4646 Utilities::int_to_string(elid_start) + ">."));
4647
4648 // https://stackoverflow.com/questions/8046357/how-do-i-check-if-a-stringstream-variable-is-empty-null
4649 if (iss.rdbuf()->in_avail() != 0)
4650 iss >> comma >> elis_step;
4651 AssertThrow(comma == ',',
4652 ExcMessage(
4653 std::string(
4654 "While reading an ABAQUS file, the reader "
4655 "expected a comma but found a <") +
4656 comma + "> in the line <" + line + ">."));
4657
4658 for (int i = elid_start; i <= elid_end; i += elis_step)
4659 elements.push_back(i);
4660 elsets_list[elset_name] = elements;
4661
4662 std::getline(input_stream, line);
4663 }
4664 else
4665 {
4666 // Option (2)
4667 while (std::getline(input_stream, line))
4668 {
4669 if (line[0] == '*')
4670 break;
4671
4672 std::istringstream iss(line);
4673 char comma;
4674 int elid;
4675 while (!iss.eof())
4676 {
4677 iss >> elid >> comma;
4679 comma == ',',
4680 ExcMessage(
4681 std::string(
4682 "While reading an ABAQUS file, the reader "
4683 "expected a comma but found a <") +
4684 comma + "> in the line <" + line + ">."));
4685
4686 elements.push_back(elid);
4687 }
4688 }
4689
4690 elsets_list[elset_name] = elements;
4691 }
4692
4693 goto cont;
4694 }
4695 else if (line.compare(0, 5, "*NSET") == 0)
4696 {
4697 // Skip nodesets; we have no use for them
4698 while (std::getline(input_stream, line))
4699 {
4700 if (line[0] == '*')
4701 goto cont;
4702 }
4703 }
4704 else if (line.compare(0, 14, "*SOLID SECTION") == 0)
4705 {
4706 // The ELSET name, which describes a section for particular
4707 // material
4708 const std::string elset_key = "ELSET=";
4709 const std::size_t elset_start =
4710 line.find("ELSET=") + elset_key.size();
4711 const std::size_t elset_end = line.find(',', elset_start + 1);
4712 const std::string elset_name =
4713 line.substr(elset_start, elset_end - elset_start);
4714
4715 // Solid material definition.
4716 // We assume that material id is taken from material name,
4717 // eg. "Material-1" -> ID=1
4718 const std::string material_key = "MATERIAL=";
4719 const std::size_t last_equal =
4720 line.find("MATERIAL=") + material_key.size();
4721 const std::size_t material_id_start = line.find('-', last_equal);
4722 int material_id = 0;
4723 from_string(material_id,
4724 line.substr(material_id_start + 1),
4725 std::dec);
4726
4727 // Assign material id to cells
4728 const std::vector<int> &elset_cells = elsets_list[elset_name];
4729 for (const int elset_cell : elset_cells)
4730 {
4731 const int cell_id = elset_cell - 1;
4732 cell_list[cell_id][0] = material_id;
4733 }
4734 }
4735 // Note: All other lines / entries are ignored
4736 }
4737 }
4738
4739 template <int dim, int spacedim>
4740 std::vector<double>
4741 Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
4742 const int face_cell_no,
4743 const int face_cell_face_no) const
4744 {
4745 std::vector<double> quad_node_list(GeometryInfo<dim>::vertices_per_face);
4746
4747 // These orderings were reverse engineered by hand and may
4748 // conceivably be erroneous.
4749 // TODO: Currently one test (2d unstructured mesh) in the test
4750 // suite fails, presumably because of an ordering issue.
4751 if (dim == 2)
4752 {
4753 if (face_cell_face_no == 1)
4754 {
4755 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4756 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4757 }
4758 else if (face_cell_face_no == 2)
4759 {
4760 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4761 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4762 }
4763 else if (face_cell_face_no == 3)
4764 {
4765 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4766 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4767 }
4768 else if (face_cell_face_no == 4)
4769 {
4770 quad_node_list[0] = cell_list[face_cell_no - 1][4];
4771 quad_node_list[1] = cell_list[face_cell_no - 1][1];
4772 }
4773 else
4774 {
4775 AssertThrow(face_cell_face_no <= 4,
4776 ExcMessage("Invalid face number in 2d"));
4777 }
4778 }
4779 else if (dim == 3)
4780 {
4781 if (face_cell_face_no == 1)
4782 {
4783 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4784 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4785 quad_node_list[2] = cell_list[face_cell_no - 1][3];
4786 quad_node_list[3] = cell_list[face_cell_no - 1][2];
4787 }
4788 else if (face_cell_face_no == 2)
4789 {
4790 quad_node_list[0] = cell_list[face_cell_no - 1][5];
4791 quad_node_list[1] = cell_list[face_cell_no - 1][8];
4792 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4793 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4794 }
4795 else if (face_cell_face_no == 3)
4796 {
4797 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4798 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4799 quad_node_list[2] = cell_list[face_cell_no - 1][6];
4800 quad_node_list[3] = cell_list[face_cell_no - 1][5];
4801 }
4802 else if (face_cell_face_no == 4)
4803 {
4804 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4805 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4806 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4807 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4808 }
4809 else if (face_cell_face_no == 5)
4810 {
4811 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4812 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4813 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4814 quad_node_list[3] = cell_list[face_cell_no - 1][7];
4815 }
4816 else if (face_cell_face_no == 6)
4817 {
4818 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4819 quad_node_list[1] = cell_list[face_cell_no - 1][5];
4820 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4821 quad_node_list[3] = cell_list[face_cell_no - 1][4];
4822 }
4823 else
4824 {
4825 AssertThrow(face_cell_no <= 6,
4826 ExcMessage("Invalid face number in 3d"));
4827 }
4828 }
4829 else
4830 {
4831 AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
4832 }
4833
4834 return quad_node_list;
4835 }
4836
4837 template <int dim, int spacedim>
4838 void
4839 Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output) const
4840 {
4841 // References:
4842 // http://www.dealii.org/developer/doxygen/deal.II/structGeometryInfo.html
4843 // http://people.scs.fsu.edu/~burkardt/data/ucd/ucd.html
4844
4845 AssertThrow(output.fail() == false, ExcIO());
4846
4847 // save old formatting options
4848 const boost::io::ios_base_all_saver formatting_saver(output);
4849
4850 // Write out title - Note: No other commented text can be inserted below
4851 // the title in a UCD file
4852 output << "# Abaqus to UCD mesh conversion" << std::endl;
4853 output << "# Mesh type: AVS UCD" << std::endl;
4854
4855 // ========================================================
4856 // ASCII UCD File Format
4857 // The input file cannot contain blank lines or lines with leading blanks.
4858 // Comments, if present, must precede all data in the file.
4859 // Comments within the data will cause read errors.
4860 // The general order of the data is as follows:
4861 // 1. Numbers defining the overall structure, including the number of
4862 // nodes,
4863 // the number of cells, and the length of the vector of data associated
4864 // with the nodes, cells, and the model.
4865 // e.g. 1:
4866 // <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
4867 // e.g. 2:
4868 // n_elements = n_hex_cells + n_bc_quads + n_quad_cells +
4869 // n_bc_edges outfile.write(str(n_nodes) + " " + str(n_elements) +
4870 // " 0 0 0\n")
4871 // 2. For each node, its node id and the coordinates of that node in
4872 // space.
4873 // Node-ids must be integers, but any number including non sequential
4874 // numbers can be used. Mid-edge nodes are treated like any other node.
4875 // 3. For each cell: its cell-id, material, cell type (hexahedral,
4876 // pyramid,
4877 // etc.), and the list of node-ids that correspond to each of the
4878 // cell's vertices. The below table specifies the different cell types
4879 // and the keyword used to represent them in the file.
4880
4881 // Write out header
4882 output << node_list.size() << "\t" << (cell_list.size() + face_list.size())
4883 << "\t0\t0\t0" << std::endl;
4884
4885 output.width(16);
4886 output.precision(8);
4887
4888 // Write out node numbers
4889 // Loop over all nodes
4890 for (const auto &node : node_list)
4891 {
4892 // Node number
4893 output << node[0] << "\t";
4894
4895 // Node coordinates
4896 output.setf(std::ios::scientific, std::ios::floatfield);
4897 for (unsigned int jj = 1; jj < spacedim + 1; ++jj)
4898 {
4899 // invoke tolerance -> set points close to zero equal to zero
4900 if (std::abs(node[jj]) > tolerance)
4901 output << static_cast<double>(node[jj]) << "\t";
4902 else
4903 output << 0.0 << "\t";
4904 }
4905 if (spacedim == 2)
4906 output << 0.0 << "\t";
4907
4908 output << std::endl;
4909 output.unsetf(std::ios::floatfield);
4910 }
4911
4912 // Write out cell node numbers
4913 for (unsigned int ii = 0; ii < cell_list.size(); ++ii)
4914 {
4915 output << ii + 1 << "\t" << cell_list[ii][0] << "\t"
4916 << (dim == 2 ? "quad" : "hex") << "\t";
4917 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
4918 ++jj)
4919 output << cell_list[ii][jj] << "\t";
4920
4921 output << std::endl;
4922 }
4923
4924 // Write out quad node numbers
4925 for (unsigned int ii = 0; ii < face_list.size(); ++ii)
4926 {
4927 output << ii + 1 << "\t" << face_list[ii][0] << "\t"
4928 << (dim == 2 ? "line" : "quad") << "\t";
4929 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
4930 ++jj)
4931 output << face_list[ii][jj] << "\t";
4932
4933 output << std::endl;
4934 }
4935 }
4936} // namespace
4937
4938
4939// explicit instantiations
4940#include "grid_in.inst"
4941
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:3961
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:3426
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:4265
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:3991
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:4300
void read_vtu(std::istream &in)
Definition grid_in.cc:587
void read_tecplot(std::istream &in)
Definition grid_in.cc:3417
ExodusIIData read_exodusii(const std::string &filename, const bool apply_all_indicators_to_manifolds=false)
Definition grid_in.cc:3796
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:4198
static void debug_output_grid(const std::vector< CellData< dim > > &cells, const std::vector< Point< spacedim > > &vertices, std::ostream &out)
Definition grid_in.cc:4017
static std::string get_format_names()
Definition grid_in.cc:4357
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:2952
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define DEAL_II_FALLTHROUGH
Definition config.h:233
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:2459
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