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