1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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
17 #include <deal.II/base/exceptions.h>
18 #include <deal.II/base/path_search.h>
19 #include <deal.II/base/utilities.h>
20
21 #include <deal.II/grid/grid_in.h>
22 #include <deal.II/grid/grid_reordering.h>
23 #include <deal.II/grid/grid_tools.h>
24 #include <deal.II/grid/tria.h>
25
26 #include <boost/archive/binary_iarchive.hpp>
27 #include <boost/io/ios_state.hpp>
28 #include <boost/property_tree/ptree.hpp>
29 #include <boost/property_tree/xml_parser.hpp>
30 #include <boost/serialization/serialization.hpp>
31
32 #include <algorithm>
33 #include <cctype>
34 #include <fstream>
35 #include <functional>
36 #include <map>
37
38 #ifdef DEAL_II_WITH_ASSIMP
39 # include <assimp/Importer.hpp> // C++ importer interface
40 # include <assimp/postprocess.h> // Post processing flags
41 # include <assimp/scene.h> // Output data structure
42 #endif
43
44
45 DEAL_II_NAMESPACE_OPEN
46
47
48 namespace
49 {
50 /**
51 * In 1d, boundary indicators are associated with vertices, but this is not
52 * currently passed through the SubcellData structure. This function sets
53 * boundary indicators on vertices after the triangulation has already been
54 * created.
55 *
56 * TODO: Fix this properly via SubcellData
57 */
58 template <int spacedim>
59 void
assign_1d_boundary_ids(const std::map<unsigned int,types::boundary_id> & boundary_ids,Triangulation<1,spacedim> & triangulation)60 assign_1d_boundary_ids(
61 const std::map<unsigned int, types::boundary_id> &boundary_ids,
62 Triangulation<1, spacedim> & triangulation)
63 {
64 if (boundary_ids.size() > 0)
65 for (const auto &cell : triangulation.active_cell_iterators())
66 for (unsigned int f : GeometryInfo<1>::face_indices())
67 if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
68 {
69 AssertThrow(
70 cell->at_boundary(f),
71 ExcMessage(
72 "You are trying to prescribe boundary ids on the face "
73 "of a 1d cell (i.e., on a vertex), but this face is not actually at "
74 "the boundary of the mesh. This is not allowed."));
75 cell->face(f)->set_boundary_id(
76 boundary_ids.find(cell->vertex_index(f))->second);
77 }
78 }
79
80
81 template <int dim, int spacedim>
82 void
assign_1d_boundary_ids(const std::map<unsigned int,types::boundary_id> &,Triangulation<dim,spacedim> &)83 assign_1d_boundary_ids(const std::map<unsigned int, types::boundary_id> &,
84 Triangulation<dim, spacedim> &)
85 {
86 // we shouldn't get here since boundary ids are not assigned to
87 // vertices except in 1d
88 Assert(dim != 1, ExcInternalError());
89 }
90 } // namespace
91
92 template <int dim, int spacedim>
GridIn()93 GridIn<dim, spacedim>::GridIn()
94 : tria(nullptr, typeid(*this).name())
95 , default_format(ucd)
96 {}
97
98
99
100 template <int dim, int spacedim>
GridIn(Triangulation<dim,spacedim> & t)101 GridIn<dim, spacedim>::GridIn(Triangulation<dim, spacedim> &t)
102 : tria(&t, typeid(*this).name())
103 , default_format(ucd)
104 {}
105
106
107
108 template <int dim, int spacedim>
109 void
attach_triangulation(Triangulation<dim,spacedim> & t)110 GridIn<dim, spacedim>::attach_triangulation(Triangulation<dim, spacedim> &t)
111 {
112 tria = &t;
113 }
114
115
116
117 template <int dim, int spacedim>
118 void
read_vtk(std::istream & in)119 GridIn<dim, spacedim>::read_vtk(std::istream &in)
120 {
121 std::string line;
122
123 // verify that the first, third and fourth lines match
124 // expectations. the second line of the file may essentially be
125 // anything the author of the file chose to identify what's in
126 // there, so we just ensure that we can read it
127 {
128 std::string text[4];
129 text[0] = "# vtk DataFile Version 3.0";
130 text[1] = "****";
131 text[2] = "ASCII";
132 text[3] = "DATASET UNSTRUCTURED_GRID";
133
134 for (unsigned int i = 0; i < 4; ++i)
135 {
136 getline(in, line);
137 if (i != 1)
138 AssertThrow(
139 line.compare(text[i]) == 0,
140 ExcMessage(
141 std::string(
142 "While reading VTK file, failed to find a header line with text <") +
143 text[i] + ">"));
144 }
145 }
146
147 ///////////////////Declaring storage and mappings//////////////////
148
149 std::vector<Point<spacedim>> vertices;
150 std::vector<CellData<dim>> cells;
151 SubCellData subcelldata;
152
153 std::string keyword;
154
155 in >> keyword;
156
157 //////////////////Processing the POINTS section///////////////
158
159 if (keyword == "POINTS")
160 {
161 unsigned int n_vertices;
162 in >> n_vertices;
163
164 in >> keyword; // float, double, int, char, etc.
165
166 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
167 {
168 // VTK format always specifies vertex coordinates with 3 components
169 Point<3> x;
170 in >> x(0) >> x(1) >> x(2);
171
172 vertices.emplace_back();
173 for (unsigned int d = 0; d < spacedim; ++d)
174 vertices.back()(d) = x(d);
175 }
176 }
177
178 else
179 AssertThrow(false,
180 ExcMessage(
181 "While reading VTK file, failed to find POINTS section"));
182
183 in >> keyword;
184
185 unsigned int n_geometric_objects = 0;
186 unsigned int n_ints;
187
188 bool is_quad_or_hex_mesh = false;
189 bool is_tria_or_tet_mesh = false;
190
191 if (keyword == "CELLS")
192 {
193 // jump to the `CELL_TYPES` section and read in cell types
194 std::vector<unsigned int> cell_types;
195 {
196 std::streampos oldpos = in.tellg();
197
198
199 while (in >> keyword)
200 if (keyword == "CELL_TYPES")
201 {
202 in >> n_ints;
203
204 cell_types.resize(n_ints);
205
206 for (unsigned int i = 0; i < n_ints; ++i)
207 in >> cell_types[i];
208
209 break;
210 }
211
212 in.seekg(oldpos);
213 }
214
215 in >> n_geometric_objects;
216 in >> n_ints; // Ignore this, since we don't need it.
217
218 if (dim == 3)
219 {
220 for (unsigned int count = 0; count < n_geometric_objects; count++)
221 {
222 unsigned int type;
223 in >> type;
224
225 if (cell_types[count] == 10 || cell_types[count] == 12)
226 {
227 if (cell_types[count] == 10)
228 is_tria_or_tet_mesh = true;
229 if (cell_types[count] == 12)
230 is_quad_or_hex_mesh = true;
231
232 // we assume that the file contains first all cells,
233 // and only then any faces or lines
234 AssertThrow(subcelldata.boundary_quads.size() == 0 &&
235 subcelldata.boundary_lines.size() == 0,
236 ExcNotImplemented());
237
238 cells.emplace_back(type);
239
240 for (unsigned int j = 0; j < type; j++) // loop to feed data
241 in >> cells.back().vertices[j];
242
243 cells.back().material_id = 0;
244 }
245
246 else if (cell_types[count] == 5 || cell_types[count] == 9)
247 {
248 if (cell_types[count] == 5)
249 is_tria_or_tet_mesh = true;
250 if (cell_types[count] == 9)
251 is_quad_or_hex_mesh = true;
252
253 // we assume that the file contains first all cells,
254 // then all faces, and finally all lines
255 AssertThrow(subcelldata.boundary_lines.size() == 0,
256 ExcNotImplemented());
257
258 subcelldata.boundary_quads.emplace_back(type);
259
260 for (unsigned int j = 0; j < type;
261 j++) // loop to feed the data to the boundary
262 in >> subcelldata.boundary_quads.back().vertices[j];
263
264 subcelldata.boundary_quads.back().material_id = 0;
265 }
266 else if (cell_types[count] == 3)
267 {
268 subcelldata.boundary_lines.emplace_back(type);
269
270 for (unsigned int j = 0; j < type;
271 j++) // loop to feed the data to the boundary
272 in >> subcelldata.boundary_lines.back().vertices[j];
273
274 subcelldata.boundary_lines.back().material_id = 0;
275 }
276
277 else
278 AssertThrow(
279 false,
280 ExcMessage(
281 "While reading VTK file, unknown file type encountered"));
282 }
283 }
284 else if (dim == 2)
285 {
286 for (unsigned int count = 0; count < n_geometric_objects; count++)
287 {
288 unsigned int type;
289 in >> type;
290
291 if (cell_types[count] == 5 || cell_types[count] == 9)
292 {
293 // we assume that the file contains first all cells,
294 // and only then any faces
295 AssertThrow(subcelldata.boundary_lines.size() == 0,
296 ExcNotImplemented());
297
298 if (cell_types[count] == 5)
299 is_tria_or_tet_mesh = true;
300 if (cell_types[count] == 9)
301 is_quad_or_hex_mesh = true;
302
303 cells.emplace_back(type);
304
305 for (unsigned int j = 0; j < type; j++) // loop to feed data
306 in >> cells.back().vertices[j];
307
308 cells.back().material_id = 0;
309 }
310
311 else if (cell_types[count] == 3)
312 {
313 // If this is encountered, the pointer comes out of the loop
314 // and starts processing boundaries.
315 subcelldata.boundary_lines.emplace_back(type);
316
317 for (unsigned int j = 0; j < type;
318 j++) // loop to feed the data to the boundary
319 {
320 in >> subcelldata.boundary_lines.back().vertices[j];
321 }
322
323 subcelldata.boundary_lines.back().material_id = 0;
324 }
325
326 else
327 AssertThrow(
328 false,
329 ExcMessage(
330 "While reading VTK file, unknown cell type encountered"));
331 }
332 }
333 else if (dim == 1)
334 {
335 for (unsigned int count = 0; count < n_geometric_objects; count++)
336 {
337 unsigned int type;
338 in >> type;
339
340 AssertThrow(
341 cell_types[count] == 3 && type == 2,
342 ExcMessage(
343 "While reading VTK file, unknown cell type encountered"));
344 cells.emplace_back(type);
345
346 for (unsigned int j = 0; j < type; j++) // loop to feed data
347 in >> cells.back().vertices[j];
348
349 cells.back().material_id = 0;
350 }
351 }
352 else
353 AssertThrow(false,
354 ExcMessage(
355 "While reading VTK file, failed to find CELLS section"));
356
357 /////////////////////Processing the CELL_TYPES
358 /// section////////////////////////
359
360 in >> keyword;
361
362 AssertThrow(
363 keyword == "CELL_TYPES",
364 ExcMessage(std::string(
365 "While reading VTK file, missing CELL_TYPES section. Found <" +
366 keyword + "> instead.")));
367
368 in >> n_ints;
369 AssertThrow(
370 n_ints == n_geometric_objects,
371 ExcMessage("The VTK reader found a CELL_DATA statement "
372 "that lists a total of " +
373 Utilities::int_to_string(n_ints) +
374 " cell data objects, but this needs to "
375 "equal the number of cells (which is " +
376 Utilities::int_to_string(cells.size()) +
377 ") plus the number of quads (" +
378 Utilities::int_to_string(subcelldata.boundary_quads.size()) +
379 " in 3d or the number of lines (" +
380 Utilities::int_to_string(subcelldata.boundary_lines.size()) +
381 ") in 2d."));
382
383 int tmp_int;
384 for (unsigned int i = 0; i < n_ints; ++i)
385 in >> tmp_int;
386
387 // Ignore everything up to CELL_DATA
388 while (in >> keyword)
389 if (keyword == "CELL_DATA")
390 {
391 unsigned int n_ids;
392 in >> n_ids;
393
394 AssertThrow(n_ids == n_geometric_objects,
395 ExcMessage("The VTK reader found a CELL_DATA statement "
396 "that lists a total of " +
397 Utilities::int_to_string(n_ids) +
398 " cell data objects, but this needs to "
399 "equal the number of cells (which is " +
400 Utilities::int_to_string(cells.size()) +
401 ") plus the number of quads (" +
402 Utilities::int_to_string(
403 subcelldata.boundary_quads.size()) +
404 " in 3d or the number of lines (" +
405 Utilities::int_to_string(
406 subcelldata.boundary_lines.size()) +
407 ") in 2d."));
408
409 const std::vector<std::string> data_sets{"MaterialID",
410 "ManifoldID"};
411
412 for (unsigned int i = 0; i < data_sets.size(); ++i)
413 {
414 // Ignore everything until we get to a SCALARS data set
415 while (in >> keyword)
416 if (keyword == "SCALARS")
417 {
418 // Now see if we know about this type of data set,
419 // if not, just ignore everything till the next SCALARS
420 // keyword
421 std::string field_name;
422 in >> field_name;
423 if (std::find(data_sets.begin(),
424 data_sets.end(),
425 field_name) == data_sets.end())
426 // The data set here is not one of the ones we know, so
427 // keep ignoring everything until the next SCALARS
428 // keyword.
429 continue;
430
431 // Now we got somewhere. Proceed from here, assert
432 // that the type of the table is int, and ignore the
433 // rest of the line.
434 // SCALARS MaterialID int 1
435 // (the last number is optional)
436 std::string line;
437 std::getline(in, line);
438 AssertThrow(
439 line.substr(1,
440 std::min(static_cast<std::size_t>(3),
441 line.size() - 1)) == "int",
442 ExcMessage(
443 "While reading VTK file, material- and manifold IDs can only have type 'int'."));
444
445 in >> keyword;
446 AssertThrow(
447 keyword == "LOOKUP_TABLE",
448 ExcMessage(
449 "While reading VTK file, missing keyword 'LOOKUP_TABLE'."));
450
451 in >> keyword;
452 AssertThrow(
453 keyword == "default",
454 ExcMessage(
455 "While reading VTK file, missing keyword 'default'."));
456
457 // read material or manifold ids first for all cells,
458 // then for all faces, and finally for all lines. the
459 // assumption that cells come before all faces and
460 // lines has been verified above via an assertion, so
461 // the order used in the following blocks makes sense
462 for (unsigned int i = 0; i < cells.size(); i++)
463 {
464 int id;
465 in >> id;
466 if (field_name == "MaterialID")
467 cells[i].material_id =
468 static_cast<types::material_id>(id);
469 else if (field_name == "ManifoldID")
470 cells[i].manifold_id =
471 static_cast<types::manifold_id>(id);
472 else
473 Assert(false, ExcInternalError());
474 }
475
476 if (dim == 3)
477 {
478 for (auto &boundary_quad : subcelldata.boundary_quads)
479 {
480 int id;
481 in >> id;
482 if (field_name == "MaterialID")
483 boundary_quad.material_id =
484 static_cast<types::material_id>(id);
485 else if (field_name == "ManifoldID")
486 boundary_quad.manifold_id =
487 static_cast<types::manifold_id>(id);
488 else
489 Assert(false, ExcInternalError());
490 }
491 for (auto &boundary_line : subcelldata.boundary_lines)
492 {
493 int id;
494 in >> id;
495 if (field_name == "MaterialID")
496 boundary_line.material_id =
497 static_cast<types::material_id>(id);
498 else if (field_name == "ManifoldID")
499 boundary_line.manifold_id =
500 static_cast<types::manifold_id>(id);
501 else
502 Assert(false, ExcInternalError());
503 }
504 }
505 else if (dim == 2)
506 {
507 for (auto &boundary_line : subcelldata.boundary_lines)
508 {
509 int id;
510 in >> id;
511 if (field_name == "MaterialID")
512 boundary_line.material_id =
513 static_cast<types::material_id>(id);
514 else if (field_name == "ManifoldID")
515 boundary_line.manifold_id =
516 static_cast<types::manifold_id>(id);
517 else
518 Assert(false, ExcInternalError());
519 }
520 }
521 }
522 }
523 }
524
525 Assert(subcelldata.check_consistency(dim), ExcInternalError());
526
527
528 // make sure that only either simplex or hypercube cells are available
529 //
530 // TODO: the functions below (GridTools::delete_unused_vertices(),
531 // GridTools::invert_all_cells_of_negative_grid(),
532 // GridReordering::reorder_cells(),
533 // Triangulation::create_triangulation_compatibility()) need to be
534 // revisited for simplex meshes
535 AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
536 ExcNotImplemented());
537
538 if (dim == 1 || is_quad_or_hex_mesh)
539 {
540 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
541
542 if (dim == spacedim)
543 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(
544 vertices, cells);
545
546 GridReordering<dim, spacedim>::reorder_cells(cells);
547 tria->create_triangulation_compatibility(vertices,
548 cells,
549 subcelldata);
550
551 return;
552 }
553 else
554 {
555 tria->create_triangulation(vertices, cells, subcelldata);
556 }
557 }
558 else
559 AssertThrow(false,
560 ExcMessage(
561 "While reading VTK file, failed to find CELLS section"));
562 }
563
564
565
566 template <int dim, int spacedim>
567 void
read_vtu(std::istream & in)568 GridIn<dim, spacedim>::read_vtu(std::istream &in)
569 {
570 namespace pt = boost::property_tree;
571 pt::ptree tree;
572 pt::read_xml(in, tree);
573 auto section = tree.get_optional<std::string>("VTKFile.dealiiData");
574
575 AssertThrow(section,
576 ExcMessage(
577 "While reading a VTU file, failed to find dealiiData section. "
578 "Notice that we can only read grid files in .vtu format that "
579 "were created by the deal.II library, using a call to "
580 "GridOut::write_vtu(), where the flag "
581 "GridOutFlags::Vtu::serialize_triangulation is set to true."));
582
583 const auto decoded =
584 Utilities::decode_base64({section->begin(), section->end() - 1});
585 const auto string_archive =
586 Utilities::decompress({decoded.begin(), decoded.end()});
587 std::istringstream in_stream(string_archive);
588 boost::archive::binary_iarchive ia(in_stream);
589 tria->load(ia, 0);
590 }
591
592
593 template <int dim, int spacedim>
594 void
read_unv(std::istream & in)595 GridIn<dim, spacedim>::read_unv(std::istream &in)
596 {
597 Assert(tria != nullptr, ExcNoTriangulationSelected());
598 Assert((dim == 2) || (dim == 3), ExcNotImplemented());
599
600 AssertThrow(in, ExcIO());
601 skip_comment_lines(in, '#'); // skip comments (if any) at beginning of file
602
603 int tmp;
604
605 AssertThrow(in, ExcIO());
606 in >> tmp;
607 AssertThrow(in, ExcIO());
608 in >> tmp;
609
610 // section 2411 describes vertices: see
611 // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2411
612 AssertThrow(tmp == 2411, ExcUnknownSectionType(tmp));
613
614 std::vector<Point<spacedim>> vertices; // vector of vertex coordinates
615 std::map<int, int>
616 vertex_indices; // # vert in unv (key) ---> # vert in deal.II (value)
617
618 int no_vertex = 0; // deal.II
619
620 while (tmp != -1) // we do until reach end of 2411
621 {
622 int no; // unv
623 int dummy;
624 double x[3];
625
626 AssertThrow(in, ExcIO());
627 in >> no;
628
629 tmp = no;
630 if (tmp == -1)
631 break;
632
633 in >> dummy >> dummy >> dummy;
634
635 AssertThrow(in, ExcIO());
636 in >> x[0] >> x[1] >> x[2];
637
638 vertices.emplace_back();
639
640 for (unsigned int d = 0; d < spacedim; d++)
641 vertices.back()(d) = x[d];
642
643 vertex_indices[no] = no_vertex;
644
645 no_vertex++;
646 }
647
648 AssertThrow(in, ExcIO());
649 in >> tmp;
650 AssertThrow(in, ExcIO());
651 in >> tmp;
652
653 // section 2412 describes elements: see
654 // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2412
655 AssertThrow(tmp == 2412, ExcUnknownSectionType(tmp));
656
657 std::vector<CellData<dim>> cells; // vector of cells
658 SubCellData subcelldata;
659
660 std::map<int, int>
661 cell_indices; // # cell in unv (key) ---> # cell in deal.II (value)
662 std::map<int, int>
663 line_indices; // # line in unv (key) ---> # line in deal.II (value)
664 std::map<int, int>
665 quad_indices; // # quad in unv (key) ---> # quad in deal.II (value)
666
667 int no_cell = 0; // deal.II
668 int no_line = 0; // deal.II
669 int no_quad = 0; // deal.II
670
671 while (tmp != -1) // we do until reach end of 2412
672 {
673 int no; // unv
674 int type;
675 int dummy;
676
677 AssertThrow(in, ExcIO());
678 in >> no;
679
680 tmp = no;
681 if (tmp == -1)
682 break;
683
684 in >> type >> dummy >> dummy >> dummy >> dummy;
685
686 AssertThrow((type == 11) || (type == 44) || (type == 94) || (type == 115),
687 ExcUnknownElementType(type));
688
689 if ((((type == 44) || (type == 94)) && (dim == 2)) ||
690 ((type == 115) && (dim == 3))) // cell
691 {
692 cells.emplace_back();
693
694 AssertThrow(in, ExcIO());
695 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
696 in >> cells.back().vertices[v];
697
698 cells.back().material_id = 0;
699
700 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
701 cells.back().vertices[v] = vertex_indices[cells.back().vertices[v]];
702
703 cell_indices[no] = no_cell;
704
705 no_cell++;
706 }
707 else if (((type == 11) && (dim == 2)) ||
708 ((type == 11) && (dim == 3))) // boundary line
709 {
710 AssertThrow(in, ExcIO());
711 in >> dummy >> dummy >> dummy;
712
713 subcelldata.boundary_lines.emplace_back();
714
715 AssertThrow(in, ExcIO());
716 for (unsigned int &vertex :
717 subcelldata.boundary_lines.back().vertices)
718 in >> vertex;
719
720 subcelldata.boundary_lines.back().material_id = 0;
721
722 for (unsigned int &vertex :
723 subcelldata.boundary_lines.back().vertices)
724 vertex = vertex_indices[vertex];
725
726 line_indices[no] = no_line;
727
728 no_line++;
729 }
730 else if (((type == 44) || (type == 94)) && (dim == 3)) // boundary quad
731 {
732 subcelldata.boundary_quads.emplace_back();
733
734 AssertThrow(in, ExcIO());
735 for (unsigned int &vertex :
736 subcelldata.boundary_quads.back().vertices)
737 in >> vertex;
738
739 subcelldata.boundary_quads.back().material_id = 0;
740
741 for (unsigned int &vertex :
742 subcelldata.boundary_quads.back().vertices)
743 vertex = vertex_indices[vertex];
744
745 quad_indices[no] = no_quad;
746
747 no_quad++;
748 }
749 else
750 AssertThrow(false,
751 ExcMessage("Unknown element label <" +
752 Utilities::int_to_string(type) +
753 "> when running in dim=" +
754 Utilities::int_to_string(dim)));
755 }
756
757 // note that so far all materials and bcs are explicitly set to 0
758 // if we do not need more info on materials and bcs - this is end of file
759 // if we do - section 2467 or 2477 comes
760
761 in >> tmp; // tmp can be either -1 or end-of-file
762
763 if (!in.eof())
764 {
765 AssertThrow(in, ExcIO());
766 in >> tmp;
767
768 // section 2467 (2477) describes (materials - first and bcs - second) or
769 // (bcs - first and materials - second) - sequence depends on which
770 // group is created first: see
771 // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2467
772 AssertThrow((tmp == 2467) || (tmp == 2477), ExcUnknownSectionType(tmp));
773
774 while (tmp != -1) // we do until reach end of 2467 or 2477
775 {
776 int n_entities; // number of entities in group
777 int id; // id is either material or bc
778 int no; // unv
779 int dummy;
780
781 AssertThrow(in, ExcIO());
782 in >> dummy;
783
784 tmp = dummy;
785 if (tmp == -1)
786 break;
787
788 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
789 n_entities;
790
791 AssertThrow(in, ExcIO());
792 in >> id;
793
794 const unsigned int n_lines =
795 (n_entities % 2 == 0) ? (n_entities / 2) : ((n_entities + 1) / 2);
796
797 for (unsigned int line = 0; line < n_lines; line++)
798 {
799 unsigned int n_fragments;
800
801 if (line == n_lines - 1)
802 n_fragments = (n_entities % 2 == 0) ? (2) : (1);
803 else
804 n_fragments = 2;
805
806 for (unsigned int no_fragment = 0; no_fragment < n_fragments;
807 no_fragment++)
808 {
809 AssertThrow(in, ExcIO());
810 in >> dummy >> no >> dummy >> dummy;
811
812 if (cell_indices.count(no) > 0) // cell - material
813 cells[cell_indices[no]].material_id = id;
814
815 if (line_indices.count(no) > 0) // boundary line - bc
816 subcelldata.boundary_lines[line_indices[no]].material_id =
817 id;
818
819 if (quad_indices.count(no) > 0) // boundary quad - bc
820 subcelldata.boundary_quads[quad_indices[no]].material_id =
821 id;
822 }
823 }
824 }
825 }
826
827 Assert(subcelldata.check_consistency(dim), ExcInternalError());
828
829 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
830
831 if (dim == spacedim)
832 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
833 cells);
834
835 GridReordering<dim, spacedim>::reorder_cells(cells);
836
837 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
838 }
839
840
841
842 template <int dim, int spacedim>
843 void
read_ucd(std::istream & in,const bool apply_all_indicators_to_manifolds)844 GridIn<dim, spacedim>::read_ucd(std::istream &in,
845 const bool apply_all_indicators_to_manifolds)
846 {
847 Assert(tria != nullptr, ExcNoTriangulationSelected());
848 AssertThrow(in, ExcIO());
849
850 // skip comments at start of file
851 skip_comment_lines(in, '#');
852
853
854 unsigned int n_vertices;
855 unsigned int n_cells;
856 int dummy;
857
858 in >> n_vertices >> n_cells >> dummy // number of data vectors
859 >> dummy // cell data
860 >> dummy; // model data
861 AssertThrow(in, ExcIO());
862
863 // set up array of vertices
864 std::vector<Point<spacedim>> vertices(n_vertices);
865 // set up mapping between numbering
866 // in ucd-file (key) and in the
867 // vertices vector
868 std::map<int, int> vertex_indices;
869
870 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
871 {
872 int vertex_number;
873 double x[3];
874
875 // read vertex
876 AssertThrow(in, ExcIO());
877 in >> vertex_number >> x[0] >> x[1] >> x[2];
878
879 // store vertex
880 for (unsigned int d = 0; d < spacedim; ++d)
881 vertices[vertex](d) = x[d];
882 // store mapping; note that
883 // vertices_indices[i] is automatically
884 // created upon first usage
885 vertex_indices[vertex_number] = vertex;
886 }
887
888 // set up array of cells
889 std::vector<CellData<dim>> cells;
890 SubCellData subcelldata;
891
892 for (unsigned int cell = 0; cell < n_cells; ++cell)
893 {
894 // note that since in the input
895 // file we found the number of
896 // cells at the top, there
897 // should still be input here,
898 // so check this:
899 AssertThrow(in, ExcIO());
900
901 std::string cell_type;
902
903 // we use an unsigned int because we
904 // fill this variable through an read-in process
905 unsigned int material_id;
906
907 in >> dummy // cell number
908 >> material_id;
909 in >> cell_type;
910
911 if (((cell_type == "line") && (dim == 1)) ||
912 ((cell_type == "quad") && (dim == 2)) ||
913 ((cell_type == "hex") && (dim == 3)))
914 // found a cell
915 {
916 // allocate and read indices
917 cells.emplace_back();
918 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
919 in >> cells.back().vertices[i];
920
921 // to make sure that the cast won't fail
922 Assert(material_id <= std::numeric_limits<types::material_id>::max(),
923 ExcIndexRange(material_id,
924 0,
925 std::numeric_limits<types::material_id>::max()));
926 // we use only material_ids in the range from 0 to
927 // numbers::invalid_material_id-1
928 AssertIndexRange(material_id, numbers::invalid_material_id);
929
930 if (apply_all_indicators_to_manifolds)
931 cells.back().manifold_id =
932 static_cast<types::manifold_id>(material_id);
933 cells.back().material_id = material_id;
934
935 // transform from ucd to
936 // consecutive numbering
937 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
938 if (vertex_indices.find(cells.back().vertices[i]) !=
939 vertex_indices.end())
940 // vertex with this index exists
941 cells.back().vertices[i] =
942 vertex_indices[cells.back().vertices[i]];
943 else
944 {
945 // no such vertex index
946 AssertThrow(false,
947 ExcInvalidVertexIndex(cell,
948 cells.back().vertices[i]));
949
950 cells.back().vertices[i] = numbers::invalid_unsigned_int;
951 }
952 }
953 else if ((cell_type == "line") && ((dim == 2) || (dim == 3)))
954 // boundary info
955 {
956 subcelldata.boundary_lines.emplace_back();
957 in >> subcelldata.boundary_lines.back().vertices[0] >>
958 subcelldata.boundary_lines.back().vertices[1];
959
960 // to make sure that the cast won't fail
961 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
962 ExcIndexRange(material_id,
963 0,
964 std::numeric_limits<types::boundary_id>::max()));
965 // we use only boundary_ids in the range from 0 to
966 // numbers::internal_face_boundary_id-1
967 AssertIndexRange(material_id, numbers::internal_face_boundary_id);
968
969 // Make sure to set both manifold id and boundary id appropriately in
970 // both cases:
971 // numbers::internal_face_boundary_id and numbers::flat_manifold_id
972 // are ignored in Triangulation::create_triangulation.
973 if (apply_all_indicators_to_manifolds)
974 {
975 subcelldata.boundary_lines.back().boundary_id =
976 numbers::internal_face_boundary_id;
977 subcelldata.boundary_lines.back().manifold_id =
978 static_cast<types::manifold_id>(material_id);
979 }
980 else
981 {
982 subcelldata.boundary_lines.back().boundary_id =
983 static_cast<types::boundary_id>(material_id);
984 subcelldata.boundary_lines.back().manifold_id =
985 numbers::flat_manifold_id;
986 }
987
988 // transform from ucd to
989 // consecutive numbering
990 for (unsigned int &vertex :
991 subcelldata.boundary_lines.back().vertices)
992 if (vertex_indices.find(vertex) != vertex_indices.end())
993 // vertex with this index exists
994 vertex = vertex_indices[vertex];
995 else
996 {
997 // no such vertex index
998 AssertThrow(false, ExcInvalidVertexIndex(cell, vertex));
999 vertex = numbers::invalid_unsigned_int;
1000 }
1001 }
1002 else if ((cell_type == "quad") && (dim == 3))
1003 // boundary info
1004 {
1005 subcelldata.boundary_quads.emplace_back();
1006 in >> subcelldata.boundary_quads.back().vertices[0] >>
1007 subcelldata.boundary_quads.back().vertices[1] >>
1008 subcelldata.boundary_quads.back().vertices[2] >>
1009 subcelldata.boundary_quads.back().vertices[3];
1010
1011 // to make sure that the cast won't fail
1012 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
1013 ExcIndexRange(material_id,
1014 0,
1015 std::numeric_limits<types::boundary_id>::max()));
1016 // we use only boundary_ids in the range from 0 to
1017 // numbers::internal_face_boundary_id-1
1018 AssertIndexRange(material_id, numbers::internal_face_boundary_id);
1019
1020 // Make sure to set both manifold id and boundary id appropriately in
1021 // both cases:
1022 // numbers::internal_face_boundary_id and numbers::flat_manifold_id
1023 // are ignored in Triangulation::create_triangulation.
1024 if (apply_all_indicators_to_manifolds)
1025 {
1026 subcelldata.boundary_quads.back().boundary_id =
1027 numbers::internal_face_boundary_id;
1028 subcelldata.boundary_quads.back().manifold_id =
1029 static_cast<types::manifold_id>(material_id);
1030 }
1031 else
1032 {
1033 subcelldata.boundary_quads.back().boundary_id =
1034 static_cast<types::boundary_id>(material_id);
1035 subcelldata.boundary_quads.back().manifold_id =
1036 numbers::flat_manifold_id;
1037 }
1038
1039 // transform from ucd to
1040 // consecutive numbering
1041 for (unsigned int &vertex :
1042 subcelldata.boundary_quads.back().vertices)
1043 if (vertex_indices.find(vertex) != vertex_indices.end())
1044 // vertex with this index exists
1045 vertex = vertex_indices[vertex];
1046 else
1047 {
1048 // no such vertex index
1049 Assert(false, ExcInvalidVertexIndex(cell, vertex));
1050 vertex = numbers::invalid_unsigned_int;
1051 }
1052 }
1053 else
1054 // cannot read this
1055 AssertThrow(false, ExcUnknownIdentifier(cell_type));
1056 }
1057
1058
1059 // check that no forbidden arrays are used
1060 Assert(subcelldata.check_consistency(dim), ExcInternalError());
1061
1062 AssertThrow(in, ExcIO());
1063
1064 // do some clean-up on vertices...
1065 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1066 // ... and cells
1067 if (dim == spacedim)
1068 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
1069 cells);
1070 GridReordering<dim, spacedim>::reorder_cells(cells);
1071 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1072 }
1073
1074 namespace
1075 {
1076 template <int dim, int spacedim>
1077 class Abaqus_to_UCD
1078 {
1079 public:
1080 Abaqus_to_UCD();
1081
1082 void
1083 read_in_abaqus(std::istream &in);
1084 void
1085 write_out_avs_ucd(std::ostream &out) const;
1086
1087 private:
1088 const double tolerance;
1089
1090 std::vector<double>
1091 get_global_node_numbers(const int face_cell_no,
1092 const int face_cell_face_no) const;
1093
1094 // NL: Stored as [ global node-id (int), x-coord, y-coord, z-coord ]
1095 std::vector<std::vector<double>> node_list;
1096 // CL: Stored as [ material-id (int), node1, node2, node3, node4, node5,
1097 // node6, node7, node8 ]
1098 std::vector<std::vector<double>> cell_list;
1099 // FL: Stored as [ sideset-id (int), node1, node2, node3, node4 ]
1100 std::vector<std::vector<double>> face_list;
1101 // ELSET: Stored as [ (std::string) elset_name = (std::vector) of cells
1102 // numbers]
1103 std::map<std::string, std::vector<int>> elsets_list;
1104 };
1105 } // namespace
1106
1107 template <int dim, int spacedim>
1108 void
read_abaqus(std::istream & in,const bool apply_all_indicators_to_manifolds)1109 GridIn<dim, spacedim>::read_abaqus(std::istream &in,
1110 const bool apply_all_indicators_to_manifolds)
1111 {
1112 Assert(tria != nullptr, ExcNoTriangulationSelected());
1113 // This implementation has only been verified for:
1114 // - 2d grids with codimension 0
1115 // - 3d grids with codimension 0
1116 // - 3d grids with codimension 1
1117 Assert((spacedim == 2 && dim == spacedim) ||
1118 (spacedim == 3 && (dim == spacedim || dim == spacedim - 1)),
1119 ExcNotImplemented());
1120 AssertThrow(in, ExcIO());
1121
1122 // Read in the Abaqus file into an intermediate object
1123 // that is to be passed along to the UCD reader
1124 Abaqus_to_UCD<dim, spacedim> abaqus_to_ucd;
1125 abaqus_to_ucd.read_in_abaqus(in);
1126
1127 std::stringstream in_ucd;
1128 abaqus_to_ucd.write_out_avs_ucd(in_ucd);
1129
1130 // This next call is wrapped in a try-catch for the following reason:
1131 // It ensures that if the Abaqus mesh is read in correctly but produces
1132 // an erroneous result then the user is alerted to the source of the problem
1133 // and doesn't think that they've somehow called the wrong function.
1134 try
1135 {
1136 read_ucd(in_ucd, apply_all_indicators_to_manifolds);
1137 }
1138 catch (std::exception &exc)
1139 {
1140 std::cerr << "Exception on processing internal UCD data: " << std::endl
1141 << exc.what() << std::endl;
1142
1143 AssertThrow(
1144 false,
1145 ExcMessage(
1146 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1147 "More information is provided in an error message printed above. "
1148 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1149 "listed in the documentation?"));
1150 }
1151 catch (...)
1152 {
1153 AssertThrow(
1154 false,
1155 ExcMessage(
1156 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1157 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1158 "listed in the documentation?"));
1159 }
1160 }
1161
1162
1163 template <int dim, int spacedim>
1164 void
read_dbmesh(std::istream & in)1165 GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
1166 {
1167 Assert(tria != nullptr, ExcNoTriangulationSelected());
1168 Assert(dim == 2, ExcNotImplemented());
1169
1170 AssertThrow(in, ExcIO());
1171
1172 // skip comments at start of file
1173 skip_comment_lines(in, '#');
1174
1175 // first read in identifier string
1176 std::string line;
1177 getline(in, line);
1178
1179 AssertThrow(line == "MeshVersionFormatted 0", ExcInvalidDBMESHInput(line));
1180
1181 skip_empty_lines(in);
1182
1183 // next read dimension
1184 getline(in, line);
1185 AssertThrow(line == "Dimension", ExcInvalidDBMESHInput(line));
1186 unsigned int dimension;
1187 in >> dimension;
1188 AssertThrow(dimension == dim, ExcDBMESHWrongDimension(dimension));
1189 skip_empty_lines(in);
1190
1191 // now there are a lot of fields of
1192 // which we don't know the exact
1193 // meaning and which are far from
1194 // being properly documented in the
1195 // manual. we skip everything until
1196 // we find a comment line with the
1197 // string "# END". at some point in
1198 // the future, someone may have the
1199 // knowledge to parse and interpret
1200 // the other fields in between as
1201 // well...
1202 while (getline(in, line), line.find("# END") == std::string::npos)
1203 ;
1204 skip_empty_lines(in);
1205
1206
1207 // now read vertices
1208 getline(in, line);
1209 AssertThrow(line == "Vertices", ExcInvalidDBMESHInput(line));
1210
1211 unsigned int n_vertices;
1212 double dummy;
1213
1214 in >> n_vertices;
1215 std::vector<Point<spacedim>> vertices(n_vertices);
1216 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1217 {
1218 // read vertex coordinates
1219 for (unsigned int d = 0; d < dim; ++d)
1220 in >> vertices[vertex][d];
1221 // read Ref phi_i, whatever that may be
1222 in >> dummy;
1223 }
1224 AssertThrow(in, ExcInvalidDBMeshFormat());
1225
1226 skip_empty_lines(in);
1227
1228 // read edges. we ignore them at
1229 // present, so just read them and
1230 // discard the input
1231 getline(in, line);
1232 AssertThrow(line == "Edges", ExcInvalidDBMESHInput(line));
1233
1234 unsigned int n_edges;
1235 in >> n_edges;
1236 for (unsigned int edge = 0; edge < n_edges; ++edge)
1237 {
1238 // read vertex indices
1239 in >> dummy >> dummy;
1240 // read Ref phi_i, whatever that may be
1241 in >> dummy;
1242 }
1243 AssertThrow(in, ExcInvalidDBMeshFormat());
1244
1245 skip_empty_lines(in);
1246
1247
1248
1249 // read cracked edges (whatever
1250 // that may be). we ignore them at
1251 // present, so just read them and
1252 // discard the input
1253 getline(in, line);
1254 AssertThrow(line == "CrackedEdges", ExcInvalidDBMESHInput(line));
1255
1256 in >> n_edges;
1257 for (unsigned int edge = 0; edge < n_edges; ++edge)
1258 {
1259 // read vertex indices
1260 in >> dummy >> dummy;
1261 // read Ref phi_i, whatever that may be
1262 in >> dummy;
1263 }
1264 AssertThrow(in, ExcInvalidDBMeshFormat());
1265
1266 skip_empty_lines(in);
1267
1268
1269 // now read cells.
1270 // set up array of cells
1271 getline(in, line);
1272 AssertThrow(line == "Quadrilaterals", ExcInvalidDBMESHInput(line));
1273
1274 std::vector<CellData<dim>> cells;
1275 SubCellData subcelldata;
1276 unsigned int n_cells;
1277 in >> n_cells;
1278 for (unsigned int cell = 0; cell < n_cells; ++cell)
1279 {
1280 // read in vertex numbers. they
1281 // are 1-based, so subtract one
1282 cells.emplace_back();
1283 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1284 {
1285 in >> cells.back().vertices[i];
1286
1287 AssertThrow((cells.back().vertices[i] >= 1) &&
1288 (static_cast<unsigned int>(cells.back().vertices[i]) <=
1289 vertices.size()),
1290 ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
1291
1292 --cells.back().vertices[i];
1293 }
1294
1295 // read and discard Ref phi_i
1296 in >> dummy;
1297 }
1298 AssertThrow(in, ExcInvalidDBMeshFormat());
1299
1300 skip_empty_lines(in);
1301
1302
1303 // then there are again a whole lot
1304 // of fields of which I have no
1305 // clue what they mean. skip them
1306 // all and leave the interpretation
1307 // to other implementors...
1308 while (getline(in, line), ((line.find("End") == std::string::npos) && (in)))
1309 ;
1310 // ok, so we are not at the end of
1311 // the file, that's it, mostly
1312
1313
1314 // check that no forbidden arrays are used
1315 Assert(subcelldata.check_consistency(dim), ExcInternalError());
1316
1317 AssertThrow(in, ExcIO());
1318
1319 // do some clean-up on vertices...
1320 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1321 // ...and cells
1322 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
1323 cells);
1324 GridReordering<dim, spacedim>::reorder_cells(cells);
1325 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1326 }
1327
1328
1329
1330 template <int dim, int spacedim>
1331 void
read_xda(std::istream &)1332 GridIn<dim, spacedim>::read_xda(std::istream &)
1333 {
1334 Assert(false, ExcNotImplemented());
1335 }
1336
1337
1338
1339 template <>
1340 void
read_xda(std::istream & in)1341 GridIn<2>::read_xda(std::istream &in)
1342 {
1343 Assert(tria != nullptr, ExcNoTriangulationSelected());
1344 AssertThrow(in, ExcIO());
1345
1346 std::string line;
1347 // skip comments at start of file
1348 getline(in, line);
1349
1350
1351 unsigned int n_vertices;
1352 unsigned int n_cells;
1353
1354 // read cells, throw away rest of line
1355 in >> n_cells;
1356 getline(in, line);
1357
1358 in >> n_vertices;
1359 getline(in, line);
1360
1361 // ignore following 8 lines
1362 for (unsigned int i = 0; i < 8; ++i)
1363 getline(in, line);
1364
1365 // set up array of cells
1366 std::vector<CellData<2>> cells(n_cells);
1367 SubCellData subcelldata;
1368
1369 for (unsigned int cell = 0; cell < n_cells; ++cell)
1370 {
1371 // note that since in the input
1372 // file we found the number of
1373 // cells at the top, there
1374 // should still be input here,
1375 // so check this:
1376 AssertThrow(in, ExcIO());
1377 Assert(GeometryInfo<2>::vertices_per_cell == 4, ExcInternalError());
1378
1379 for (unsigned int &vertex : cells[cell].vertices)
1380 in >> vertex;
1381 }
1382
1383
1384
1385 // set up array of vertices
1386 std::vector<Point<2>> vertices(n_vertices);
1387 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1388 {
1389 double x[3];
1390
1391 // read vertex
1392 in >> x[0] >> x[1] >> x[2];
1393
1394 // store vertex
1395 for (unsigned int d = 0; d < 2; ++d)
1396 vertices[vertex](d) = x[d];
1397 }
1398 AssertThrow(in, ExcIO());
1399
1400 // do some clean-up on vertices...
1401 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1402 // ... and cells
1403 GridReordering<2>::invert_all_cells_of_negative_grid(vertices, cells);
1404 GridReordering<2>::reorder_cells(cells);
1405 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1406 }
1407
1408
1409
1410 template <>
1411 void
read_xda(std::istream & in)1412 GridIn<3>::read_xda(std::istream &in)
1413 {
1414 Assert(tria != nullptr, ExcNoTriangulationSelected());
1415 AssertThrow(in, ExcIO());
1416
1417 static const unsigned int xda_to_dealII_map[] = {0, 1, 5, 4, 3, 2, 6, 7};
1418
1419 std::string line;
1420 // skip comments at start of file
1421 getline(in, line);
1422
1423
1424 unsigned int n_vertices;
1425 unsigned int n_cells;
1426
1427 // read cells, throw away rest of line
1428 in >> n_cells;
1429 getline(in, line);
1430
1431 in >> n_vertices;
1432 getline(in, line);
1433
1434 // ignore following 8 lines
1435 for (unsigned int i = 0; i < 8; ++i)
1436 getline(in, line);
1437
1438 // set up array of cells
1439 std::vector<CellData<3>> cells(n_cells);
1440 SubCellData subcelldata;
1441
1442 for (unsigned int cell = 0; cell < n_cells; ++cell)
1443 {
1444 // note that since in the input
1445 // file we found the number of
1446 // cells at the top, there
1447 // should still be input here,
1448 // so check this:
1449 AssertThrow(in, ExcIO());
1450 Assert(GeometryInfo<3>::vertices_per_cell == 8, ExcInternalError());
1451
1452 unsigned int xda_ordered_nodes[8];
1453
1454 for (unsigned int &xda_ordered_node : xda_ordered_nodes)
1455 in >> xda_ordered_node;
1456
1457 for (unsigned int i = 0; i < 8; i++)
1458 cells[cell].vertices[i] = xda_ordered_nodes[xda_to_dealII_map[i]];
1459 }
1460
1461
1462
1463 // set up array of vertices
1464 std::vector<Point<3>> vertices(n_vertices);
1465 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1466 {
1467 double x[3];
1468
1469 // read vertex
1470 in >> x[0] >> x[1] >> x[2];
1471
1472 // store vertex
1473 for (unsigned int d = 0; d < 3; ++d)
1474 vertices[vertex](d) = x[d];
1475 }
1476 AssertThrow(in, ExcIO());
1477
1478 // do some clean-up on vertices...
1479 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
1480 // ... and cells
1481 GridReordering<3>::invert_all_cells_of_negative_grid(vertices, cells);
1482 GridReordering<3>::reorder_cells(cells);
1483 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
1484 }
1485
1486
1487
1488 template <int dim, int spacedim>
1489 void
read_msh(std::istream & in)1490 GridIn<dim, spacedim>::read_msh(std::istream &in)
1491 {
1492 Assert(tria != nullptr, ExcNoTriangulationSelected());
1493 AssertThrow(in, ExcIO());
1494
1495 unsigned int n_vertices;
1496 unsigned int n_cells;
1497 unsigned int dummy;
1498 std::string line;
1499 // This array stores maps from the 'entities' to the 'physical tags' for
1500 // points, curves, surfaces and volumes. We use this information later to
1501 // assign boundary ids.
1502 std::array<std::map<int, int>, 4> tag_maps;
1503
1504 in >> line;
1505
1506 // first determine file format
1507 unsigned int gmsh_file_format = 0;
1508 if (line == "$NOD")
1509 gmsh_file_format = 10;
1510 else if (line == "$MeshFormat")
1511 gmsh_file_format = 20;
1512 else
1513 AssertThrow(false, ExcInvalidGMSHInput(line));
1514
1515 // if file format is 2.0 or greater then we also have to read the rest of the
1516 // header
1517 if (gmsh_file_format == 20)
1518 {
1519 double version;
1520 unsigned int file_type, data_size;
1521
1522 in >> version >> file_type >> data_size;
1523
1524 Assert((version >= 2.0) && (version <= 4.1), ExcNotImplemented());
1525 gmsh_file_format = static_cast<unsigned int>(version * 10);
1526
1527 Assert(file_type == 0, ExcNotImplemented());
1528 Assert(data_size == sizeof(double), ExcNotImplemented());
1529
1530 // read the end of the header and the first line of the nodes description
1531 // to synch ourselves with the format 1 handling above
1532 in >> line;
1533 AssertThrow(line == "$EndMeshFormat", ExcInvalidGMSHInput(line));
1534
1535 in >> line;
1536 // if the next block is of kind $PhysicalNames, ignore it
1537 if (line == "$PhysicalNames")
1538 {
1539 do
1540 {
1541 in >> line;
1542 }
1543 while (line != "$EndPhysicalNames");
1544 in >> line;
1545 }
1546
1547 // if the next block is of kind $Entities, parse it
1548 if (line == "$Entities")
1549 {
1550 unsigned long n_points, n_curves, n_surfaces, n_volumes;
1551
1552 in >> n_points >> n_curves >> n_surfaces >> n_volumes;
1553 for (unsigned int i = 0; i < n_points; ++i)
1554 {
1555 // parse point ids
1556 int tag;
1557 unsigned int n_physicals;
1558 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1559 box_max_z;
1560
1561 // we only care for 'tag' as key for tag_maps[0]
1562 if (gmsh_file_format > 40)
1563 {
1564 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
1565 n_physicals;
1566 box_max_x = box_min_x;
1567 box_max_y = box_min_y;
1568 box_max_z = box_min_z;
1569 }
1570 else
1571 {
1572 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
1573 box_max_x >> box_max_y >> box_max_z >> n_physicals;
1574 }
1575 // if there is a physical tag, we will use it as boundary id below
1576 AssertThrow(n_physicals < 2,
1577 ExcMessage("More than one tag is not supported!"));
1578 // if there is no physical tag, use 0 as default
1579 int physical_tag = 0;
1580 for (unsigned int j = 0; j < n_physicals; ++j)
1581 in >> physical_tag;
1582 tag_maps[0][tag] = physical_tag;
1583 }
1584 for (unsigned int i = 0; i < n_curves; ++i)
1585 {
1586 // parse curve ids
1587 int tag;
1588 unsigned int n_physicals;
1589 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1590 box_max_z;
1591
1592 // we only care for 'tag' as key for tag_maps[1]
1593 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1594 box_max_y >> box_max_z >> n_physicals;
1595 // if there is a physical tag, we will use it as boundary id below
1596 AssertThrow(n_physicals < 2,
1597 ExcMessage("More than one tag is not supported!"));
1598 // if there is no physical tag, use 0 as default
1599 int physical_tag = 0;
1600 for (unsigned int j = 0; j < n_physicals; ++j)
1601 in >> physical_tag;
1602 tag_maps[1][tag] = physical_tag;
1603 // we don't care about the points associated to a curve, but have
1604 // to parse them anyway because their format is unstructured
1605 in >> n_points;
1606 for (unsigned int j = 0; j < n_points; ++j)
1607 in >> tag;
1608 }
1609
1610 for (unsigned int i = 0; i < n_surfaces; ++i)
1611 {
1612 // parse surface ids
1613 int tag;
1614 unsigned int n_physicals;
1615 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1616 box_max_z;
1617
1618 // we only care for 'tag' as key for tag_maps[2]
1619 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1620 box_max_y >> box_max_z >> n_physicals;
1621 // if there is a physical tag, we will use it as boundary id below
1622 AssertThrow(n_physicals < 2,
1623 ExcMessage("More than one tag is not supported!"));
1624 // if there is no physical tag, use 0 as default
1625 int physical_tag = 0;
1626 for (unsigned int j = 0; j < n_physicals; ++j)
1627 in >> physical_tag;
1628 tag_maps[2][tag] = physical_tag;
1629 // we don't care about the curves associated to a surface, but
1630 // have to parse them anyway because their format is unstructured
1631 in >> n_curves;
1632 for (unsigned int j = 0; j < n_curves; ++j)
1633 in >> tag;
1634 }
1635 for (unsigned int i = 0; i < n_volumes; ++i)
1636 {
1637 // parse volume ids
1638 int tag;
1639 unsigned int n_physicals;
1640 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
1641 box_max_z;
1642
1643 // we only care for 'tag' as key for tag_maps[3]
1644 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
1645 box_max_y >> box_max_z >> n_physicals;
1646 // if there is a physical tag, we will use it as boundary id below
1647 AssertThrow(n_physicals < 2,
1648 ExcMessage("More than one tag is not supported!"));
1649 // if there is no physical tag, use 0 as default
1650 int physical_tag = 0;
1651 for (unsigned int j = 0; j < n_physicals; ++j)
1652 in >> physical_tag;
1653 tag_maps[3][tag] = physical_tag;
1654 // we don't care about the surfaces associated to a volume, but
1655 // have to parse them anyway because their format is unstructured
1656 in >> n_surfaces;
1657 for (unsigned int j = 0; j < n_surfaces; ++j)
1658 in >> tag;
1659 }
1660 in >> line;
1661 AssertThrow(line == "$EndEntities", ExcInvalidGMSHInput(line));
1662 in >> line;
1663 }
1664
1665 // if the next block is of kind $PartitionedEntities, ignore it
1666 if (line == "$PartitionedEntities")
1667 {
1668 do
1669 {
1670 in >> line;
1671 }
1672 while (line != "$EndPartitionedEntities");
1673 in >> line;
1674 }
1675
1676 // but the next thing should,
1677 // in any case, be the list of
1678 // nodes:
1679 AssertThrow(line == "$Nodes", ExcInvalidGMSHInput(line));
1680 }
1681
1682 // now read the nodes list
1683 int n_entity_blocks = 1;
1684 if (gmsh_file_format > 40)
1685 {
1686 int min_node_tag;
1687 int max_node_tag;
1688 in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag;
1689 }
1690 else if (gmsh_file_format == 40)
1691 {
1692 in >> n_entity_blocks >> n_vertices;
1693 }
1694 else
1695 in >> n_vertices;
1696 std::vector<Point<spacedim>> vertices(n_vertices);
1697 // set up mapping between numbering
1698 // in msh-file (nod) and in the
1699 // vertices vector
1700 std::map<int, int> vertex_indices;
1701
1702 {
1703 unsigned int global_vertex = 0;
1704 for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
1705 {
1706 int parametric;
1707 unsigned long numNodes;
1708
1709 if (gmsh_file_format < 40)
1710 {
1711 numNodes = n_vertices;
1712 parametric = 0;
1713 }
1714 else
1715 {
1716 // for gmsh_file_format 4.1 the order of tag and dim is reversed,
1717 // but we are ignoring both anyway.
1718 int tagEntity, dimEntity;
1719 in >> tagEntity >> dimEntity >> parametric >> numNodes;
1720 }
1721
1722 std::vector<int> vertex_numbers;
1723 int vertex_number;
1724 if (gmsh_file_format > 40)
1725 for (unsigned long vertex_per_entity = 0;
1726 vertex_per_entity < numNodes;
1727 ++vertex_per_entity)
1728 {
1729 in >> vertex_number;
1730 vertex_numbers.push_back(vertex_number);
1731 }
1732
1733 for (unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes;
1734 ++vertex_per_entity, ++global_vertex)
1735 {
1736 int vertex_number;
1737 double x[3];
1738
1739 // read vertex
1740 if (gmsh_file_format > 40)
1741 {
1742 vertex_number = vertex_numbers[vertex_per_entity];
1743 in >> x[0] >> x[1] >> x[2];
1744 }
1745 else
1746 in >> vertex_number >> x[0] >> x[1] >> x[2];
1747
1748 for (unsigned int d = 0; d < spacedim; ++d)
1749 vertices[global_vertex](d) = x[d];
1750 // store mapping
1751 vertex_indices[vertex_number] = global_vertex;
1752
1753 // ignore parametric coordinates
1754 if (parametric != 0)
1755 {
1756 double u = 0.;
1757 double v = 0.;
1758 in >> u >> v;
1759 (void)u;
1760 (void)v;
1761 }
1762 }
1763 }
1764 AssertDimension(global_vertex, n_vertices);
1765 }
1766
1767 // Assert we reached the end of the block
1768 in >> line;
1769 static const std::string end_nodes_marker[] = {"$ENDNOD", "$EndNodes"};
1770 AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1],
1771 ExcInvalidGMSHInput(line));
1772
1773 // Now read in next bit
1774 in >> line;
1775 static const std::string begin_elements_marker[] = {"$ELM", "$Elements"};
1776 AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1],
1777 ExcInvalidGMSHInput(line));
1778
1779 // now read the cell list
1780 if (gmsh_file_format > 40)
1781 {
1782 int min_node_tag;
1783 int max_node_tag;
1784 in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag;
1785 }
1786 else if (gmsh_file_format == 40)
1787 {
1788 in >> n_entity_blocks >> n_cells;
1789 }
1790 else
1791 {
1792 n_entity_blocks = 1;
1793 in >> n_cells;
1794 }
1795
1796 // set up array of cells and subcells (faces). In 1d, there is currently no
1797 // standard way in deal.II to pass boundary indicators attached to individual
1798 // vertices, so do this by hand via the boundary_ids_1d array
1799 std::vector<CellData<dim>> cells;
1800 SubCellData subcelldata;
1801 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
1802 bool is_quad_or_hex_mesh = false;
1803 bool is_tria_or_tet_mesh = false;
1804
1805 {
1806 unsigned int global_cell = 0;
1807 for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
1808 {
1809 unsigned int material_id;
1810 unsigned long numElements;
1811 int cell_type;
1812
1813 if (gmsh_file_format < 40)
1814 {
1815 material_id = 0;
1816 cell_type = 0;
1817 numElements = n_cells;
1818 }
1819 else if (gmsh_file_format == 40)
1820 {
1821 int tagEntity, dimEntity;
1822 in >> tagEntity >> dimEntity >> cell_type >> numElements;
1823 material_id = tag_maps[dimEntity][tagEntity];
1824 }
1825 else
1826 {
1827 // for gmsh_file_format 4.1 the order of tag and dim is reversed,
1828 int tagEntity, dimEntity;
1829 in >> dimEntity >> tagEntity >> cell_type >> numElements;
1830 material_id = tag_maps[dimEntity][tagEntity];
1831 }
1832
1833 for (unsigned int cell_per_entity = 0; cell_per_entity < numElements;
1834 ++cell_per_entity, ++global_cell)
1835 {
1836 // note that since in the input
1837 // file we found the number of
1838 // cells at the top, there
1839 // should still be input here,
1840 // so check this:
1841 AssertThrow(in, ExcIO());
1842
1843 unsigned int nod_num;
1844
1845 /*
1846 For file format version 1, the format of each cell is as follows:
1847 elm-number elm-type reg-phys reg-elem number-of-nodes
1848 node-number-list
1849
1850 However, for version 2, the format reads like this:
1851 elm-number elm-type number-of-tags < tag > ... node-number-list
1852
1853 For version 4, we have:
1854 tag(int) numVert(int) ...
1855
1856 In the following, we will ignore the element number (we simply
1857 enumerate them in the order in which we read them, and we will
1858 take reg-phys (version 1) or the first tag (version 2, if any tag
1859 is given at all) as material id. For version 4, we already read
1860 the material and the cell type in above.
1861 */
1862
1863 unsigned int elm_number = 0;
1864 if (gmsh_file_format < 40)
1865 {
1866 in >> elm_number // ELM-NUMBER
1867 >> cell_type; // ELM-TYPE
1868 }
1869
1870 if (gmsh_file_format < 20)
1871 {
1872 in >> material_id // REG-PHYS
1873 >> dummy // reg_elm
1874 >> nod_num;
1875 }
1876 else if (gmsh_file_format < 40)
1877 {
1878 // read the tags; ignore all but the first one which we will
1879 // interpret as the material_id (for cells) or boundary_id
1880 // (for faces)
1881 unsigned int n_tags;
1882 in >> n_tags;
1883 if (n_tags > 0)
1884 in >> material_id;
1885 else
1886 material_id = 0;
1887
1888 for (unsigned int i = 1; i < n_tags; ++i)
1889 in >> dummy;
1890
1891 if (cell_type == 1) // line
1892 nod_num = 2;
1893 else if (cell_type == 2) // tri
1894 nod_num = 3;
1895 else if (cell_type == 3) // quad
1896 nod_num = 4;
1897 else if (cell_type == 4) // tet
1898 nod_num = 4;
1899 else if (cell_type == 5) // hex
1900 nod_num = 8;
1901 }
1902 else
1903 {
1904 // ignore tag
1905 int tag;
1906 in >> tag;
1907
1908 if (cell_type == 1) // line
1909 nod_num = 2;
1910 else if (cell_type == 2) // tri
1911 nod_num = 3;
1912 else if (cell_type == 3) // quad
1913 nod_num = 4;
1914 else if (cell_type == 4) // tet
1915 nod_num = 4;
1916 else if (cell_type == 5) // hex
1917 nod_num = 8;
1918 }
1919
1920
1921 /* `ELM-TYPE'
1922 defines the geometrical type of the N-th element:
1923 `1'
1924 Line (2 nodes, 1 edge).
1925
1926 `2'
1927 Triangle (3 nodes, 3 edges).
1928
1929 `3'
1930 Quadrangle (4 nodes, 4 edges).
1931
1932 `4'
1933 Tetrahedron (4 nodes, 6 edges, 6 faces).
1934
1935 `5'
1936 Hexahedron (8 nodes, 12 edges, 6 faces).
1937
1938 `15'
1939 Point (1 node).
1940 */
1941
1942 if (((cell_type == 1) && (dim == 1)) ||
1943 ((cell_type == 2) && (dim == 2)) ||
1944 ((cell_type == 3) && (dim == 2)) ||
1945 ((cell_type == 4) && (dim == 3)) ||
1946 ((cell_type == 5) && (dim == 3)))
1947 // found a cell
1948 {
1949 unsigned int vertices_per_cell = 0;
1950 if (cell_type == 1) // line
1951 vertices_per_cell = 2;
1952 else if (cell_type == 2) // tri
1953 {
1954 vertices_per_cell = 3;
1955 is_tria_or_tet_mesh = true;
1956 }
1957 else if (cell_type == 3) // quad
1958 {
1959 vertices_per_cell = 4;
1960 is_quad_or_hex_mesh = true;
1961 }
1962 else if (cell_type == 4) // tet
1963 {
1964 vertices_per_cell = 4;
1965 is_tria_or_tet_mesh = true;
1966 }
1967 else if (cell_type == 5) // hex
1968 {
1969 vertices_per_cell = 8;
1970 is_quad_or_hex_mesh = true;
1971 }
1972
1973 AssertThrow(nod_num == vertices_per_cell,
1974 ExcMessage(
1975 "Number of nodes does not coincide with the "
1976 "number required for this object"));
1977
1978 // allocate and read indices
1979 cells.emplace_back();
1980 cells.back().vertices.resize(vertices_per_cell);
1981 for (unsigned int i = 0; i < vertices_per_cell; ++i)
1982 in >> cells.back().vertices[i];
1983
1984 // to make sure that the cast won't fail
1985 Assert(material_id <=
1986 std::numeric_limits<types::material_id>::max(),
1987 ExcIndexRange(
1988 material_id,
1989 0,
1990 std::numeric_limits<types::material_id>::max()));
1991 // we use only material_ids in the range from 0 to
1992 // numbers::invalid_material_id-1
1993 AssertIndexRange(material_id, numbers::invalid_material_id);
1994
1995 cells.back().material_id = material_id;
1996
1997 // transform from ucd to
1998 // consecutive numbering
1999 for (unsigned int i = 0; i < vertices_per_cell; ++i)
2000 {
2001 AssertThrow(
2002 vertex_indices.find(cells.back().vertices[i]) !=
2003 vertex_indices.end(),
2004 ExcInvalidVertexIndexGmsh(cell_per_entity,
2005 elm_number,
2006 cells.back().vertices[i]));
2007
2008 // vertex with this index exists
2009 cells.back().vertices[i] =
2010 vertex_indices[cells.back().vertices[i]];
2011 }
2012 }
2013 else if ((cell_type == 1) && ((dim == 2) || (dim == 3)))
2014 // boundary info
2015 {
2016 subcelldata.boundary_lines.emplace_back();
2017 in >> subcelldata.boundary_lines.back().vertices[0] >>
2018 subcelldata.boundary_lines.back().vertices[1];
2019
2020 // to make sure that the cast won't fail
2021 Assert(material_id <=
2022 std::numeric_limits<types::boundary_id>::max(),
2023 ExcIndexRange(
2024 material_id,
2025 0,
2026 std::numeric_limits<types::boundary_id>::max()));
2027 // we use only boundary_ids in the range from 0 to
2028 // numbers::internal_face_boundary_id-1
2029 AssertIndexRange(material_id,
2030 numbers::internal_face_boundary_id);
2031
2032 subcelldata.boundary_lines.back().boundary_id =
2033 static_cast<types::boundary_id>(material_id);
2034
2035 // transform from ucd to
2036 // consecutive numbering
2037 for (unsigned int &vertex :
2038 subcelldata.boundary_lines.back().vertices)
2039 if (vertex_indices.find(vertex) != vertex_indices.end())
2040 // vertex with this index exists
2041 vertex = vertex_indices[vertex];
2042 else
2043 {
2044 // no such vertex index
2045 AssertThrow(false,
2046 ExcInvalidVertexIndex(cell_per_entity,
2047 vertex));
2048 vertex = numbers::invalid_unsigned_int;
2049 }
2050 }
2051 else if ((cell_type == 2 || cell_type == 3) && (dim == 3))
2052 // boundary info
2053 {
2054 unsigned int vertices_per_cell = 0;
2055 // check cell type
2056 if (cell_type == 2) // tri
2057 {
2058 vertices_per_cell = 3;
2059 is_tria_or_tet_mesh = true;
2060 }
2061 else if (cell_type == 3) // quad
2062 {
2063 vertices_per_cell = 4;
2064 is_quad_or_hex_mesh = true;
2065 }
2066
2067 subcelldata.boundary_quads.emplace_back();
2068
2069 // resize vertices
2070 subcelldata.boundary_quads.back().vertices.resize(
2071 vertices_per_cell);
2072 // for loop
2073 for (unsigned int i = 0; i < vertices_per_cell; ++i)
2074 in >> subcelldata.boundary_quads.back().vertices[i];
2075
2076 // to make sure that the cast won't fail
2077 Assert(material_id <=
2078 std::numeric_limits<types::boundary_id>::max(),
2079 ExcIndexRange(
2080 material_id,
2081 0,
2082 std::numeric_limits<types::boundary_id>::max()));
2083 // we use only boundary_ids in the range from 0 to
2084 // numbers::internal_face_boundary_id-1
2085 AssertIndexRange(material_id,
2086 numbers::internal_face_boundary_id);
2087
2088 subcelldata.boundary_quads.back().boundary_id =
2089 static_cast<types::boundary_id>(material_id);
2090
2091 // transform from gmsh to
2092 // consecutive numbering
2093 for (unsigned int &vertex :
2094 subcelldata.boundary_quads.back().vertices)
2095 if (vertex_indices.find(vertex) != vertex_indices.end())
2096 // vertex with this index exists
2097 vertex = vertex_indices[vertex];
2098 else
2099 {
2100 // no such vertex index
2101 Assert(false,
2102 ExcInvalidVertexIndex(cell_per_entity, vertex));
2103 vertex = numbers::invalid_unsigned_int;
2104 }
2105 }
2106 else if (cell_type == 15)
2107 {
2108 // read the indices of nodes given
2109 unsigned int node_index = 0;
2110 if (gmsh_file_format < 20)
2111 {
2112 // For points (cell_type==15), we can only ever
2113 // list one node index.
2114 AssertThrow(nod_num == 1, ExcInternalError());
2115 in >> node_index;
2116 }
2117 else
2118 {
2119 in >> node_index;
2120 }
2121
2122 // we only care about boundary indicators assigned to individual
2123 // vertices in 1d (because otherwise the vertices are not faces)
2124 if (dim == 1)
2125 boundary_ids_1d[vertex_indices[node_index]] = material_id;
2126 }
2127 else
2128 {
2129 AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type));
2130 }
2131 }
2132 }
2133 AssertDimension(global_cell, n_cells);
2134 }
2135 // Assert we reached the end of the block
2136 in >> line;
2137 static const std::string end_elements_marker[] = {"$ENDELM", "$EndElements"};
2138 AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2139 ExcInvalidGMSHInput(line));
2140
2141 // check that no forbidden arrays are used
2142 Assert(subcelldata.check_consistency(dim), ExcInternalError());
2143
2144 AssertThrow(in, ExcIO());
2145
2146 // check that we actually read some cells.
2147 AssertThrow(cells.size() > 0, ExcGmshNoCellInformation());
2148
2149 // make sure that only either simplex or hypercube cells are available
2150 //
2151 // TODO: the functions below (GridTools::delete_unused_vertices(),
2152 // GridTools::invert_all_cells_of_negative_grid(),
2153 // GridReordering::reorder_cells(),
2154 // Triangulation::create_triangulation_compatibility()) need to be revisited
2155 // for simplex meshes
2156 AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
2157 ExcNotImplemented());
2158
2159 if (dim == 1 || is_quad_or_hex_mesh)
2160 {
2161 // do some clean-up on vertices...
2162 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2163 // ... and cells
2164 if (dim == spacedim)
2165 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(
2166 vertices, cells);
2167 GridReordering<dim, spacedim>::reorder_cells(cells);
2168 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2169 }
2170 else
2171 {
2172 tria->create_triangulation(vertices, cells, subcelldata);
2173 }
2174
2175 // in 1d, we also have to attach boundary ids to vertices, which does not
2176 // currently work through the call above
2177 if (dim == 1)
2178 assign_1d_boundary_ids(boundary_ids_1d, *tria);
2179 }
2180
2181
2182
2183 template <int dim, int spacedim>
2184 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)2185 GridIn<dim, spacedim>::parse_tecplot_header(
2186 std::string & header,
2187 std::vector<unsigned int> &tecplot2deal,
2188 unsigned int & n_vars,
2189 unsigned int & n_vertices,
2190 unsigned int & n_cells,
2191 std::vector<unsigned int> &IJK,
2192 bool & structured,
2193 bool & blocked)
2194 {
2195 Assert(tecplot2deal.size() == dim, ExcInternalError());
2196 Assert(IJK.size() == dim, ExcInternalError());
2197 // initialize the output variables
2198 n_vars = 0;
2199 n_vertices = 0;
2200 n_cells = 0;
2201 switch (dim)
2202 {
2203 case 3:
2204 IJK[2] = 0;
2205 DEAL_II_FALLTHROUGH;
2206 case 2:
2207 IJK[1] = 0;
2208 DEAL_II_FALLTHROUGH;
2209 case 1:
2210 IJK[0] = 0;
2211 }
2212 structured = true;
2213 blocked = false;
2214
2215 // convert the string to upper case
2216 std::transform(header.begin(), header.end(), header.begin(), ::toupper);
2217
2218 // replace all tabs, commas, newlines by
2219 // whitespaces
2220 std::replace(header.begin(), header.end(), '\t', ' ');
2221 std::replace(header.begin(), header.end(), ',', ' ');
2222 std::replace(header.begin(), header.end(), '\n', ' ');
2223
2224 // now remove whitespace in front of and
2225 // after '='
2226 std::string::size_type pos = header.find('=');
2227
2228 while (pos != static_cast<std::string::size_type>(std::string::npos))
2229 if (header[pos + 1] == ' ')
2230 header.erase(pos + 1, 1);
2231 else if (header[pos - 1] == ' ')
2232 {
2233 header.erase(pos - 1, 1);
2234 --pos;
2235 }
2236 else
2237 pos = header.find('=', ++pos);
2238
2239 // split the string into individual entries
2240 std::vector<std::string> entries =
2241 Utilities::break_text_into_lines(header, 1, ' ');
2242
2243 // now go through the list and try to extract
2244 for (unsigned int i = 0; i < entries.size(); ++i)
2245 {
2246 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\""))
2247 {
2248 ++n_vars;
2249 // we assume, that the first variable
2250 // is x or no coordinate at all (not y or z)
2251 if (Utilities::match_at_string_start(entries[i], "VARIABLES=\"X\""))
2252 {
2253 tecplot2deal[0] = 0;
2254 }
2255 ++i;
2256 while (entries[i][0] == '"')
2257 {
2258 if (entries[i] == "\"X\"")
2259 tecplot2deal[0] = n_vars;
2260 else if (entries[i] == "\"Y\"")
2261 {
2262 // we assume, that y contains
2263 // zero data in 1d, so do
2264 // nothing
2265 if (dim > 1)
2266 tecplot2deal[1] = n_vars;
2267 }
2268 else if (entries[i] == "\"Z\"")
2269 {
2270 // we assume, that z contains
2271 // zero data in 1d and 2d, so
2272 // do nothing
2273 if (dim > 2)
2274 tecplot2deal[2] = n_vars;
2275 }
2276 ++n_vars;
2277 ++i;
2278 }
2279 // set i back, so that the next
2280 // string is treated correctly
2281 --i;
2282
2283 AssertThrow(
2284 n_vars >= dim,
2285 ExcMessage(
2286 "Tecplot file must contain at least one variable for each dimension"));
2287 for (unsigned int d = 1; d < dim; ++d)
2288 AssertThrow(
2289 tecplot2deal[d] > 0,
2290 ExcMessage(
2291 "Tecplot file must contain at least one variable for each dimension."));
2292 }
2293 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE=ORDERED"))
2294 structured = true;
2295 else if (Utilities::match_at_string_start(entries[i],
2296 "ZONETYPE=FELINESEG") &&
2297 dim == 1)
2298 structured = false;
2299 else if (Utilities::match_at_string_start(entries[i],
2300 "ZONETYPE=FEQUADRILATERAL") &&
2301 dim == 2)
2302 structured = false;
2303 else if (Utilities::match_at_string_start(entries[i],
2304 "ZONETYPE=FEBRICK") &&
2305 dim == 3)
2306 structured = false;
2307 else if (Utilities::match_at_string_start(entries[i], "ZONETYPE="))
2308 // unsupported ZONETYPE
2309 {
2310 AssertThrow(false,
2311 ExcMessage(
2312 "The tecplot file contains an unsupported ZONETYPE."));
2313 }
2314 else if (Utilities::match_at_string_start(entries[i],
2315 "DATAPACKING=POINT"))
2316 blocked = false;
2317 else if (Utilities::match_at_string_start(entries[i],
2318 "DATAPACKING=BLOCK"))
2319 blocked = true;
2320 else if (Utilities::match_at_string_start(entries[i], "F=POINT"))
2321 {
2322 structured = true;
2323 blocked = false;
2324 }
2325 else if (Utilities::match_at_string_start(entries[i], "F=BLOCK"))
2326 {
2327 structured = true;
2328 blocked = true;
2329 }
2330 else if (Utilities::match_at_string_start(entries[i], "F=FEPOINT"))
2331 {
2332 structured = false;
2333 blocked = false;
2334 }
2335 else if (Utilities::match_at_string_start(entries[i], "F=FEBLOCK"))
2336 {
2337 structured = false;
2338 blocked = true;
2339 }
2340 else if (Utilities::match_at_string_start(entries[i],
2341 "ET=QUADRILATERAL") &&
2342 dim == 2)
2343 structured = false;
2344 else if (Utilities::match_at_string_start(entries[i], "ET=BRICK") &&
2345 dim == 3)
2346 structured = false;
2347 else if (Utilities::match_at_string_start(entries[i], "ET="))
2348 // unsupported ElementType
2349 {
2350 AssertThrow(
2351 false,
2352 ExcMessage(
2353 "The tecplot file contains an unsupported ElementType."));
2354 }
2355 else if (Utilities::match_at_string_start(entries[i], "I="))
2356 IJK[0] = Utilities::get_integer_at_position(entries[i], 2).first;
2357 else if (Utilities::match_at_string_start(entries[i], "J="))
2358 {
2359 IJK[1] = Utilities::get_integer_at_position(entries[i], 2).first;
2360 AssertThrow(
2361 dim > 1 || IJK[1] == 1,
2362 ExcMessage(
2363 "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
2364 }
2365 else if (Utilities::match_at_string_start(entries[i], "K="))
2366 {
2367 IJK[2] = Utilities::get_integer_at_position(entries[i], 2).first;
2368 AssertThrow(
2369 dim > 2 || IJK[2] == 1,
2370 ExcMessage(
2371 "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
2372 }
2373 else if (Utilities::match_at_string_start(entries[i], "N="))
2374 n_vertices = Utilities::get_integer_at_position(entries[i], 2).first;
2375 else if (Utilities::match_at_string_start(entries[i], "E="))
2376 n_cells = Utilities::get_integer_at_position(entries[i], 2).first;
2377 }
2378
2379 // now we have read all the fields we are
2380 // interested in. do some checks and
2381 // calculate the variables
2382 if (structured)
2383 {
2384 n_vertices = 1;
2385 n_cells = 1;
2386 for (unsigned int d = 0; d < dim; ++d)
2387 {
2388 AssertThrow(
2389 IJK[d] > 0,
2390 ExcMessage(
2391 "Tecplot file does not contain a complete and consistent set of parameters"));
2392 n_vertices *= IJK[d];
2393 n_cells *= (IJK[d] - 1);
2394 }
2395 }
2396 else
2397 {
2398 AssertThrow(
2399 n_vertices > 0,
2400 ExcMessage(
2401 "Tecplot file does not contain a complete and consistent set of parameters"));
2402 if (n_cells == 0)
2403 // this means an error, although
2404 // tecplot itself accepts entries like
2405 // 'J=20' instead of 'E=20'. therefore,
2406 // take the max of IJK
2407 n_cells = *std::max_element(IJK.begin(), IJK.end());
2408 AssertThrow(
2409 n_cells > 0,
2410 ExcMessage(
2411 "Tecplot file does not contain a complete and consistent set of parameters"));
2412 }
2413 }
2414
2415
2416
2417 template <>
2418 void
read_tecplot(std::istream & in)2419 GridIn<2>::read_tecplot(std::istream &in)
2420 {
2421 const unsigned int dim = 2;
2422 const unsigned int spacedim = 2;
2423 Assert(tria != nullptr, ExcNoTriangulationSelected());
2424 AssertThrow(in, ExcIO());
2425
2426 // skip comments at start of file
2427 skip_comment_lines(in, '#');
2428
2429 // some strings for parsing the header
2430 std::string line, header;
2431
2432 // first, concatenate all header lines
2433 // create a searchstring with almost all
2434 // letters. exclude e and E from the letters
2435 // to search, as they might appear in
2436 // exponential notation
2437 std::string letters = "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
2438
2439 getline(in, line);
2440 while (line.find_first_of(letters) != std::string::npos)
2441 {
2442 header += " " + line;
2443 getline(in, line);
2444 }
2445
2446 // now create some variables holding
2447 // important information on the mesh, get
2448 // this information from the header string
2449 std::vector<unsigned int> tecplot2deal(dim);
2450 std::vector<unsigned int> IJK(dim);
2451 unsigned int n_vars, n_vertices, n_cells;
2452 bool structured, blocked;
2453
2454 parse_tecplot_header(header,
2455 tecplot2deal,
2456 n_vars,
2457 n_vertices,
2458 n_cells,
2459 IJK,
2460 structured,
2461 blocked);
2462
2463 // reserve space for vertices. note, that in
2464 // tecplot vertices are ordered beginning
2465 // with 1, whereas in deal all indices start
2466 // with 0. in order not to use -1 for all the
2467 // connectivity information, a 0th vertex
2468 // (unused) is inserted at the origin.
2469 std::vector<Point<spacedim>> vertices(n_vertices + 1);
2470 vertices[0] = Point<spacedim>();
2471 // reserve space for cells
2472 std::vector<CellData<dim>> cells(n_cells);
2473 SubCellData subcelldata;
2474
2475 if (blocked)
2476 {
2477 // blocked data format. first we get all
2478 // the values of the first variable for
2479 // all points, after that we get all
2480 // values for the second variable and so
2481 // on.
2482
2483 // dummy variable to read in all the info
2484 // we do not want to use
2485 double dummy;
2486 // which is the first index to read in
2487 // the loop (see below)
2488 unsigned int next_index = 0;
2489
2490 // note, that we have already read the
2491 // first line containing the first variable
2492 if (tecplot2deal[0] == 0)
2493 {
2494 // we need the information in this
2495 // line, so extract it
2496 std::vector<std::string> first_var =
2497 Utilities::break_text_into_lines(line, 1);
2498 char *endptr;
2499 for (unsigned int i = 1; i < first_var.size() + 1; ++i)
2500 vertices[i](0) = std::strtod(first_var[i - 1].c_str(), &endptr);
2501
2502 // if there are many points, the data
2503 // for this var might continue in the
2504 // next line(s)
2505 for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
2506 in >> vertices[j](next_index);
2507 // now we got all values of the first
2508 // variable, so increase the counter
2509 next_index = 1;
2510 }
2511
2512 // main loop over all variables
2513 for (unsigned int i = 1; i < n_vars; ++i)
2514 {
2515 // if we read all the important
2516 // variables and do not want to
2517 // read further, because we are
2518 // using a structured grid, we can
2519 // stop here (and skip, for
2520 // example, a whole lot of solution
2521 // variables)
2522 if (next_index == dim && structured)
2523 break;
2524
2525 if ((next_index < dim) && (i == tecplot2deal[next_index]))
2526 {
2527 // we need this line, read it in
2528 for (unsigned int j = 1; j < n_vertices + 1; ++j)
2529 in >> vertices[j](next_index);
2530 ++next_index;
2531 }
2532 else
2533 {
2534 // we do not need this line, read
2535 // it in and discard it
2536 for (unsigned int j = 1; j < n_vertices + 1; ++j)
2537 in >> dummy;
2538 }
2539 }
2540 Assert(next_index == dim, ExcInternalError());
2541 }
2542 else
2543 {
2544 // the data is not blocked, so we get all
2545 // the variables for one point, then the
2546 // next and so on. create a vector to
2547 // hold these components
2548 std::vector<double> vars(n_vars);
2549
2550 // now fill the first vertex. note, that we
2551 // have already read the first line
2552 // containing the first vertex
2553 std::vector<std::string> first_vertex =
2554 Utilities::break_text_into_lines(line, 1);
2555 char *endptr;
2556 for (unsigned int d = 0; d < dim; ++d)
2557 vertices[1](d) =
2558 std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
2559
2560 // read the remaining vertices from the
2561 // list
2562 for (unsigned int v = 2; v < n_vertices + 1; ++v)
2563 {
2564 for (unsigned int i = 0; i < n_vars; ++i)
2565 in >> vars[i];
2566 // fill the vertex
2567 // coordinates. respect the position
2568 // of coordinates in the list of
2569 // variables
2570 for (unsigned int i = 0; i < dim; ++i)
2571 vertices[v](i) = vars[tecplot2deal[i]];
2572 }
2573 }
2574
2575 if (structured)
2576 {
2577 // this is the part of the code that only
2578 // works in 2d
2579 unsigned int I = IJK[0], J = IJK[1];
2580
2581 unsigned int cell = 0;
2582 // set up array of cells
2583 for (unsigned int j = 0; j < J - 1; ++j)
2584 for (unsigned int i = 1; i < I; ++i)
2585 {
2586 cells[cell].vertices[0] = i + j * I;
2587 cells[cell].vertices[1] = i + 1 + j * I;
2588 cells[cell].vertices[2] = i + 1 + (j + 1) * I;
2589 cells[cell].vertices[3] = i + (j + 1) * I;
2590 ++cell;
2591 }
2592 Assert(cell == n_cells, ExcInternalError());
2593 std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
2594 unsigned int k = 0;
2595 for (unsigned int i = 1; i < I + 1; ++i)
2596 {
2597 boundary_vertices[k] = i;
2598 ++k;
2599 boundary_vertices[k] = i + (J - 1) * I;
2600 ++k;
2601 }
2602 for (unsigned int j = 1; j < J - 1; ++j)
2603 {
2604 boundary_vertices[k] = 1 + j * I;
2605 ++k;
2606 boundary_vertices[k] = I + j * I;
2607 ++k;
2608 }
2609 Assert(k == boundary_vertices.size(), ExcInternalError());
2610 // delete the duplicated vertices at the
2611 // boundary, which occur, e.g. in c-type
2612 // or o-type grids around a body
2613 // (airfoil). this automatically deletes
2614 // unused vertices as well.
2615 GridTools::delete_duplicated_vertices(vertices,
2616 cells,
2617 subcelldata,
2618 boundary_vertices);
2619 }
2620 else
2621 {
2622 // set up array of cells, unstructured
2623 // mode, so the connectivity is
2624 // explicitly given
2625 for (unsigned int i = 0; i < n_cells; ++i)
2626 {
2627 // note that since in the input file
2628 // we found the number of cells at
2629 // the top, there should still be
2630 // input here, so check this:
2631 AssertThrow(in, ExcIO());
2632
2633 // get the connectivity from the
2634 // input file. the vertices are
2635 // ordered like in the ucd format
2636 for (unsigned int &vertex : cells[i].vertices)
2637 in >> vertex;
2638 }
2639 // do some clean-up on vertices
2640 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2641 }
2642
2643 // check that no forbidden arrays are
2644 // used. as we do not read in any
2645 // subcelldata, nothing should happen here.
2646 Assert(subcelldata.check_consistency(dim), ExcInternalError());
2647 AssertThrow(in, ExcIO());
2648
2649 // do some cleanup on cells
2650 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
2651 cells);
2652 GridReordering<dim, spacedim>::reorder_cells(cells);
2653 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2654 }
2655
2656
2657
2658 template <int dim, int spacedim>
2659 void
read_tecplot(std::istream &)2660 GridIn<dim, spacedim>::read_tecplot(std::istream &)
2661 {
2662 Assert(false, ExcNotImplemented());
2663 }
2664
2665
2666
2667 template <int dim, int spacedim>
2668 void
read_assimp(const std::string & filename,const unsigned int mesh_index,const bool remove_duplicates,const double tol,const bool ignore_unsupported_types)2669 GridIn<dim, spacedim>::read_assimp(const std::string &filename,
2670 const unsigned int mesh_index,
2671 const bool remove_duplicates,
2672 const double tol,
2673 const bool ignore_unsupported_types)
2674 {
2675 #ifdef DEAL_II_WITH_ASSIMP
2676 // Only good for surface grids.
2677 AssertThrow(dim < 3, ExcImpossibleInDim(dim));
2678
2679 // Create an instance of the Importer class
2680 Assimp::Importer importer;
2681
2682 // And have it read the given file with some postprocessing
2683 const aiScene *scene =
2684 importer.ReadFile(filename.c_str(),
2685 aiProcess_RemoveComponent |
2686 aiProcess_JoinIdenticalVertices |
2687 aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
2688 aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
2689
2690 // If the import failed, report it
2691 AssertThrow(scene != nullptr, ExcMessage(importer.GetErrorString()));
2692
2693 AssertThrow(scene->mNumMeshes != 0,
2694 ExcMessage("Input file contains no meshes."));
2695
2696 AssertThrow((mesh_index == numbers::invalid_unsigned_int) ||
2697 (mesh_index < scene->mNumMeshes),
2698 ExcMessage("Too few meshes in the file."));
2699
2700 unsigned int start_mesh =
2701 (mesh_index == numbers::invalid_unsigned_int ? 0 : mesh_index);
2702 unsigned int end_mesh =
2703 (mesh_index == numbers::invalid_unsigned_int ? scene->mNumMeshes :
2704 mesh_index + 1);
2705
2706 // Deal.II objects are created empty, and then filled with imported file.
2707 std::vector<Point<spacedim>> vertices;
2708 std::vector<CellData<dim>> cells;
2709 SubCellData subcelldata;
2710
2711 // A series of counters to merge cells.
2712 unsigned int v_offset = 0;
2713 unsigned int c_offset = 0;
2714
2715 // The index of the mesh will be used as a material index.
2716 for (unsigned int m = start_mesh; m < end_mesh; ++m)
2717 {
2718 const aiMesh *mesh = scene->mMeshes[m];
2719
2720 // Check that we know what to do with this mesh, otherwise just
2721 // ignore it
2722 if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
2723 {
2724 AssertThrow(ignore_unsupported_types,
2725 ExcMessage("Incompatible mesh " + std::to_string(m) +
2726 "/" + std::to_string(scene->mNumMeshes)));
2727 continue;
2728 }
2729 else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
2730 {
2731 AssertThrow(ignore_unsupported_types,
2732 ExcMessage("Incompatible mesh " + std::to_string(m) +
2733 "/" + std::to_string(scene->mNumMeshes)));
2734 continue;
2735 }
2736 // Vertices
2737 const unsigned int n_vertices = mesh->mNumVertices;
2738 const aiVector3D * mVertices = mesh->mVertices;
2739
2740 // Faces
2741 const unsigned int n_faces = mesh->mNumFaces;
2742 const aiFace * mFaces = mesh->mFaces;
2743
2744 vertices.resize(v_offset + n_vertices);
2745 cells.resize(c_offset + n_faces);
2746
2747 for (unsigned int i = 0; i < n_vertices; ++i)
2748 for (unsigned int d = 0; d < spacedim; ++d)
2749 vertices[i + v_offset][d] = mVertices[i][d];
2750
2751 unsigned int valid_cell = c_offset;
2752 for (unsigned int i = 0; i < n_faces; ++i)
2753 {
2754 if (mFaces[i].mNumIndices == GeometryInfo<dim>::vertices_per_cell)
2755 {
2756 for (const unsigned int f : GeometryInfo<dim>::vertex_indices())
2757 {
2758 cells[valid_cell].vertices[f] =
2759 mFaces[i].mIndices[f] + v_offset;
2760 }
2761 cells[valid_cell].material_id = m;
2762 ++valid_cell;
2763 }
2764 else
2765 {
2766 AssertThrow(ignore_unsupported_types,
2767 ExcMessage("Face " + std::to_string(i) + " of mesh " +
2768 std::to_string(m) + " has " +
2769 std::to_string(mFaces[i].mNumIndices) +
2770 " vertices. We expected only " +
2771 std::to_string(
2772 GeometryInfo<dim>::vertices_per_cell)));
2773 }
2774 }
2775 cells.resize(valid_cell);
2776
2777 // The vertices are added all at once. Cells are checked for
2778 // validity, so only valid_cells are now present in the deal.II
2779 // list of cells.
2780 v_offset += n_vertices;
2781 c_offset = valid_cell;
2782 }
2783
2784 // No cells were read
2785 if (cells.size() == 0)
2786 return;
2787
2788 if (remove_duplicates)
2789 {
2790 // The function delete_duplicated_vertices() needs to be called more
2791 // than once if a vertex is duplicated more than once. So we keep
2792 // calling it until the number of vertices does not change any more.
2793 unsigned int n_verts = 0;
2794 while (n_verts != vertices.size())
2795 {
2796 n_verts = vertices.size();
2797 std::vector<unsigned int> considered_vertices;
2798 GridTools::delete_duplicated_vertices(
2799 vertices, cells, subcelldata, considered_vertices, tol);
2800 }
2801 }
2802
2803 GridTools::delete_unused_vertices(vertices, cells, subcelldata);
2804 if (dim == spacedim)
2805 GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
2806 cells);
2807
2808 GridReordering<dim, spacedim>::reorder_cells(cells);
2809 if (dim == 2)
2810 tria->create_triangulation_compatibility(vertices, cells, subcelldata);
2811 else
2812 tria->create_triangulation(vertices, cells, subcelldata);
2813 #else
2814 (void)filename;
2815 (void)mesh_index;
2816 (void)remove_duplicates;
2817 (void)tol;
2818 (void)ignore_unsupported_types;
2819 AssertThrow(false, ExcNeedsAssimp());
2820 #endif
2821 }
2822
2823
2824 template <int dim, int spacedim>
2825 void
skip_empty_lines(std::istream & in)2826 GridIn<dim, spacedim>::skip_empty_lines(std::istream &in)
2827 {
2828 std::string line;
2829 while (in)
2830 {
2831 // get line
2832 getline(in, line);
2833
2834 // check if this is a line that
2835 // consists only of spaces, and
2836 // if not put the whole thing
2837 // back and return
2838 if (std::find_if(line.begin(), line.end(), [](const char c) {
2839 return c != ' ';
2840 }) != line.end())
2841 {
2842 in.putback('\n');
2843 for (int i = line.length() - 1; i >= 0; --i)
2844 in.putback(line[i]);
2845 return;
2846 }
2847
2848 // else: go on with next line
2849 }
2850 }
2851
2852
2853
2854 template <int dim, int spacedim>
2855 void
skip_comment_lines(std::istream & in,const char comment_start)2856 GridIn<dim, spacedim>::skip_comment_lines(std::istream &in,
2857 const char comment_start)
2858 {
2859 char c;
2860 // loop over the following comment
2861 // lines
2862 while (in.get(c) && c == comment_start)
2863 // loop over the characters after
2864 // the comment starter
2865 while (in.get() != '\n')
2866 ;
2867
2868
2869 // put back first character of
2870 // first non-comment line
2871 if (in)
2872 in.putback(c);
2873
2874 // at last: skip additional empty lines, if present
2875 skip_empty_lines(in);
2876 }
2877
2878
2879
2880 template <int dim, int spacedim>
2881 void
debug_output_grid(const std::vector<CellData<dim>> &,const std::vector<Point<spacedim>> &,std::ostream &)2882 GridIn<dim, spacedim>::debug_output_grid(
2883 const std::vector<CellData<dim>> & /*cells*/,
2884 const std::vector<Point<spacedim>> & /*vertices*/,
2885 std::ostream & /*out*/)
2886 {
2887 Assert(false, ExcNotImplemented());
2888 }
2889
2890
2891
2892 template <>
2893 void
debug_output_grid(const std::vector<CellData<2>> & cells,const std::vector<Point<2>> & vertices,std::ostream & out)2894 GridIn<2>::debug_output_grid(const std::vector<CellData<2>> &cells,
2895 const std::vector<Point<2>> & vertices,
2896 std::ostream & out)
2897 {
2898 double min_x = vertices[cells[0].vertices[0]](0),
2899 max_x = vertices[cells[0].vertices[0]](0),
2900 min_y = vertices[cells[0].vertices[0]](1),
2901 max_y = vertices[cells[0].vertices[0]](1);
2902
2903 for (unsigned int i = 0; i < cells.size(); ++i)
2904 {
2905 for (const auto vertex : cells[i].vertices)
2906 {
2907 const Point<2> &p = vertices[vertex];
2908
2909 if (p(0) < min_x)
2910 min_x = p(0);
2911 if (p(0) > max_x)
2912 max_x = p(0);
2913 if (p(1) < min_y)
2914 min_y = p(1);
2915 if (p(1) > max_y)
2916 max_y = p(1);
2917 }
2918
2919 out << "# cell " << i << std::endl;
2920 Point<2> center;
2921 for (const auto vertex : cells[i].vertices)
2922 center += vertices[vertex];
2923 center /= 4;
2924
2925 out << "set label \"" << i << "\" at " << center(0) << ',' << center(1)
2926 << " center" << std::endl;
2927
2928 // first two line right direction
2929 for (unsigned int f = 0; f < 2; ++f)
2930 out << "set arrow from " << vertices[cells[i].vertices[f]](0) << ','
2931 << vertices[cells[i].vertices[f]](1) << " to "
2932 << vertices[cells[i].vertices[(f + 1) % 4]](0) << ','
2933 << vertices[cells[i].vertices[(f + 1) % 4]](1) << std::endl;
2934 // other two lines reverse direction
2935 for (unsigned int f = 2; f < 4; ++f)
2936 out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]](0)
2937 << ',' << vertices[cells[i].vertices[(f + 1) % 4]](1) << " to "
2938 << vertices[cells[i].vertices[f]](0) << ','
2939 << vertices[cells[i].vertices[f]](1) << std::endl;
2940 out << std::endl;
2941 }
2942
2943
2944 out << std::endl
2945 << "set nokey" << std::endl
2946 << "pl [" << min_x << ':' << max_x << "][" << min_y << ':' << max_y
2947 << "] " << min_y << std::endl
2948 << "pause -1" << std::endl;
2949 }
2950
2951
2952
2953 template <>
2954 void
debug_output_grid(const std::vector<CellData<3>> & cells,const std::vector<Point<3>> & vertices,std::ostream & out)2955 GridIn<3>::debug_output_grid(const std::vector<CellData<3>> &cells,
2956 const std::vector<Point<3>> & vertices,
2957 std::ostream & out)
2958 {
2959 for (const auto &cell : cells)
2960 {
2961 // line 0
2962 out << vertices[cell.vertices[0]] << std::endl
2963 << vertices[cell.vertices[1]] << std::endl
2964 << std::endl
2965 << std::endl;
2966 // line 1
2967 out << vertices[cell.vertices[1]] << std::endl
2968 << vertices[cell.vertices[2]] << std::endl
2969 << std::endl
2970 << std::endl;
2971 // line 2
2972 out << vertices[cell.vertices[3]] << std::endl
2973 << vertices[cell.vertices[2]] << std::endl
2974 << std::endl
2975 << std::endl;
2976 // line 3
2977 out << vertices[cell.vertices[0]] << std::endl
2978 << vertices[cell.vertices[3]] << std::endl
2979 << std::endl
2980 << std::endl;
2981 // line 4
2982 out << vertices[cell.vertices[4]] << std::endl
2983 << vertices[cell.vertices[5]] << std::endl
2984 << std::endl
2985 << std::endl;
2986 // line 5
2987 out << vertices[cell.vertices[5]] << std::endl
2988 << vertices[cell.vertices[6]] << std::endl
2989 << std::endl
2990 << std::endl;
2991 // line 6
2992 out << vertices[cell.vertices[7]] << std::endl
2993 << vertices[cell.vertices[6]] << std::endl
2994 << std::endl
2995 << std::endl;
2996 // line 7
2997 out << vertices[cell.vertices[4]] << std::endl
2998 << vertices[cell.vertices[7]] << std::endl
2999 << std::endl
3000 << std::endl;
3001 // line 8
3002 out << vertices[cell.vertices[0]] << std::endl
3003 << vertices[cell.vertices[4]] << std::endl
3004 << std::endl
3005 << std::endl;
3006 // line 9
3007 out << vertices[cell.vertices[1]] << std::endl
3008 << vertices[cell.vertices[5]] << std::endl
3009 << std::endl
3010 << std::endl;
3011 // line 10
3012 out << vertices[cell.vertices[2]] << std::endl
3013 << vertices[cell.vertices[6]] << std::endl
3014 << std::endl
3015 << std::endl;
3016 // line 11
3017 out << vertices[cell.vertices[3]] << std::endl
3018 << vertices[cell.vertices[7]] << std::endl
3019 << std::endl
3020 << std::endl;
3021 }
3022 }
3023
3024
3025
3026 template <int dim, int spacedim>
3027 void
read(const std::string & filename,Format format)3028 GridIn<dim, spacedim>::read(const std::string &filename, Format format)
3029 {
3030 // Search file class for meshes
3031 PathSearch search("MESH");
3032 std::string name;
3033 // Open the file and remember its name
3034 if (format == Default)
3035 name = search.find(filename);
3036 else
3037 name = search.find(filename, default_suffix(format));
3038
3039
3040 if (format == Default)
3041 {
3042 const std::string::size_type slashpos = name.find_last_of('/');
3043 const std::string::size_type dotpos = name.find_last_of('.');
3044 if (dotpos < name.length() &&
3045 (dotpos > slashpos || slashpos == std::string::npos))
3046 {
3047 std::string ext = name.substr(dotpos + 1);
3048 format = parse_format(ext);
3049 }
3050 }
3051
3052 if (format == assimp)
3053 {
3054 read_assimp(name);
3055 }
3056 else
3057 {
3058 std::ifstream in(name.c_str());
3059 read(in, format);
3060 }
3061 }
3062
3063
3064 template <int dim, int spacedim>
3065 void
read(std::istream & in,Format format)3066 GridIn<dim, spacedim>::read(std::istream &in, Format format)
3067 {
3068 if (format == Default)
3069 format = default_format;
3070
3071 switch (format)
3072 {
3073 case dbmesh:
3074 read_dbmesh(in);
3075 return;
3076
3077 case msh:
3078 read_msh(in);
3079 return;
3080
3081 case vtk:
3082 read_vtk(in);
3083 return;
3084
3085 case vtu:
3086 read_vtu(in);
3087 return;
3088
3089 case unv:
3090 read_unv(in);
3091 return;
3092
3093 case ucd:
3094 read_ucd(in);
3095 return;
3096
3097 case abaqus:
3098 read_abaqus(in);
3099 return;
3100
3101 case xda:
3102 read_xda(in);
3103 return;
3104
3105 case tecplot:
3106 read_tecplot(in);
3107 return;
3108
3109 case assimp:
3110 Assert(false,
3111 ExcMessage("There is no read_assimp(istream &) function. "
3112 "Use the read_assimp(string &filename, ...) "
3113 "functions, instead."));
3114 return;
3115
3116 case Default:
3117 break;
3118 }
3119 Assert(false, ExcInternalError());
3120 }
3121
3122
3123
3124 template <int dim, int spacedim>
3125 std::string
default_suffix(const Format format)3126 GridIn<dim, spacedim>::default_suffix(const Format format)
3127 {
3128 switch (format)
3129 {
3130 case dbmesh:
3131 return ".dbmesh";
3132 case msh:
3133 return ".msh";
3134 case vtk:
3135 return ".vtk";
3136 case vtu:
3137 return ".vtu";
3138 case unv:
3139 return ".unv";
3140 case ucd:
3141 return ".inp";
3142 case abaqus:
3143 return ".inp"; // Typical suffix for Abaqus mesh files conflicts with
3144 // UCD.
3145 case xda:
3146 return ".xda";
3147 case tecplot:
3148 return ".dat";
3149 default:
3150 Assert(false, ExcNotImplemented());
3151 return ".unknown_format";
3152 }
3153 }
3154
3155
3156
3157 template <int dim, int spacedim>
3158 typename GridIn<dim, spacedim>::Format
parse_format(const std::string & format_name)3159 GridIn<dim, spacedim>::parse_format(const std::string &format_name)
3160 {
3161 if (format_name == "dbmesh")
3162 return dbmesh;
3163
3164 if (format_name == "msh")
3165 return msh;
3166
3167 if (format_name == "unv")
3168 return unv;
3169
3170 if (format_name == "vtk")
3171 return vtk;
3172
3173 if (format_name == "vtu")
3174 return vtu;
3175
3176 // This is also the typical extension of Abaqus input files.
3177 if (format_name == "inp")
3178 return ucd;
3179
3180 if (format_name == "ucd")
3181 return ucd;
3182
3183 if (format_name == "xda")
3184 return xda;
3185
3186 if (format_name == "tecplot")
3187 return tecplot;
3188
3189 if (format_name == "dat")
3190 return tecplot;
3191
3192 if (format_name == "plt")
3193 // Actually, this is the extension for the
3194 // tecplot binary format, which we do not
3195 // support right now. However, some people
3196 // tend to create tecplot ascii files with
3197 // the extension 'plt' instead of
3198 // 'dat'. Thus, include this extension
3199 // here. If it actually is a binary file,
3200 // the read_tecplot() function will fail
3201 // and throw an exception, anyway.
3202 return tecplot;
3203
3204 AssertThrow(false, ExcInvalidState());
3205 // return something weird
3206 return Format(Default);
3207 }
3208
3209
3210
3211 template <int dim, int spacedim>
3212 std::string
get_format_names()3213 GridIn<dim, spacedim>::get_format_names()
3214 {
3215 return "dbmesh|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
3216 }
3217
3218
3219
3220 namespace
3221 {
3222 template <int dim, int spacedim>
Abaqus_to_UCD()3223 Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
3224 : tolerance(5e-16) // Used to offset Cubit tolerance error when outputting
3225 // value close to zero
3226 {
3227 AssertThrow(spacedim == 2 || spacedim == 3, ExcNotImplemented());
3228 }
3229
3230
3231
3232 // Convert from a string to some other data type
3233 // Reference: http://www.codeguru.com/forum/showthread.php?t=231054
3234 template <class T>
3235 bool
from_string(T & t,const std::string & s,std::ios_base & (* f)(std::ios_base &))3236 from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &))
3237 {
3238 std::istringstream iss(s);
3239 return !(iss >> f >> t).fail();
3240 }
3241
3242
3243
3244 // Extract an integer from a string
3245 int
extract_int(const std::string & s)3246 extract_int(const std::string &s)
3247 {
3248 std::string tmp;
3249 for (const char c : s)
3250 {
3251 if (isdigit(c))
3252 {
3253 tmp += c;
3254 }
3255 }
3256
3257 int number = 0;
3258 from_string(number, tmp, std::dec);
3259 return number;
3260 }
3261
3262
3263
3264 template <int dim, int spacedim>
3265 void
read_in_abaqus(std::istream & input_stream)3266 Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
3267 {
3268 // References:
3269 // http://www.egr.msu.edu/software/abaqus/Documentation/docs/v6.7/books/usb/default.htm?startat=pt01ch02.html
3270 // http://www.cprogramming.com/tutorial/string.html
3271
3272 AssertThrow(input_stream, ExcIO());
3273 std::string line;
3274
3275 while (std::getline(input_stream, line))
3276 {
3277 cont:
3278 std::transform(line.begin(), line.end(), line.begin(), ::toupper);
3279
3280 if (line.compare("*HEADING") == 0 || line.compare(0, 2, "**") == 0 ||
3281 line.compare(0, 5, "*PART") == 0)
3282 {
3283 // Skip header and comments
3284 while (std::getline(input_stream, line))
3285 {
3286 if (line[0] == '*')
3287 goto cont; // My eyes, they burn!
3288 }
3289 }
3290 else if (line.compare(0, 5, "*NODE") == 0)
3291 {
3292 // Extract list of vertices
3293 // Header line might be:
3294 // *NODE, NSET=ALLNODES
3295 // *NODE
3296
3297 // Contains lines in the form:
3298 // Index, x, y, z
3299 while (std::getline(input_stream, line))
3300 {
3301 if (line[0] == '*')
3302 goto cont;
3303
3304 std::vector<double> node(spacedim + 1);
3305
3306 std::istringstream iss(line);
3307 char comma;
3308 for (unsigned int i = 0; i < spacedim + 1; ++i)
3309 iss >> node[i] >> comma;
3310
3311 node_list.push_back(node);
3312 }
3313 }
3314 else if (line.compare(0, 8, "*ELEMENT") == 0)
3315 {
3316 // Element construction.
3317 // There are different header formats, the details
3318 // of which we're not particularly interested in except
3319 // whether they represent quads or hexahedrals.
3320 // *ELEMENT, TYPE=S4R, ELSET=EB<material id>
3321 // *ELEMENT, TYPE=C3D8R, ELSET=EB<material id>
3322 // *ELEMENT, TYPE=C3D8
3323 // Elements itself (n=4 or n=8):
3324 // Index, i[0], ..., i[n]
3325
3326 int material = 0;
3327 // Scan for material id
3328 {
3329 const std::string before_material = "ELSET=EB";
3330 const std::size_t idx = line.find(before_material);
3331 if (idx != std::string::npos)
3332 {
3333 from_string(material,
3334 line.substr(idx + before_material.size()),
3335 std::dec);
3336 }
3337 }
3338
3339 // Read ELEMENT definition
3340 while (std::getline(input_stream, line))
3341 {
3342 if (line[0] == '*')
3343 goto cont;
3344
3345 std::istringstream iss(line);
3346 char comma;
3347
3348 // We will store the material id in the zeroth entry of the
3349 // vector and the rest of the elements represent the global
3350 // node numbers
3351 const unsigned int n_data_per_cell =
3352 1 + GeometryInfo<dim>::vertices_per_cell;
3353 std::vector<double> cell(n_data_per_cell);
3354 for (unsigned int i = 0; i < n_data_per_cell; ++i)
3355 iss >> cell[i] >> comma;
3356
3357 // Overwrite cell index from file by material
3358 cell[0] = static_cast<double>(material);
3359 cell_list.push_back(cell);
3360 }
3361 }
3362 else if (line.compare(0, 8, "*SURFACE") == 0)
3363 {
3364 // Extract the definitions of boundary surfaces
3365 // Old format from Cubit:
3366 // *SURFACE, NAME=SS<boundary indicator>
3367 // <element index>, S<face number>
3368 // Abaqus default format:
3369 // *SURFACE, TYPE=ELEMENT, NAME=SURF-<indicator>
3370
3371 // Get name of the surface and extract id from it;
3372 // this will be the boundary indicator
3373 const std::string name_key = "NAME=";
3374 const std::size_t name_idx_start =
3375 line.find(name_key) + name_key.size();
3376 std::size_t name_idx_end = line.find(',', name_idx_start);
3377 if (name_idx_end == std::string::npos)
3378 {
3379 name_idx_end = line.size();
3380 }
3381 const int b_indicator = extract_int(
3382 line.substr(name_idx_start, name_idx_end - name_idx_start));
3383
3384 // Read SURFACE definition
3385 // Note that the orientation of the faces is embedded within the
3386 // definition of each "set" of faces that comprise the surface
3387 // These are either marked by an "S" or "E" in 3d or 2d
3388 // respectively.
3389 while (std::getline(input_stream, line))
3390 {
3391 if (line[0] == '*')
3392 goto cont;
3393
3394 // Change all characters to upper case
3395 std::transform(line.begin(),
3396 line.end(),
3397 line.begin(),
3398 ::toupper);
3399
3400 // Surface can be created from ELSET, or directly from cells
3401 // If elsets_list contains a key with specific name - refers to
3402 // that ELSET, otherwise refers to cell
3403 std::istringstream iss(line);
3404 int el_idx;
3405 int face_number;
3406 char temp;
3407
3408 // Get relevant faces, taking into account the element
3409 // orientation
3410 std::vector<double> quad_node_list;
3411 const std::string elset_name = line.substr(0, line.find(','));
3412 if (elsets_list.count(elset_name) != 0)
3413 {
3414 // Surface refers to ELSET
3415 std::string stmp;
3416 iss >> stmp >> temp >> face_number;
3417
3418 const std::vector<int> cells = elsets_list[elset_name];
3419 for (const int cell : cells)
3420 {
3421 el_idx = cell;
3422 quad_node_list =
3423 get_global_node_numbers(el_idx, face_number);
3424 quad_node_list.insert(quad_node_list.begin(),
3425 b_indicator);
3426
3427 face_list.push_back(quad_node_list);
3428 }
3429 }
3430 else
3431 {
3432 // Surface refers directly to elements
3433 char comma;
3434 iss >> el_idx >> comma >> temp >> face_number;
3435 quad_node_list =
3436 get_global_node_numbers(el_idx, face_number);
3437 quad_node_list.insert(quad_node_list.begin(), b_indicator);
3438
3439 face_list.push_back(quad_node_list);
3440 }
3441 }
3442 }
3443 else if (line.compare(0, 6, "*ELSET") == 0)
3444 {
3445 // Get ELSET name.
3446 // Materials are attached to elsets with specific name
3447 std::string elset_name;
3448 {
3449 const std::string elset_key = "*ELSET, ELSET=";
3450 const std::size_t idx = line.find(elset_key);
3451 if (idx != std::string::npos)
3452 {
3453 const std::string comma = ",";
3454 const std::size_t first_comma = line.find(comma);
3455 const std::size_t second_comma =
3456 line.find(comma, first_comma + 1);
3457 const std::size_t elset_name_start =
3458 line.find(elset_key) + elset_key.size();
3459 elset_name = line.substr(elset_name_start,
3460 second_comma - elset_name_start);
3461 }
3462 }
3463
3464 // There are two possibilities of storing cells numbers in ELSET:
3465 // 1. If the header contains the 'GENERATE' keyword, then the next
3466 // line describes range of cells as:
3467 // cell_id_start, cell_id_end, cell_step
3468 // 2. If the header does not contain the 'GENERATE' keyword, then
3469 // the next lines contain cells numbers
3470 std::vector<int> elements;
3471 const std::size_t generate_idx = line.find("GENERATE");
3472 if (generate_idx != std::string::npos)
3473 {
3474 // Option (1)
3475 std::getline(input_stream, line);
3476 std::istringstream iss(line);
3477 char comma;
3478 int elid_start;
3479 int elid_end;
3480 int elis_step = 1; // Default if case stride not provided
3481
3482 // Some files don't have the stride size
3483 // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp to
3484 // ./grids/abaqus/2d/2d_test_abaqus.inp
3485 iss >> elid_start >> comma >> elid_end;
3486 AssertThrow(comma == ',',
3487 ExcMessage(
3488 std::string(
3489 "While reading an ABAQUS file, the reader "
3490 "expected a comma but found a <") +
3491 comma + "> in the line <" + line + ">."));
3492 AssertThrow(
3493 elid_start <= elid_end,
3494 ExcMessage(
3495 std::string(
3496 "While reading an ABAQUS file, the reader encountered "
3497 "a GENERATE statement in which the upper bound <") +
3498 Utilities::int_to_string(elid_end) +
3499 "> for the element numbers is not larger or equal "
3500 "than the lower bound <" +
3501 Utilities::int_to_string(elid_start) + ">."));
3502
3503 // https://stackoverflow.com/questions/8046357/how-do-i-check-if-a-stringstream-variable-is-empty-null
3504 if (iss.rdbuf()->in_avail() != 0)
3505 iss >> comma >> elis_step;
3506 AssertThrow(comma == ',',
3507 ExcMessage(
3508 std::string(
3509 "While reading an ABAQUS file, the reader "
3510 "expected a comma but found a <") +
3511 comma + "> in the line <" + line + ">."));
3512
3513 for (int i = elid_start; i <= elid_end; i += elis_step)
3514 elements.push_back(i);
3515 elsets_list[elset_name] = elements;
3516
3517 std::getline(input_stream, line);
3518 }
3519 else
3520 {
3521 // Option (2)
3522 while (std::getline(input_stream, line))
3523 {
3524 if (line[0] == '*')
3525 break;
3526
3527 std::istringstream iss(line);
3528 char comma;
3529 int elid;
3530 while (!iss.eof())
3531 {
3532 iss >> elid >> comma;
3533 AssertThrow(
3534 comma == ',',
3535 ExcMessage(
3536 std::string(
3537 "While reading an ABAQUS file, the reader "
3538 "expected a comma but found a <") +
3539 comma + "> in the line <" + line + ">."));
3540
3541 elements.push_back(elid);
3542 }
3543 }
3544
3545 elsets_list[elset_name] = elements;
3546 }
3547
3548 goto cont;
3549 }
3550 else if (line.compare(0, 5, "*NSET") == 0)
3551 {
3552 // Skip nodesets; we have no use for them
3553 while (std::getline(input_stream, line))
3554 {
3555 if (line[0] == '*')
3556 goto cont;
3557 }
3558 }
3559 else if (line.compare(0, 14, "*SOLID SECTION") == 0)
3560 {
3561 // The ELSET name, which describes a section for particular material
3562 const std::string elset_key = "ELSET=";
3563 const std::size_t elset_start =
3564 line.find("ELSET=") + elset_key.size();
3565 const std::size_t elset_end = line.find(',', elset_start + 1);
3566 const std::string elset_name =
3567 line.substr(elset_start, elset_end - elset_start);
3568
3569 // Solid material definition.
3570 // We assume that material id is taken from material name,
3571 // eg. "Material-1" -> ID=1
3572 const std::string material_key = "MATERIAL=";
3573 const std::size_t last_equal =
3574 line.find("MATERIAL=") + material_key.size();
3575 const std::size_t material_id_start = line.find('-', last_equal);
3576 int material_id = 0;
3577 from_string(material_id,
3578 line.substr(material_id_start + 1),
3579 std::dec);
3580
3581 // Assign material id to cells
3582 const std::vector<int> &elset_cells = elsets_list[elset_name];
3583 for (const int elset_cell : elset_cells)
3584 {
3585 const int cell_id = elset_cell - 1;
3586 cell_list[cell_id][0] = material_id;
3587 }
3588 }
3589 // Note: All other lines / entries are ignored
3590 }
3591 }
3592
3593 template <int dim, int spacedim>
3594 std::vector<double>
get_global_node_numbers(const int face_cell_no,const int face_cell_face_no) const3595 Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
3596 const int face_cell_no,
3597 const int face_cell_face_no) const
3598 {
3599 std::vector<double> quad_node_list(GeometryInfo<dim>::vertices_per_face);
3600
3601 // These orderings were reverse engineered by hand and may
3602 // conceivably be erroneous.
3603 // TODO: Currently one test (2d unstructured mesh) in the test
3604 // suite fails, presumably because of an ordering issue.
3605 if (dim == 2)
3606 {
3607 if (face_cell_face_no == 1)
3608 {
3609 quad_node_list[0] = cell_list[face_cell_no - 1][1];
3610 quad_node_list[1] = cell_list[face_cell_no - 1][2];
3611 }
3612 else if (face_cell_face_no == 2)
3613 {
3614 quad_node_list[0] = cell_list[face_cell_no - 1][2];
3615 quad_node_list[1] = cell_list[face_cell_no - 1][3];
3616 }
3617 else if (face_cell_face_no == 3)
3618 {
3619 quad_node_list[0] = cell_list[face_cell_no - 1][3];
3620 quad_node_list[1] = cell_list[face_cell_no - 1][4];
3621 }
3622 else if (face_cell_face_no == 4)
3623 {
3624 quad_node_list[0] = cell_list[face_cell_no - 1][4];
3625 quad_node_list[1] = cell_list[face_cell_no - 1][1];
3626 }
3627 else
3628 {
3629 AssertThrow(face_cell_face_no <= 4,
3630 ExcMessage("Invalid face number in 2d"));
3631 }
3632 }
3633 else if (dim == 3)
3634 {
3635 if (face_cell_face_no == 1)
3636 {
3637 quad_node_list[0] = cell_list[face_cell_no - 1][1];
3638 quad_node_list[1] = cell_list[face_cell_no - 1][4];
3639 quad_node_list[2] = cell_list[face_cell_no - 1][3];
3640 quad_node_list[3] = cell_list[face_cell_no - 1][2];
3641 }
3642 else if (face_cell_face_no == 2)
3643 {
3644 quad_node_list[0] = cell_list[face_cell_no - 1][5];
3645 quad_node_list[1] = cell_list[face_cell_no - 1][8];
3646 quad_node_list[2] = cell_list[face_cell_no - 1][7];
3647 quad_node_list[3] = cell_list[face_cell_no - 1][6];
3648 }
3649 else if (face_cell_face_no == 3)
3650 {
3651 quad_node_list[0] = cell_list[face_cell_no - 1][1];
3652 quad_node_list[1] = cell_list[face_cell_no - 1][2];
3653 quad_node_list[2] = cell_list[face_cell_no - 1][6];
3654 quad_node_list[3] = cell_list[face_cell_no - 1][5];
3655 }
3656 else if (face_cell_face_no == 4)
3657 {
3658 quad_node_list[0] = cell_list[face_cell_no - 1][2];
3659 quad_node_list[1] = cell_list[face_cell_no - 1][3];
3660 quad_node_list[2] = cell_list[face_cell_no - 1][7];
3661 quad_node_list[3] = cell_list[face_cell_no - 1][6];
3662 }
3663 else if (face_cell_face_no == 5)
3664 {
3665 quad_node_list[0] = cell_list[face_cell_no - 1][3];
3666 quad_node_list[1] = cell_list[face_cell_no - 1][4];
3667 quad_node_list[2] = cell_list[face_cell_no - 1][8];
3668 quad_node_list[3] = cell_list[face_cell_no - 1][7];
3669 }
3670 else if (face_cell_face_no == 6)
3671 {
3672 quad_node_list[0] = cell_list[face_cell_no - 1][1];
3673 quad_node_list[1] = cell_list[face_cell_no - 1][5];
3674 quad_node_list[2] = cell_list[face_cell_no - 1][8];
3675 quad_node_list[3] = cell_list[face_cell_no - 1][4];
3676 }
3677 else
3678 {
3679 AssertThrow(face_cell_no <= 6,
3680 ExcMessage("Invalid face number in 3d"));
3681 }
3682 }
3683 else
3684 {
3685 AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
3686 }
3687
3688 return quad_node_list;
3689 }
3690
3691 template <int dim, int spacedim>
3692 void
write_out_avs_ucd(std::ostream & output) const3693 Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output) const
3694 {
3695 // References:
3696 // http://www.dealii.org/developer/doxygen/deal.II/structGeometryInfo.html
3697 // http://people.scs.fsu.edu/~burkardt/data/ucd/ucd.html
3698
3699 AssertThrow(output, ExcIO());
3700
3701 // save old formatting options
3702 const boost::io::ios_base_all_saver formatting_saver(output);
3703
3704 // Write out title - Note: No other commented text can be inserted below the
3705 // title in a UCD file
3706 output << "# Abaqus to UCD mesh conversion" << std::endl;
3707 output << "# Mesh type: AVS UCD" << std::endl;
3708
3709 // ========================================================
3710 // ASCII UCD File Format
3711 // The input file cannot contain blank lines or lines with leading blanks.
3712 // Comments, if present, must precede all data in the file.
3713 // Comments within the data will cause read errors.
3714 // The general order of the data is as follows:
3715 // 1. Numbers defining the overall structure, including the number of nodes,
3716 // the number of cells, and the length of the vector of data associated
3717 // with the nodes, cells, and the model.
3718 // e.g. 1:
3719 // <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
3720 // e.g. 2:
3721 // n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
3722 // outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
3723 // 2. For each node, its node id and the coordinates of that node in space.
3724 // Node-ids must be integers, but any number including non sequential
3725 // numbers can be used. Mid-edge nodes are treated like any other node.
3726 // 3. For each cell: its cell-id, material, cell type (hexahedral, pyramid,
3727 // etc.), and the list of node-ids that correspond to each of the cell's
3728 // vertices. The below table specifies the different cell types and the
3729 // keyword used to represent them in the file.
3730
3731 // Write out header
3732 output << node_list.size() << "\t" << (cell_list.size() + face_list.size())
3733 << "\t0\t0\t0" << std::endl;
3734
3735 output.width(16);
3736 output.precision(8);
3737
3738 // Write out node numbers
3739 // Loop over all nodes
3740 for (const auto &node : node_list)
3741 {
3742 // Node number
3743 output << node[0] << "\t";
3744
3745 // Node coordinates
3746 output.setf(std::ios::scientific, std::ios::floatfield);
3747 for (unsigned int jj = 1; jj < spacedim + 1; ++jj)
3748 {
3749 // invoke tolerance -> set points close to zero equal to zero
3750 if (std::abs(node[jj]) > tolerance)
3751 output << static_cast<double>(node[jj]) << "\t";
3752 else
3753 output << 0.0 << "\t";
3754 }
3755 if (spacedim == 2)
3756 output << 0.0 << "\t";
3757
3758 output << std::endl;
3759 output.unsetf(std::ios::floatfield);
3760 }
3761
3762 // Write out cell node numbers
3763 for (unsigned int ii = 0; ii < cell_list.size(); ++ii)
3764 {
3765 output << ii + 1 << "\t" << cell_list[ii][0] << "\t"
3766 << (dim == 2 ? "quad" : "hex") << "\t";
3767 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
3768 ++jj)
3769 output << cell_list[ii][jj] << "\t";
3770
3771 output << std::endl;
3772 }
3773
3774 // Write out quad node numbers
3775 for (unsigned int ii = 0; ii < face_list.size(); ++ii)
3776 {
3777 output << ii + 1 << "\t" << face_list[ii][0] << "\t"
3778 << (dim == 2 ? "line" : "quad") << "\t";
3779 for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
3780 ++jj)
3781 output << face_list[ii][jj] << "\t";
3782
3783 output << std::endl;
3784 }
3785 }
3786 } // namespace
3787
3788
3789 // explicit instantiations
3790 #include "grid_in.inst"
3791
3792 DEAL_II_NAMESPACE_CLOSE
3793