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