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