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