MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh_readers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "mesh_headers.hpp"
13 #include "../fem/fem.hpp"
14 #include "../general/text.hpp"
15 
16 #include <iostream>
17 #include <cstdio>
18 
19 #ifdef MFEM_USE_NETCDF
20 #include "netcdf.h"
21 #endif
22 
23 using namespace std;
24 
25 namespace mfem
26 {
27 
28 bool Mesh::remove_unused_vertices = true;
29 
30 void Mesh::ReadMFEMMesh(std::istream &input, bool mfem_v11, int &curved)
31 {
32  // Read MFEM mesh v1.0 format
33  string ident;
34 
35  // read lines beginning with '#' (comments)
36  skip_comment_lines(input, '#');
37  input >> ident; // 'dimension'
38 
39  MFEM_VERIFY(ident == "dimension", "invalid mesh file");
40  input >> Dim;
41 
42  skip_comment_lines(input, '#');
43  input >> ident; // 'elements'
44 
45  MFEM_VERIFY(ident == "elements", "invalid mesh file");
46  input >> NumOfElements;
47  elements.SetSize(NumOfElements);
48  for (int j = 0; j < NumOfElements; j++)
49  {
50  elements[j] = ReadElement(input);
51  }
52 
53  skip_comment_lines(input, '#');
54  input >> ident; // 'boundary'
55 
56  MFEM_VERIFY(ident == "boundary", "invalid mesh file");
57  input >> NumOfBdrElements;
58  boundary.SetSize(NumOfBdrElements);
59  for (int j = 0; j < NumOfBdrElements; j++)
60  {
61  boundary[j] = ReadElement(input);
62  }
63 
64  skip_comment_lines(input, '#');
65  input >> ident;
66 
67  if (mfem_v11 && ident == "vertex_parents")
68  {
69  ncmesh = new NCMesh(this, &input);
70  // NOTE: the constructor above will call LoadVertexParents
71 
72  skip_comment_lines(input, '#');
73  input >> ident;
74 
75  if (ident == "coarse_elements")
76  {
77  ncmesh->LoadCoarseElements(input);
78 
79  skip_comment_lines(input, '#');
80  input >> ident;
81  }
82  }
83 
84  MFEM_VERIFY(ident == "vertices", "invalid mesh file");
85  input >> NumOfVertices;
86  vertices.SetSize(NumOfVertices);
87 
88  input >> ws >> ident;
89  if (ident != "nodes")
90  {
91  // read the vertices
92  spaceDim = atoi(ident.c_str());
93  for (int j = 0; j < NumOfVertices; j++)
94  {
95  for (int i = 0; i < spaceDim; i++)
96  {
97  input >> vertices[j](i);
98  }
99  }
100 
101  // initialize vertex positions in NCMesh
102  if (ncmesh) { ncmesh->SetVertexPositions(vertices); }
103  }
104  else
105  {
106  // prepare to read the nodes
107  input >> ws;
108  curved = 1;
109  }
110 
111  // When visualizing solutions on non-conforming grids, PETSc
112  // may dump additional vertices
113  if (remove_unused_vertices) { RemoveUnusedVertices(); }
114 }
115 
116 void Mesh::ReadLineMesh(std::istream &input)
117 {
118  int j,p1,p2,a;
119 
120  Dim = 1;
121 
122  input >> NumOfVertices;
123  vertices.SetSize(NumOfVertices);
124  // Sets vertices and the corresponding coordinates
125  for (j = 0; j < NumOfVertices; j++)
126  {
127  input >> vertices[j](0);
128  }
129 
130  input >> NumOfElements;
131  elements.SetSize(NumOfElements);
132  // Sets elements and the corresponding indices of vertices
133  for (j = 0; j < NumOfElements; j++)
134  {
135  input >> a >> p1 >> p2;
136  elements[j] = new Segment(p1-1, p2-1, a);
137  }
138 
139  int ind[1];
140  input >> NumOfBdrElements;
141  boundary.SetSize(NumOfBdrElements);
142  for (j = 0; j < NumOfBdrElements; j++)
143  {
144  input >> a >> ind[0];
145  ind[0]--;
146  boundary[j] = new Point(ind,a);
147  }
148 }
149 
150 void Mesh::ReadNetgen2DMesh(std::istream &input, int &curved)
151 {
152  int ints[32], attr, n;
153 
154  // Read planar mesh in Netgen format.
155  Dim = 2;
156 
157  // Read the boundary elements.
158  input >> NumOfBdrElements;
159  boundary.SetSize(NumOfBdrElements);
160  for (int i = 0; i < NumOfBdrElements; i++)
161  {
162  input >> attr
163  >> ints[0] >> ints[1];
164  ints[0]--; ints[1]--;
165  boundary[i] = new Segment(ints, attr);
166  }
167 
168  // Read the elements.
169  input >> NumOfElements;
170  elements.SetSize(NumOfElements);
171  for (int i = 0; i < NumOfElements; i++)
172  {
173  input >> attr >> n;
174  for (int j = 0; j < n; j++)
175  {
176  input >> ints[j];
177  ints[j]--;
178  }
179  switch (n)
180  {
181  case 2:
182  elements[i] = new Segment(ints, attr);
183  break;
184  case 3:
185  elements[i] = new Triangle(ints, attr);
186  break;
187  case 4:
188  elements[i] = new Quadrilateral(ints, attr);
189  break;
190  }
191  }
192 
193  if (!curved)
194  {
195  // Read the vertices.
196  input >> NumOfVertices;
197  vertices.SetSize(NumOfVertices);
198  for (int i = 0; i < NumOfVertices; i++)
199  for (int j = 0; j < Dim; j++)
200  {
201  input >> vertices[i](j);
202  }
203  }
204  else
205  {
206  input >> NumOfVertices;
207  vertices.SetSize(NumOfVertices);
208  input >> ws;
209  }
210 }
211 
212 void Mesh::ReadNetgen3DMesh(std::istream &input)
213 {
214  int ints[32], attr;
215 
216  // Read a Netgen format mesh of tetrahedra.
217  Dim = 3;
218 
219  // Read the vertices
220  input >> NumOfVertices;
221 
222  vertices.SetSize(NumOfVertices);
223  for (int i = 0; i < NumOfVertices; i++)
224  for (int j = 0; j < Dim; j++)
225  {
226  input >> vertices[i](j);
227  }
228 
229  // Read the elements
230  input >> NumOfElements;
231  elements.SetSize(NumOfElements);
232  for (int i = 0; i < NumOfElements; i++)
233  {
234  input >> attr;
235  for (int j = 0; j < 4; j++)
236  {
237  input >> ints[j];
238  ints[j]--;
239  }
240 #ifdef MFEM_USE_MEMALLOC
241  Tetrahedron *tet;
242  tet = TetMemory.Alloc();
243  tet->SetVertices(ints);
244  tet->SetAttribute(attr);
245  elements[i] = tet;
246 #else
247  elements[i] = new Tetrahedron(ints, attr);
248 #endif
249  }
250 
251  // Read the boundary information.
252  input >> NumOfBdrElements;
253  boundary.SetSize(NumOfBdrElements);
254  for (int i = 0; i < NumOfBdrElements; i++)
255  {
256  input >> attr;
257  for (int j = 0; j < 3; j++)
258  {
259  input >> ints[j];
260  ints[j]--;
261  }
262  boundary[i] = new Triangle(ints, attr);
263  }
264 }
265 
266 void Mesh::ReadTrueGridMesh(std::istream &input)
267 {
268  int i, j, ints[32], attr;
269  const int buflen = 1024;
270  char buf[buflen];
271 
272  // TODO: find the actual dimension
273  Dim = 3;
274 
275  if (Dim == 2)
276  {
277  int vari;
278  double varf;
279 
280  input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
281  input.getline(buf, buflen);
282  input.getline(buf, buflen);
283  input >> vari;
284  input.getline(buf, buflen);
285  input.getline(buf, buflen);
286  input.getline(buf, buflen);
287 
288  // Read the vertices.
289  vertices.SetSize(NumOfVertices);
290  for (i = 0; i < NumOfVertices; i++)
291  {
292  input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
293  input.getline(buf, buflen);
294  }
295 
296  // Read the elements.
297  elements.SetSize(NumOfElements);
298  for (i = 0; i < NumOfElements; i++)
299  {
300  input >> vari >> attr;
301  for (j = 0; j < 4; j++)
302  {
303  input >> ints[j];
304  ints[j]--;
305  }
306  input.getline(buf, buflen);
307  input.getline(buf, buflen);
308  elements[i] = new Quadrilateral(ints, attr);
309  }
310  }
311  else if (Dim == 3)
312  {
313  int vari;
314  double varf;
315  input >> vari >> NumOfVertices >> NumOfElements;
316  input.getline(buf, buflen);
317  input.getline(buf, buflen);
318  input >> vari >> vari >> NumOfBdrElements;
319  input.getline(buf, buflen);
320  input.getline(buf, buflen);
321  input.getline(buf, buflen);
322  // Read the vertices.
323  vertices.SetSize(NumOfVertices);
324  for (i = 0; i < NumOfVertices; i++)
325  {
326  input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
327  >> vertices[i](2);
328  input.getline(buf, buflen);
329  }
330  // Read the elements.
331  elements.SetSize(NumOfElements);
332  for (i = 0; i < NumOfElements; i++)
333  {
334  input >> vari >> attr;
335  for (j = 0; j < 8; j++)
336  {
337  input >> ints[j];
338  ints[j]--;
339  }
340  input.getline(buf, buflen);
341  elements[i] = new Hexahedron(ints, attr);
342  }
343  // Read the boundary elements.
344  boundary.SetSize(NumOfBdrElements);
345  for (i = 0; i < NumOfBdrElements; i++)
346  {
347  input >> attr;
348  for (j = 0; j < 4; j++)
349  {
350  input >> ints[j];
351  ints[j]--;
352  }
353  input.getline(buf, buflen);
354  boundary[i] = new Quadrilateral(ints, attr);
355  }
356  }
357 }
358 
359 // see Tetrahedron::edges
360 const int Mesh::vtk_quadratic_tet[10] =
361 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
362 
363 // see Wedge::edges & Mesh::GenerateFaces
364 // https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
365 const int Mesh::vtk_quadratic_wedge[18] =
366 { 0, 2, 1, 3, 5, 4, 8, 7, 6, 11, 10, 9, 12, 14, 13, 17, 16, 15};
367 
368 // see Hexahedron::edges & Mesh::GenerateFaces
369 const int Mesh::vtk_quadratic_hex[27] =
370 {
371  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
372  24, 22, 21, 23, 20, 25, 26
373 };
374 
375 void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
376  bool &finalize_topo)
377 {
378  // VTK resources:
379  // * https://www.vtk.org/doc/nightly/html/vtkCellType_8h_source.html
380  // * https://www.vtk.org/doc/nightly/html/classvtkCell.html
381  // * https://lorensen.github.io/VTKExamples/site/VTKFileFormats
382  // * https://www.kitware.com/products/books/VTKUsersGuide.pdf
383 
384  int i, j, n, attr;
385 
386  string buff;
387  getline(input, buff); // comment line
388  getline(input, buff);
389  filter_dos(buff);
390  if (buff != "ASCII")
391  {
392  MFEM_ABORT("VTK mesh is not in ASCII format!");
393  return;
394  }
395  getline(input, buff);
396  filter_dos(buff);
397  if (buff != "DATASET UNSTRUCTURED_GRID")
398  {
399  MFEM_ABORT("VTK mesh is not UNSTRUCTURED_GRID!");
400  return;
401  }
402 
403  // Read the points, skipping optional sections such as the FIELD data from
404  // VisIt's VTK export (or from Mesh::PrintVTK with field_data==1).
405  do
406  {
407  input >> buff;
408  if (!input.good())
409  {
410  MFEM_ABORT("VTK mesh does not have POINTS data!");
411  }
412  }
413  while (buff != "POINTS");
414  int np = 0;
415  Vector points;
416  {
417  input >> np >> ws;
418  points.SetSize(3*np);
419  getline(input, buff); // "double"
420  for (i = 0; i < points.Size(); i++)
421  {
422  input >> points(i);
423  }
424  }
425 
426  // Read the cells
427  NumOfElements = n = 0;
428  Array<int> cells_data;
429  input >> ws >> buff;
430  if (buff == "CELLS")
431  {
432  input >> NumOfElements >> n >> ws;
433  cells_data.SetSize(n);
434  for (i = 0; i < n; i++)
435  {
436  input >> cells_data[i];
437  }
438  }
439 
440  // Read the cell types
441  Dim = -1;
442  int order = -1;
443  input >> ws >> buff;
444  if (buff == "CELL_TYPES")
445  {
446  input >> NumOfElements;
447  elements.SetSize(NumOfElements);
448  for (j = i = 0; i < NumOfElements; i++)
449  {
450  int ct, elem_dim, elem_order = 1;
451  input >> ct;
452  switch (ct)
453  {
454  case 5: // triangle
455  elem_dim = 2;
456  elements[i] = new Triangle(&cells_data[j+1]);
457  break;
458  case 9: // quadrilateral
459  elem_dim = 2;
460  elements[i] = new Quadrilateral(&cells_data[j+1]);
461  break;
462  case 10: // tetrahedron
463  elem_dim = 3;
464 #ifdef MFEM_USE_MEMALLOC
465  elements[i] = TetMemory.Alloc();
466  elements[i]->SetVertices(&cells_data[j+1]);
467 #else
468  elements[i] = new Tetrahedron(&cells_data[j+1]);
469 #endif
470  break;
471  case 12: // hexahedron
472  elem_dim = 3;
473  elements[i] = new Hexahedron(&cells_data[j+1]);
474  break;
475  case 13: // wedge
476  elem_dim = 3;
477  // switch between vtk vertex ordering and mfem vertex ordering:
478  // swap vertices (1,2) and (4,5)
479  elements[i] =
480  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
481  cells_data[j+4], cells_data[j+6], cells_data[j+5]);
482  break;
483 
484  case 22: // quadratic triangle
485  elem_dim = 2;
486  elem_order = 2;
487  elements[i] = new Triangle(&cells_data[j+1]);
488  break;
489  case 28: // biquadratic quadrilateral
490  elem_dim = 2;
491  elem_order = 2;
492  elements[i] = new Quadrilateral(&cells_data[j+1]);
493  break;
494  case 24: // quadratic tetrahedron
495  elem_dim = 3;
496  elem_order = 2;
497 #ifdef MFEM_USE_MEMALLOC
498  elements[i] = TetMemory.Alloc();
499  elements[i]->SetVertices(&cells_data[j+1]);
500 #else
501  elements[i] = new Tetrahedron(&cells_data[j+1]);
502 #endif
503  break;
504  case 32: // biquadratic-quadratic wedge
505  elem_dim = 3;
506  elem_order = 2;
507  // switch between vtk vertex ordering and mfem vertex ordering:
508  // swap vertices (1,2) and (4,5)
509  elements[i] =
510  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
511  cells_data[j+4], cells_data[j+6], cells_data[j+5]);
512  break;
513  case 29: // triquadratic hexahedron
514  elem_dim = 3;
515  elem_order = 2;
516  elements[i] = new Hexahedron(&cells_data[j+1]);
517  break;
518  default:
519  MFEM_ABORT("VTK mesh : cell type " << ct << " is not supported!");
520  return;
521  }
522  MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
523  "elements with different dimensions are not supported");
524  MFEM_VERIFY(order == -1 || order == elem_order,
525  "elements with different orders are not supported");
526  Dim = elem_dim;
527  order = elem_order;
528  j += cells_data[j] + 1;
529  }
530  }
531 
532  // Read attributes
533  streampos sp = input.tellg();
534  input >> ws >> buff;
535  if (buff == "CELL_DATA")
536  {
537  input >> n >> ws;
538  getline(input, buff);
539  filter_dos(buff);
540  // "SCALARS material dataType numComp"
541  if (!strncmp(buff.c_str(), "SCALARS material", 16))
542  {
543  getline(input, buff); // "LOOKUP_TABLE default"
544  for (i = 0; i < NumOfElements; i++)
545  {
546  input >> attr;
547  elements[i]->SetAttribute(attr);
548  }
549  }
550  else
551  {
552  input.seekg(sp);
553  }
554  }
555  else
556  {
557  input.seekg(sp);
558  }
559 
560  if (order == 1)
561  {
562  cells_data.DeleteAll();
563  NumOfVertices = np;
564  vertices.SetSize(np);
565  for (i = 0; i < np; i++)
566  {
567  vertices[i](0) = points(3*i+0);
568  vertices[i](1) = points(3*i+1);
569  vertices[i](2) = points(3*i+2);
570  }
571  points.Destroy();
572 
573  // No boundary is defined in a VTK mesh
574  NumOfBdrElements = 0;
575  }
576  else if (order == 2)
577  {
578  curved = 1;
579 
580  // generate new enumeration for the vertices
581  Array<int> pts_dof(np);
582  pts_dof = -1;
583  for (n = i = 0; i < NumOfElements; i++)
584  {
585  int *v = elements[i]->GetVertices();
586  int nv = elements[i]->GetNVertices();
587  for (j = 0; j < nv; j++)
588  if (pts_dof[v[j]] == -1)
589  {
590  pts_dof[v[j]] = n++;
591  }
592  }
593  // keep the original ordering of the vertices
594  for (n = i = 0; i < np; i++)
595  if (pts_dof[i] != -1)
596  {
597  pts_dof[i] = n++;
598  }
599  // update the element vertices
600  for (i = 0; i < NumOfElements; i++)
601  {
602  int *v = elements[i]->GetVertices();
603  int nv = elements[i]->GetNVertices();
604  for (j = 0; j < nv; j++)
605  {
606  v[j] = pts_dof[v[j]];
607  }
608  }
609  // Define the 'vertices' from the 'points' through the 'pts_dof' map
610  NumOfVertices = n;
611  vertices.SetSize(n);
612  for (i = 0; i < np; i++)
613  {
614  if ((j = pts_dof[i]) != -1)
615  {
616  vertices[j](0) = points(3*i+0);
617  vertices[j](1) = points(3*i+1);
618  vertices[j](2) = points(3*i+2);
619  }
620  }
621 
622  // No boundary is defined in a VTK mesh
623  NumOfBdrElements = 0;
624 
625  // Generate faces and edges so that we can define quadratic
626  // FE space on the mesh
627  FinalizeTopology();
628  finalize_topo = false;
629 
630  // Define quadratic FE space
632  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim);
633  Nodes = new GridFunction(fes);
634  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
635  own_nodes = 1;
636 
637  // Map vtk points to edge/face/element dofs
638  Array<int> dofs;
639  for (n = i = 0; i < NumOfElements; i++)
640  {
641  fes->GetElementDofs(i, dofs);
642  const int *vtk_mfem;
643  switch (elements[i]->GetGeometryType())
644  {
645  case Geometry::TRIANGLE:
646  case Geometry::SQUARE:
647  vtk_mfem = vtk_quadratic_hex; break; // identity map
648  case Geometry::TETRAHEDRON:
649  vtk_mfem = vtk_quadratic_tet; break;
650  case Geometry::CUBE:
651  vtk_mfem = vtk_quadratic_hex; break;
652  case Geometry::PRISM:
653  vtk_mfem = vtk_quadratic_wedge; break;
654  default:
655  vtk_mfem = NULL; // suppress a warning
656  break;
657  }
658 
659  for (n++, j = 0; j < dofs.Size(); j++, n++)
660  {
661  if (pts_dof[cells_data[n]] == -1)
662  {
663  pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
664  }
665  else
666  {
667  if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
668  {
669  MFEM_ABORT("VTK mesh : inconsistent quadratic mesh!");
670  }
671  }
672  }
673  }
674 
675  // Define the 'Nodes' from the 'points' through the 'pts_dof' map
676  for (i = 0; i < np; i++)
677  {
678  dofs.SetSize(1);
679  if ((dofs[0] = pts_dof[i]) != -1)
680  {
681  fes->DofsToVDofs(dofs);
682  for (j = 0; j < dofs.Size(); j++)
683  {
684  (*Nodes)(dofs[j]) = points(3*i+j);
685  }
686  }
687  }
688 
689  read_gf = 0;
690  }
691 }
692 
693 void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
694 {
695  NURBSext = new NURBSExtension(input);
696 
697  Dim = NURBSext->Dimension();
698  NumOfVertices = NURBSext->GetNV();
699  NumOfElements = NURBSext->GetNE();
700  NumOfBdrElements = NURBSext->GetNBE();
701 
702  NURBSext->GetElementTopo(elements);
703  NURBSext->GetBdrElementTopo(boundary);
704 
705  vertices.SetSize(NumOfVertices);
706  curved = 1;
707  if (NURBSext->HavePatches())
708  {
709  NURBSFECollection *fec = new NURBSFECollection(NURBSext->GetOrder());
710  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
711  Ordering::byVDIM);
712  Nodes = new GridFunction(fes);
713  Nodes->MakeOwner(fec);
714  NURBSext->SetCoordsFromPatches(*Nodes);
715  own_nodes = 1;
716  read_gf = 0;
717  int vd = Nodes->VectorDim();
718  for (int i = 0; i < vd; i++)
719  {
720  Vector vert_val;
721  Nodes->GetNodalValues(vert_val, i+1);
722  for (int j = 0; j < NumOfVertices; j++)
723  {
724  vertices[j](i) = vert_val(j);
725  }
726  }
727  }
728  else
729  {
730  read_gf = 1;
731  }
732 }
733 
734 void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
735 {
736  // Initialize to negative numbers so that we know if they've been set. We're
737  // using Element::POINT as our flag, since we're not going to make a 0D mesh,
738  // ever.
739  int nx = -1;
740  int ny = -1;
741  int nz = -1;
742  double sx = -1.0;
743  double sy = -1.0;
744  double sz = -1.0;
745  Element::Type type = Element::POINT;
746 
747  while (true)
748  {
749  skip_comment_lines(input, '#');
750  // Break out if we reached the end of the file after gobbling up the
751  // whitespace and comments after the last keyword.
752  if (!input.good())
753  {
754  break;
755  }
756 
757  // Read the next keyword
758  std::string name;
759  input >> name;
760  input >> std::ws;
761  // Make sure there's an equal sign
762  MFEM_VERIFY(input.get() == '=',
763  "Inline mesh expected '=' after keyword " << name);
764  input >> std::ws;
765 
766  if (name == "nx")
767  {
768  input >> nx;
769  }
770  else if (name == "ny")
771  {
772  input >> ny;
773  }
774  else if (name == "nz")
775  {
776  input >> nz;
777  }
778  else if (name == "sx")
779  {
780  input >> sx;
781  }
782  else if (name == "sy")
783  {
784  input >> sy;
785  }
786  else if (name == "sz")
787  {
788  input >> sz;
789  }
790  else if (name == "type")
791  {
792  std::string eltype;
793  input >> eltype;
794  if (eltype == "segment")
795  {
796  type = Element::SEGMENT;
797  }
798  else if (eltype == "quad")
799  {
800  type = Element::QUADRILATERAL;
801  }
802  else if (eltype == "tri")
803  {
804  type = Element::TRIANGLE;
805  }
806  else if (eltype == "hex")
807  {
808  type = Element::HEXAHEDRON;
809  }
810  else if (eltype == "wedge")
811  {
812  type = Element::WEDGE;
813  }
814  else if (eltype == "tet")
815  {
816  type = Element::TETRAHEDRON;
817  }
818  else
819  {
820  MFEM_ABORT("unrecognized element type (read '" << eltype
821  << "') in inline mesh format. "
822  "Allowed: segment, tri, quad, tet, hex, wedge");
823  }
824  }
825  else
826  {
827  MFEM_ABORT("unrecognized keyword (" << name
828  << ") in inline mesh format. "
829  "Allowed: nx, ny, nz, type, sx, sy, sz");
830  }
831 
832  input >> std::ws;
833  // Allow an optional semi-colon at the end of each line.
834  if (input.peek() == ';')
835  {
836  input.get();
837  }
838 
839  // Done reading file
840  if (!input)
841  {
842  break;
843  }
844  }
845 
846  // Now make the mesh.
847  if (type == Element::SEGMENT)
848  {
849  MFEM_VERIFY(nx > 0 && sx > 0.0,
850  "invalid 1D inline mesh format, all values must be "
851  "positive\n"
852  << " nx = " << nx << "\n"
853  << " sx = " << sx << "\n");
854  Make1D(nx, sx);
855  }
856  else if (type == Element::TRIANGLE || type == Element::QUADRILATERAL)
857  {
858  MFEM_VERIFY(nx > 0 && ny > 0 && sx > 0.0 && sy > 0.0,
859  "invalid 2D inline mesh format, all values must be "
860  "positive\n"
861  << " nx = " << nx << "\n"
862  << " ny = " << ny << "\n"
863  << " sx = " << sx << "\n"
864  << " sy = " << sy << "\n");
865  Make2D(nx, ny, type, sx, sy, generate_edges, true);
866  }
867  else if (type == Element::TETRAHEDRON || type == Element::WEDGE ||
868  type == Element::HEXAHEDRON)
869  {
870  MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
871  sx > 0.0 && sy > 0.0 && sz > 0.0,
872  "invalid 3D inline mesh format, all values must be "
873  "positive\n"
874  << " nx = " << nx << "\n"
875  << " ny = " << ny << "\n"
876  << " nz = " << nz << "\n"
877  << " sx = " << sx << "\n"
878  << " sy = " << sy << "\n"
879  << " sz = " << sz << "\n");
880  Make3D(nx, ny, nz, type, sx, sy, sz, true);
881  // TODO: maybe have an option in the file to control ordering?
882  }
883  else
884  {
885  MFEM_ABORT("For inline mesh, must specify an element type ="
886  " [segment, tri, quad, tet, hex, wedge]");
887  }
888 }
889 
890 void Mesh::ReadGmshMesh(std::istream &input)
891 {
892  string buff;
893  double version;
894  int binary, dsize;
895  input >> version >> binary >> dsize;
896  if (version < 2.2)
897  {
898  MFEM_ABORT("Gmsh file version < 2.2");
899  }
900  if (dsize != sizeof(double))
901  {
902  MFEM_ABORT("Gmsh file : dsize != sizeof(double)");
903  }
904  getline(input, buff);
905  // There is a number 1 in binary format
906  if (binary)
907  {
908  int one;
909  input.read(reinterpret_cast<char*>(&one), sizeof(one));
910  if (one != 1)
911  {
912  MFEM_ABORT("Gmsh file : wrong binary format");
913  }
914  }
915 
916  // A map between a serial number of the vertex and its number in the file
917  // (there may be gaps in the numbering, and also Gmsh enumerates vertices
918  // starting from 1, not 0)
919  map<int, int> vertices_map;
920  // Read the lines of the mesh file. If we face specific keyword, we'll treat
921  // the section.
922  while (input >> buff)
923  {
924  if (buff == "$Nodes") // reading mesh vertices
925  {
926  input >> NumOfVertices;
927  getline(input, buff);
928  vertices.SetSize(NumOfVertices);
929  int serial_number;
930  const int gmsh_dim = 3; // Gmsh always outputs 3 coordinates
931  double coord[gmsh_dim];
932  for (int ver = 0; ver < NumOfVertices; ++ver)
933  {
934  if (binary)
935  {
936  input.read(reinterpret_cast<char*>(&serial_number), sizeof(int));
937  input.read(reinterpret_cast<char*>(coord), gmsh_dim*sizeof(double));
938  }
939  else // ASCII
940  {
941  input >> serial_number;
942  for (int ci = 0; ci < gmsh_dim; ++ci)
943  {
944  input >> coord[ci];
945  }
946  }
947  vertices[ver] = Vertex(coord, gmsh_dim);
948  vertices_map[serial_number] = ver;
949  }
950  if (static_cast<int>(vertices_map.size()) != NumOfVertices)
951  {
952  MFEM_ABORT("Gmsh file : vertices indices are not unique");
953  }
954  } // section '$Nodes'
955  else if (buff == "$Elements") // reading mesh elements
956  {
957  int num_of_all_elements;
958  input >> num_of_all_elements;
959  // = NumOfElements + NumOfBdrElements + (maybe, PhysicalPoints)
960  getline(input, buff);
961 
962  int serial_number; // serial number of an element
963  int type_of_element; // ID describing a type of a mesh element
964  int n_tags; // number of different tags describing an element
965  int phys_domain; // element's attribute
966  int elem_domain; // another element's attribute (rarely used)
967  int n_partitions; // number of partitions where an element takes place
968 
969  // number of nodes for each type of Gmsh elements, type is the index of
970  // the array + 1
971  int nodes_of_gmsh_element[] =
972  {
973  2, // 2-node line.
974  3, // 3-node triangle.
975  4, // 4-node quadrangle.
976  4, // 4-node tetrahedron.
977  8, // 8-node hexahedron.
978  6, // 6-node prism.
979  5, // 5-node pyramid.
980  3, /* 3-node second order line (2 nodes associated with the vertices
981  and 1 with the edge). */
982  6, /* 6-node second order triangle (3 nodes associated with the
983  vertices and 3 with the edges). */
984  9, /* 9-node second order quadrangle (4 nodes associated with the
985  vertices, 4 with the edges and 1 with the face). */
986  10,/* 10-node second order tetrahedron (4 nodes associated with the
987  vertices and 6 with the edges). */
988  27,/* 27-node second order hexahedron (8 nodes associated with the
989  vertices, 12 with the edges, 6 with the faces and 1 with
990  the volume). */
991  18,/* 18-node second order prism (6 nodes associated with the
992  vertices, 9 with the edges and 3 with the quadrangular
993  faces). */
994  14,/* 14-node second order pyramid (5 nodes associated with the
995  vertices, 8 with the edges and 1 with the quadrangular
996  face). */
997  1, // 1-node point.
998  8, /* 8-node second order quadrangle (4 nodes associated with the
999  vertices and 4 with the edges). */
1000  20,/* 20-node second order hexahedron (8 nodes associated with the
1001  vertices and 12 with the edges). */
1002  15,/* 15-node second order prism (6 nodes associated with the
1003  vertices and 9 with the edges). */
1004  13,/* 13-node second order pyramid (5 nodes associated with the
1005  vertices and 8 with the edges). */
1006  9, /* 9-node third order incomplete triangle (3 nodes associated
1007  with the vertices, 6 with the edges) */
1008  10,/* 10-node third order triangle (3 nodes associated with the
1009  vertices, 6 with the edges, 1 with the face) */
1010  12,/* 12-node fourth order incomplete triangle (3 nodes associated
1011  with the vertices, 9 with the edges) */
1012  15,/* 15-node fourth order triangle (3 nodes associated with the
1013  vertices, 9 with the edges, 3 with the face) */
1014  15,/* 15-node fifth order incomplete triangle (3 nodes associated
1015  with the vertices, 12 with the edges) */
1016  21,/* 21-node fifth order complete triangle (3 nodes associated with
1017  the vertices, 12 with the edges, 6 with the face) */
1018  4, /* 4-node third order edge (2 nodes associated with the vertices,
1019  2 internal to the edge) */
1020  5, /* 5-node fourth order edge (2 nodes associated with the
1021  vertices, 3 internal to the edge) */
1022  6, /* 6-node fifth order edge (2 nodes associated with the vertices,
1023  4 internal to the edge) */
1024  20 /* 20-node third order tetrahedron (4 nodes associated with the
1025  vertices, 12 with the edges, 4 with the faces) */
1026  };
1027 
1028  vector<Element*> elements_0D, elements_1D, elements_2D, elements_3D;
1029  elements_0D.reserve(num_of_all_elements);
1030  elements_1D.reserve(num_of_all_elements);
1031  elements_2D.reserve(num_of_all_elements);
1032  elements_3D.reserve(num_of_all_elements);
1033 
1034  if (binary)
1035  {
1036  int n_elem_part = 0; // partial sum of elements that are read
1037  const int header_size = 3;
1038  // header consists of 3 numbers: type of the element, number of
1039  // elements of this type, and number of tags
1040  int header[header_size];
1041  int n_elem_one_type; // number of elements of a specific type
1042 
1043  while (n_elem_part < num_of_all_elements)
1044  {
1045  input.read(reinterpret_cast<char*>(header),
1046  header_size*sizeof(int));
1047  type_of_element = header[0];
1048  n_elem_one_type = header[1];
1049  n_tags = header[2];
1050 
1051  n_elem_part += n_elem_one_type;
1052 
1053  const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
1054  vector<int> data(1+n_tags+n_elem_nodes);
1055  for (int el = 0; el < n_elem_one_type; ++el)
1056  {
1057  input.read(reinterpret_cast<char*>(&data[0]),
1058  data.size()*sizeof(int));
1059  int dd = 0; // index for data array
1060  serial_number = data[dd++];
1061  // physical domain - the most important value (to distinguish
1062  // materials with different properties)
1063  phys_domain = (n_tags > 0) ? data[dd++] : 1;
1064  // elementary domain - to distinguish different geometrical
1065  // domains (typically, it's used rarely)
1066  elem_domain = (n_tags > 1) ? data[dd++] : 0;
1067  // the number of tags is bigger than 2 if there are some
1068  // partitions (domain decompositions)
1069  n_partitions = (n_tags > 2) ? data[dd++] : 0;
1070  // we currently just skip the partitions if they exist, and go
1071  // directly to vertices describing the mesh element
1072  vector<int> vert_indices(n_elem_nodes);
1073  for (int vi = 0; vi < n_elem_nodes; ++vi)
1074  {
1075  map<int, int>::const_iterator it =
1076  vertices_map.find(data[1+n_tags+vi]);
1077  if (it == vertices_map.end())
1078  {
1079  MFEM_ABORT("Gmsh file : vertex index doesn't exist");
1080  }
1081  vert_indices[vi] = it->second;
1082  }
1083 
1084  // non-positive attributes are not allowed in MFEM
1085  if (phys_domain <= 0)
1086  {
1087  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!");
1088  }
1089 
1090  // initialize the mesh element
1091  switch (type_of_element)
1092  {
1093  case 1: // 2-node line
1094  {
1095  elements_1D.push_back(
1096  new Segment(&vert_indices[0], phys_domain));
1097  break;
1098  }
1099  case 2: // 3-node triangle
1100  {
1101  elements_2D.push_back(
1102  new Triangle(&vert_indices[0], phys_domain));
1103  break;
1104  }
1105  case 3: // 4-node quadrangle
1106  {
1107  elements_2D.push_back(
1108  new Quadrilateral(&vert_indices[0], phys_domain));
1109  break;
1110  }
1111  case 4: // 4-node tetrahedron
1112  {
1113  elements_3D.push_back(
1114  new Tetrahedron(&vert_indices[0], phys_domain));
1115  break;
1116  }
1117  case 5: // 8-node hexahedron
1118  {
1119  elements_3D.push_back(
1120  new Hexahedron(&vert_indices[0], phys_domain));
1121  break;
1122  }
1123  case 15: // 1-node point
1124  {
1125  elements_0D.push_back(
1126  new Point(&vert_indices[0], phys_domain));
1127  break;
1128  }
1129  default: // any other element
1130  MFEM_WARNING("Unsupported Gmsh element type.");
1131  break;
1132 
1133  } // switch (type_of_element)
1134  } // el (elements of one type)
1135  } // all elements
1136  } // if binary
1137  else // ASCII
1138  {
1139  for (int el = 0; el < num_of_all_elements; ++el)
1140  {
1141  input >> serial_number >> type_of_element >> n_tags;
1142  vector<int> data(n_tags);
1143  for (int i = 0; i < n_tags; ++i) { input >> data[i]; }
1144  // physical domain - the most important value (to distinguish
1145  // materials with different properties)
1146  phys_domain = (n_tags > 0) ? data[0] : 1;
1147  // elementary domain - to distinguish different geometrical
1148  // domains (typically, it's used rarely)
1149  elem_domain = (n_tags > 1) ? data[1] : 0;
1150  // the number of tags is bigger than 2 if there are some
1151  // partitions (domain decompositions)
1152  n_partitions = (n_tags > 2) ? data[2] : 0;
1153  // we currently just skip the partitions if they exist, and go
1154  // directly to vertices describing the mesh element
1155  const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
1156  vector<int> vert_indices(n_elem_nodes);
1157  int index;
1158  for (int vi = 0; vi < n_elem_nodes; ++vi)
1159  {
1160  input >> index;
1161  map<int, int>::const_iterator it = vertices_map.find(index);
1162  if (it == vertices_map.end())
1163  {
1164  MFEM_ABORT("Gmsh file : vertex index doesn't exist");
1165  }
1166  vert_indices[vi] = it->second;
1167  }
1168 
1169  // non-positive attributes are not allowed in MFEM
1170  if (phys_domain <= 0)
1171  {
1172  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!");
1173  }
1174 
1175  // initialize the mesh element
1176  switch (type_of_element)
1177  {
1178  case 1: // 2-node line
1179  {
1180  elements_1D.push_back(
1181  new Segment(&vert_indices[0], phys_domain));
1182  break;
1183  }
1184  case 2: // 3-node triangle
1185  {
1186  elements_2D.push_back(
1187  new Triangle(&vert_indices[0], phys_domain));
1188  break;
1189  }
1190  case 3: // 4-node quadrangle
1191  {
1192  elements_2D.push_back(
1193  new Quadrilateral(&vert_indices[0], phys_domain));
1194  break;
1195  }
1196  case 4: // 4-node tetrahedron
1197  {
1198  elements_3D.push_back(
1199  new Tetrahedron(&vert_indices[0], phys_domain));
1200  break;
1201  }
1202  case 5: // 8-node hexahedron
1203  {
1204  elements_3D.push_back(
1205  new Hexahedron(&vert_indices[0], phys_domain));
1206  break;
1207  }
1208  case 15: // 1-node point
1209  {
1210  elements_0D.push_back(
1211  new Point(&vert_indices[0], phys_domain));
1212  break;
1213  }
1214  default: // any other element
1215  MFEM_WARNING("Unsupported Gmsh element type.");
1216  break;
1217 
1218  } // switch (type_of_element)
1219  } // el (all elements)
1220  } // if ASCII
1221 
1222  if (!elements_3D.empty())
1223  {
1224  Dim = 3;
1225  NumOfElements = elements_3D.size();
1226  elements.SetSize(NumOfElements);
1227  for (int el = 0; el < NumOfElements; ++el)
1228  {
1229  elements[el] = elements_3D[el];
1230  }
1231  NumOfBdrElements = elements_2D.size();
1232  boundary.SetSize(NumOfBdrElements);
1233  for (int el = 0; el < NumOfBdrElements; ++el)
1234  {
1235  boundary[el] = elements_2D[el];
1236  }
1237  // discard other elements
1238  for (size_t el = 0; el < elements_1D.size(); ++el)
1239  {
1240  delete elements_1D[el];
1241  }
1242  for (size_t el = 0; el < elements_0D.size(); ++el)
1243  {
1244  delete elements_0D[el];
1245  }
1246  }
1247  else if (!elements_2D.empty())
1248  {
1249  Dim = 2;
1250  NumOfElements = elements_2D.size();
1251  elements.SetSize(NumOfElements);
1252  for (int el = 0; el < NumOfElements; ++el)
1253  {
1254  elements[el] = elements_2D[el];
1255  }
1256  NumOfBdrElements = elements_1D.size();
1257  boundary.SetSize(NumOfBdrElements);
1258  for (int el = 0; el < NumOfBdrElements; ++el)
1259  {
1260  boundary[el] = elements_1D[el];
1261  }
1262  // discard other elements
1263  for (size_t el = 0; el < elements_0D.size(); ++el)
1264  {
1265  delete elements_0D[el];
1266  }
1267  }
1268  else if (!elements_1D.empty())
1269  {
1270  Dim = 1;
1271  NumOfElements = elements_1D.size();
1272  elements.SetSize(NumOfElements);
1273  for (int el = 0; el < NumOfElements; ++el)
1274  {
1275  elements[el] = elements_1D[el];
1276  }
1277  NumOfBdrElements = elements_0D.size();
1278  boundary.SetSize(NumOfBdrElements);
1279  for (int el = 0; el < NumOfBdrElements; ++el)
1280  {
1281  boundary[el] = elements_0D[el];
1282  }
1283  }
1284  else
1285  {
1286  MFEM_ABORT("Gmsh file : no elements found");
1287  return;
1288  }
1289 
1290  MFEM_CONTRACT_VAR(n_partitions);
1291  MFEM_CONTRACT_VAR(elem_domain);
1292 
1293  } // section '$Elements'
1294  } // we reach the end of the file
1295 }
1296 
1297 
1298 #ifdef MFEM_USE_NETCDF
1299 void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
1300 {
1301  read_gf = 0;
1302 
1303  // curved set to zero will change if mesh is indeed curved
1304  curved = 0;
1305 
1306  const int sideMapTri3[3][2] =
1307  {
1308  {1,2},
1309  {2,3},
1310  {3,1},
1311  };
1312 
1313  const int sideMapQuad4[4][2] =
1314  {
1315  {1,2},
1316  {2,3},
1317  {3,4},
1318  {4,1},
1319  };
1320 
1321  const int sideMapTri6[3][3] =
1322  {
1323  {1,2,4},
1324  {2,3,5},
1325  {3,1,6},
1326  };
1327 
1328  const int sideMapQuad9[4][3] =
1329  {
1330  {1,2,5},
1331  {2,3,6},
1332  {3,4,7},
1333  {4,1,8},
1334  };
1335 
1336  const int sideMapTet4[4][3] =
1337  {
1338  {1,2,4},
1339  {2,3,4},
1340  {1,4,3},
1341  {1,3,2}
1342  };
1343 
1344  const int sideMapTet10[4][6] =
1345  {
1346  {1,2,4,5,9,8},
1347  {2,3,4,6,10,9},
1348  {1,4,3,8,10,7},
1349  {1,3,2,7,6,5}
1350  };
1351 
1352  const int sideMapHex8[6][4] =
1353  {
1354  {1,2,6,5},
1355  {2,3,7,6},
1356  {4,3,7,8},
1357  {1,4,8,5},
1358  {1,4,3,2},
1359  {5,8,7,6}
1360  };
1361 
1362  const int sideMapHex27[6][9] =
1363  {
1364  {1,2,6,5,9,14,17,13,26},
1365  {2,3,7,6,10,15,18,14,25},
1366  {4,3,7,8,11,15,19,16,27},
1367  {1,4,8,5,12,16,20,13,24},
1368  {1,4,3,2,12,11,10,9,22},
1369  {5,8,7,6,20,19,18,17,23}
1370  };
1371 
1372 
1373  // 1,2,3,4,5,6,7,8,9,10
1374  const int mfemToGenesisTet10[10] = {1,2,3,4,5,7,8,6,9,10};
1375 
1376  // 1,2,3,4,5,6,7,8,9,10,11,
1377  const int mfemToGenesisHex27[27] = {1,2,3,4,5,6,7,8,9,10,11,
1378  // 12,13,14,15,16,17,18,19
1379  12,17,18,19,20,13,14,15,
1380  // 20,21,22,23,24,25,26,27
1381  16,22,26,25,27,24,23,21
1382  };
1383 
1384  const int mfemToGenesisTri6[6] = {1,2,3,4,5,6};
1385  const int mfemToGenesisQuad9[9] = {1,2,3,4,5,6,7,8,9};
1386 
1387 
1388  // error handling.
1389  int retval;
1390 
1391  // dummy string
1392  char str_dummy[256];
1393 
1394  char temp_str[256];
1395  int temp_id;
1396 
1397  // open the file.
1398  int ncid;
1399  if ((retval = nc_open(filename, NC_NOWRITE, &ncid)))
1400  {
1401  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1402  }
1403 
1404  // read important dimensions
1405 
1406  int id;
1407  size_t num_dim=0, num_nodes=0, num_elem=0, num_el_blk=0, num_side_sets=0;
1408 
1409  if ((retval = nc_inq_dimid(ncid, "num_dim", &id)) ||
1410  (retval = nc_inq_dim(ncid, id, str_dummy, &num_dim)) ||
1411 
1412  (retval = nc_inq_dimid(ncid, "num_nodes", &id)) ||
1413  (retval = nc_inq_dim(ncid, id, str_dummy, &num_nodes)) ||
1414 
1415  (retval = nc_inq_dimid(ncid, "num_elem", &id)) ||
1416  (retval = nc_inq_dim(ncid, id, str_dummy, &num_elem)) ||
1417 
1418  (retval = nc_inq_dimid(ncid, "num_el_blk", &id)) ||
1419  (retval = nc_inq_dim(ncid, id, str_dummy, &num_el_blk)))
1420  {
1421  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1422  }
1423  if ((retval = nc_inq_dimid(ncid, "num_side_sets", &id)) ||
1424  (retval = nc_inq_dim(ncid, id, str_dummy, &num_side_sets)))
1425  {
1426  num_side_sets = 0;
1427  }
1428 
1429  Dim = num_dim;
1430 
1431  // create arrays for element blocks
1432  size_t *num_el_in_blk = new size_t[num_el_blk];
1433  size_t num_node_per_el;
1434 
1435  int previous_num_node_per_el = 0;
1436  for (int i = 0; i < (int) num_el_blk; i++)
1437  {
1438  sprintf(temp_str, "num_el_in_blk%d", i+1);
1439  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1440  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_el_in_blk[i])))
1441  {
1442  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1443  }
1444 
1445  sprintf(temp_str, "num_nod_per_el%d", i+1);
1446  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1447  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_node_per_el)))
1448  {
1449  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1450  }
1451 
1452  // check for different element types in each block
1453  // which is not currently supported
1454  if (i != 0)
1455  {
1456  if ((int) num_node_per_el != previous_num_node_per_el)
1457  {
1458  MFEM_ABORT("Element blocks of different element types not supported");
1459  }
1460  }
1461  previous_num_node_per_el = num_node_per_el;
1462  }
1463 
1464  // Determine CUBIT element and face type
1465  enum CubitElementType
1466  {
1467  ELEMENT_TRI3,
1468  ELEMENT_TRI6,
1469  ELEMENT_QUAD4,
1470  ELEMENT_QUAD9,
1471  ELEMENT_TET4,
1472  ELEMENT_TET10,
1473  ELEMENT_HEX8,
1474  ELEMENT_HEX27
1475  };
1476 
1477  enum CubitFaceType
1478  {
1479  FACE_EDGE2,
1480  FACE_EDGE3,
1481  FACE_TRI3,
1482  FACE_TRI6,
1483  FACE_QUAD4,
1484  FACE_QUAD9
1485  };
1486 
1487  CubitElementType cubit_element_type = ELEMENT_TRI3; // suppress a warning
1488  CubitFaceType cubit_face_type = FACE_EDGE2; // suppress a warning
1489  int num_element_linear_nodes = 0; // initialize to suppress a warning
1490 
1491  if (num_dim == 2)
1492  {
1493  switch (num_node_per_el)
1494  {
1495  case (3) :
1496  {
1497  cubit_element_type = ELEMENT_TRI3;
1498  cubit_face_type = FACE_EDGE2;
1499  num_element_linear_nodes = 3;
1500  break;
1501  }
1502  case (6) :
1503  {
1504  cubit_element_type = ELEMENT_TRI6;
1505  cubit_face_type = FACE_EDGE3;
1506  num_element_linear_nodes = 3;
1507  break;
1508  }
1509  case (4) :
1510  {
1511  cubit_element_type = ELEMENT_QUAD4;
1512  cubit_face_type = FACE_EDGE2;
1513  num_element_linear_nodes = 4;
1514  break;
1515  }
1516  case (9) :
1517  {
1518  cubit_element_type = ELEMENT_QUAD9;
1519  cubit_face_type = FACE_EDGE3;
1520  num_element_linear_nodes = 4;
1521  break;
1522  }
1523  default :
1524  {
1525  MFEM_ABORT("Don't know what to do with a " << num_node_per_el <<
1526  " node 2D element\n");
1527  }
1528  }
1529  }
1530  else if (num_dim == 3)
1531  {
1532  switch (num_node_per_el)
1533  {
1534  case (4) :
1535  {
1536  cubit_element_type = ELEMENT_TET4;
1537  cubit_face_type = FACE_TRI3;
1538  num_element_linear_nodes = 4;
1539  break;
1540  }
1541  case (10) :
1542  {
1543  cubit_element_type = ELEMENT_TET10;
1544  cubit_face_type = FACE_TRI6;
1545  num_element_linear_nodes = 4;
1546  break;
1547  }
1548  case (8) :
1549  {
1550  cubit_element_type = ELEMENT_HEX8;
1551  cubit_face_type = FACE_QUAD4;
1552  num_element_linear_nodes = 8;
1553  break;
1554  }
1555  case (27) :
1556  {
1557  cubit_element_type = ELEMENT_HEX27;
1558  cubit_face_type = FACE_QUAD9;
1559  num_element_linear_nodes = 8;
1560  break;
1561  }
1562  default :
1563  {
1564  MFEM_ABORT("Don't know what to do with a " << num_node_per_el <<
1565  " node 3D element\n");
1566  }
1567  }
1568  }
1569  else
1570  {
1571  MFEM_ABORT("Invalid dimension: num_dim = " << num_dim);
1572  }
1573 
1574  // Determine order of elements
1575  int order = 0;
1576  if (cubit_element_type == ELEMENT_TRI3 || cubit_element_type == ELEMENT_QUAD4 ||
1577  cubit_element_type == ELEMENT_TET4 || cubit_element_type == ELEMENT_HEX8)
1578  {
1579  order = 1;
1580  }
1581  else if (cubit_element_type == ELEMENT_TRI6 ||
1582  cubit_element_type == ELEMENT_QUAD9 ||
1583  cubit_element_type == ELEMENT_TET10 || cubit_element_type == ELEMENT_HEX27)
1584  {
1585  order = 2;
1586  }
1587 
1588  // create array for number of sides in side sets
1589  size_t *num_side_in_ss = new size_t[num_side_sets];
1590  for (int i = 0; i < (int) num_side_sets; i++)
1591  {
1592  sprintf(temp_str, "num_side_ss%d", i+1);
1593  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1594  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_side_in_ss[i])))
1595  {
1596  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1597  }
1598  }
1599 
1600  // read the coordinates
1601  double *coordx = new double[num_nodes];
1602  double *coordy = new double[num_nodes];
1603  double *coordz = new double[num_nodes];
1604 
1605  if ((retval = nc_inq_varid(ncid, "coordx", &id)) ||
1606  (retval = nc_get_var_double(ncid, id, coordx)) ||
1607  (retval = nc_inq_varid(ncid, "coordy", &id)) ||
1608  (retval = nc_get_var_double(ncid, id, coordy)))
1609  {
1610  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1611  }
1612 
1613  if (num_dim == 3)
1614  {
1615  if ((retval = nc_inq_varid(ncid, "coordz", &id)) ||
1616  (retval = nc_get_var_double(ncid, id, coordz)))
1617  {
1618  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1619  }
1620  }
1621 
1622  // read the element blocks
1623  int **elem_blk = new int*[num_el_blk];
1624  for (int i = 0; i < (int) num_el_blk; i++)
1625  {
1626  elem_blk[i] = new int[num_el_in_blk[i] * num_node_per_el];
1627  sprintf(temp_str, "connect%d", i+1);
1628  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1629  (retval = nc_get_var_int(ncid, temp_id, elem_blk[i])))
1630  {
1631  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1632  }
1633  }
1634  int *ebprop = new int[num_el_blk];
1635  if ((retval = nc_inq_varid(ncid, "eb_prop1", &id)) ||
1636  (retval = nc_get_var_int(ncid, id, ebprop)))
1637  {
1638  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1639  }
1640 
1641  // read the side sets, a side is is given by (element, face) pairs
1642 
1643  int **elem_ss = new int*[num_side_sets];
1644  int **side_ss = new int*[num_side_sets];
1645 
1646  for (int i = 0; i < (int) num_side_sets; i++)
1647  {
1648  elem_ss[i] = new int[num_side_in_ss[i]];
1649  side_ss[i] = new int[num_side_in_ss[i]];
1650 
1651  sprintf(temp_str, "elem_ss%d", i+1);
1652  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1653  (retval = nc_get_var_int(ncid, temp_id, elem_ss[i])))
1654  {
1655  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1656  }
1657 
1658  sprintf(temp_str,"side_ss%d",i+1);
1659  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1660  (retval = nc_get_var_int(ncid, temp_id, side_ss[i])))
1661  {
1662  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1663  }
1664  }
1665 
1666  int *ssprop = new int[num_side_sets];
1667  if ((num_side_sets > 0) &&
1668  ((retval = nc_inq_varid(ncid, "ss_prop1", &id)) ||
1669  (retval = nc_get_var_int(ncid, id, ssprop))))
1670  {
1671  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
1672  }
1673 
1674  // convert (elem,side) pairs to 2D elements
1675 
1676 
1677  int num_face_nodes = 0;
1678  int num_face_linear_nodes = 0;
1679 
1680  switch (cubit_face_type)
1681  {
1682  case (FACE_EDGE2):
1683  {
1684  num_face_nodes = 2;
1685  num_face_linear_nodes = 2;
1686  break;
1687  }
1688  case (FACE_EDGE3):
1689  {
1690  num_face_nodes = 3;
1691  num_face_linear_nodes = 2;
1692  break;
1693  }
1694  case (FACE_TRI3):
1695  {
1696  num_face_nodes = 3;
1697  num_face_linear_nodes = 3;
1698  break;
1699  }
1700  case (FACE_TRI6):
1701  {
1702  num_face_nodes = 6;
1703  num_face_linear_nodes = 3;
1704  break;
1705  }
1706  case (FACE_QUAD4):
1707  {
1708  num_face_nodes = 4;
1709  num_face_linear_nodes = 4;
1710  break;
1711  }
1712  case (FACE_QUAD9):
1713  {
1714  num_face_nodes = 9;
1715  num_face_linear_nodes = 4;
1716  break;
1717  }
1718  }
1719 
1720  // given a global element number, determine the element block and local
1721  // element number
1722  int *start_of_block = new int[num_el_blk+1];
1723  start_of_block[0] = 0;
1724  for (int i = 1; i < (int) num_el_blk+1; i++)
1725  {
1726  start_of_block[i] = start_of_block[i-1] + num_el_in_blk[i-1];
1727  }
1728 
1729  int **ss_node_id = new int*[num_side_sets];
1730 
1731  for (int i = 0; i < (int) num_side_sets; i++)
1732  {
1733  ss_node_id[i] = new int[num_side_in_ss[i]*num_face_nodes];
1734  for (int j = 0; j < (int) num_side_in_ss[i]; j++)
1735  {
1736  int glob_ind = elem_ss[i][j]-1;
1737  int iblk = 0;
1738  int loc_ind;
1739  while (iblk < (int) num_el_blk && glob_ind >= start_of_block[iblk+1])
1740  {
1741  iblk++;
1742  }
1743  if (iblk >= (int) num_el_blk)
1744  {
1745  MFEM_ABORT("Sideset element does not exist");
1746  }
1747  loc_ind = glob_ind - start_of_block[iblk];
1748  int this_side = side_ss[i][j];
1749  int ielem = loc_ind*num_node_per_el;
1750 
1751  for (int k = 0; k < num_face_nodes; k++)
1752  {
1753  int inode;
1754  switch (cubit_element_type)
1755  {
1756  case (ELEMENT_TRI3):
1757  {
1758  inode = sideMapTri3[this_side-1][k];
1759  break;
1760  }
1761  case (ELEMENT_TRI6):
1762  {
1763  inode = sideMapTri6[this_side-1][k];
1764  break;
1765  }
1766  case (ELEMENT_QUAD4):
1767  {
1768  inode = sideMapQuad4[this_side-1][k];
1769  break;
1770  }
1771  case (ELEMENT_QUAD9):
1772  {
1773  inode = sideMapQuad9[this_side-1][k];
1774  break;
1775  }
1776  case (ELEMENT_TET4):
1777  {
1778  inode = sideMapTet4[this_side-1][k];
1779  break;
1780  }
1781  case (ELEMENT_TET10):
1782  {
1783  inode = sideMapTet10[this_side-1][k];
1784  break;
1785  }
1786  case (ELEMENT_HEX8):
1787  {
1788  inode = sideMapHex8[this_side-1][k];
1789  break;
1790  }
1791  case (ELEMENT_HEX27):
1792  {
1793  inode = sideMapHex27[this_side-1][k];
1794  break;
1795  }
1796  }
1797  ss_node_id[i][j*num_face_nodes+k] =
1798  elem_blk[iblk][ielem + inode - 1];
1799  }
1800  }
1801  }
1802 
1803  // we need another node ID mapping since MFEM needs contiguous vertex IDs
1804  std::vector<int> uniqueVertexID;
1805 
1806  for (int iblk = 0; iblk < (int) num_el_blk; iblk++)
1807  {
1808  for (int i = 0; i < (int) num_el_in_blk[iblk]; i++)
1809  {
1810  for (int j = 0; j < num_element_linear_nodes; j++)
1811  {
1812  uniqueVertexID.push_back(elem_blk[iblk][i*num_node_per_el + j]);
1813  }
1814  }
1815  }
1816  std::sort(uniqueVertexID.begin(), uniqueVertexID.end());
1817  std::vector<int>::iterator newEnd;
1818  newEnd = std::unique(uniqueVertexID.begin(), uniqueVertexID.end());
1819  uniqueVertexID.resize(std::distance(uniqueVertexID.begin(), newEnd));
1820 
1821  // OK at this point uniqueVertexID contains a list of all the nodes that are
1822  // actually used by the mesh, 1-based, and sorted. We need to invert this
1823  // list, the inverse is a map
1824 
1825  std::map<int,int> cubitToMFEMVertMap;
1826  for (int i = 0; i < (int) uniqueVertexID.size(); i++)
1827  {
1828  cubitToMFEMVertMap[uniqueVertexID[i]] = i+1;
1829  }
1830  MFEM_ASSERT(cubitToMFEMVertMap.size() == uniqueVertexID.size(),
1831  "This should never happen\n");
1832 
1833  // OK now load up the MFEM mesh structures
1834 
1835  // load up the vertices
1836 
1837  NumOfVertices = uniqueVertexID.size();
1838  vertices.SetSize(NumOfVertices);
1839  for (int i = 0; i < (int) uniqueVertexID.size(); i++)
1840  {
1841  vertices[i](0) = coordx[uniqueVertexID[i] - 1];
1842  vertices[i](1) = coordy[uniqueVertexID[i] - 1];
1843  if (Dim == 3)
1844  {
1845  vertices[i](2) = coordz[uniqueVertexID[i] - 1];
1846  }
1847  }
1848 
1849  NumOfElements = num_elem;
1850  elements.SetSize(num_elem);
1851  int elcount = 0;
1852  int renumberedVertID[8];
1853  for (int iblk = 0; iblk < (int) num_el_blk; iblk++)
1854  {
1855  int NumNodePerEl = num_node_per_el;
1856  for (int i = 0; i < (int) num_el_in_blk[iblk]; i++)
1857  {
1858  for (int j = 0; j < num_element_linear_nodes; j++)
1859  {
1860  renumberedVertID[j] =
1861  cubitToMFEMVertMap[elem_blk[iblk][i*NumNodePerEl+j]]-1;
1862  }
1863 
1864  switch (cubit_element_type)
1865  {
1866  case (ELEMENT_TRI3):
1867  case (ELEMENT_TRI6):
1868  {
1869  elements[elcount] = new Triangle(renumberedVertID,ebprop[iblk]);
1870  break;
1871  }
1872  case (ELEMENT_QUAD4):
1873  case (ELEMENT_QUAD9):
1874  {
1875  elements[elcount] = new Quadrilateral(renumberedVertID,ebprop[iblk]);
1876  break;
1877  }
1878  case (ELEMENT_TET4):
1879  case (ELEMENT_TET10):
1880  {
1881  elements[elcount] = new Tetrahedron(renumberedVertID,ebprop[iblk]);
1882  break;
1883  }
1884  case (ELEMENT_HEX8):
1885  case (ELEMENT_HEX27):
1886  {
1887  elements[elcount] = new Hexahedron(renumberedVertID,ebprop[iblk]);
1888  break;
1889  }
1890  }
1891  elcount++;
1892  }
1893  }
1894 
1895  // load up the boundary elements
1896 
1897  NumOfBdrElements = 0;
1898  for (int iss = 0; iss < (int) num_side_sets; iss++)
1899  {
1900  NumOfBdrElements += num_side_in_ss[iss];
1901  }
1902  boundary.SetSize(NumOfBdrElements);
1903  int sidecount = 0;
1904  for (int iss = 0; iss < (int) num_side_sets; iss++)
1905  {
1906  for (int i = 0; i < (int) num_side_in_ss[iss]; i++)
1907  {
1908  for (int j = 0; j < num_face_linear_nodes; j++)
1909  {
1910  renumberedVertID[j] =
1911  cubitToMFEMVertMap[ss_node_id[iss][i*num_face_nodes+j]] - 1;
1912  }
1913  switch (cubit_face_type)
1914  {
1915  case (FACE_EDGE2):
1916  case (FACE_EDGE3):
1917  {
1918  boundary[sidecount] = new Segment(renumberedVertID,ssprop[iss]);
1919  break;
1920  }
1921  case (FACE_TRI3):
1922  case (FACE_TRI6):
1923  {
1924  boundary[sidecount] = new Triangle(renumberedVertID,ssprop[iss]);
1925  break;
1926  }
1927  case (FACE_QUAD4):
1928  case (FACE_QUAD9):
1929  {
1930  boundary[sidecount] = new Quadrilateral(renumberedVertID,ssprop[iss]);
1931  break;
1932  }
1933  }
1934  sidecount++;
1935  }
1936  }
1937 
1938  if (order == 2)
1939  {
1940  curved = 1;
1941  int *mymap = NULL;
1942 
1943  switch (cubit_element_type)
1944  {
1945  case (ELEMENT_TRI6):
1946  {
1947  mymap = (int *) mfemToGenesisTri6;
1948  break;
1949  }
1950  case (ELEMENT_QUAD9):
1951  {
1952  mymap = (int *) mfemToGenesisQuad9;
1953  break;
1954  }
1955  case (ELEMENT_TET10):
1956  {
1957  mymap = (int *) mfemToGenesisTet10;
1958  break;
1959  }
1960  case (ELEMENT_HEX27):
1961  {
1962  mymap = (int *) mfemToGenesisHex27;
1963  break;
1964  }
1965  case (ELEMENT_TRI3):
1966  case (ELEMENT_QUAD4):
1967  case (ELEMENT_TET4):
1968  case (ELEMENT_HEX8):
1969  {
1970  MFEM_ABORT("Something went wrong. Linear elements detected when order is 2.");
1971  break;
1972  }
1973  }
1974 
1975  FinalizeTopology();
1976 
1977  // Define quadratic FE space
1978  FiniteElementCollection *fec = new H1_FECollection(2,3);
1979  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
1980  Ordering::byVDIM);
1981  Nodes = new GridFunction(fes);
1982  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
1983  own_nodes = 1;
1984 
1985  // int nTotDofs = fes->GetNDofs();
1986  // int nTotVDofs = fes->GetVSize();
1987  // mfem::out << endl << "nTotDofs = " << nTotDofs << " nTotVDofs "
1988  // << nTotVDofs << endl << endl;
1989 
1990  for (int i = 0; i < NumOfElements; i++)
1991  {
1992  Array<int> dofs;
1993 
1994  fes->GetElementDofs(i, dofs);
1995  Array<int> vdofs;
1996  vdofs.SetSize(dofs.Size());
1997  for (int l = 0; l < dofs.Size(); l++) { vdofs[l] = dofs[l]; }
1998  fes->DofsToVDofs(vdofs);
1999  int iblk = 0;
2000  int loc_ind;
2001  while (iblk < (int) num_el_blk && i >= start_of_block[iblk+1]) { iblk++; }
2002  loc_ind = i - start_of_block[iblk];
2003  for (int j = 0; j < dofs.Size(); j++)
2004  {
2005  int point_id = elem_blk[iblk][loc_ind*num_node_per_el + mymap[j] - 1] - 1;
2006  (*Nodes)(vdofs[j]) = coordx[point_id];
2007  (*Nodes)(vdofs[j]+1) = coordy[point_id];
2008  if (Dim == 3)
2009  {
2010  (*Nodes)(vdofs[j]+2) = coordz[point_id];
2011  }
2012  }
2013  }
2014  }
2015 
2016  // clean up all netcdf stuff
2017 
2018  nc_close(ncid);
2019 
2020  for (int i = 0; i < (int) num_side_sets; i++)
2021  {
2022  delete [] elem_ss[i];
2023  delete [] side_ss[i];
2024  }
2025 
2026  delete [] elem_ss;
2027  delete [] side_ss;
2028  delete [] num_el_in_blk;
2029  delete [] num_side_in_ss;
2030  delete [] coordx;
2031  delete [] coordy;
2032  delete [] coordz;
2033 
2034  for (int i = 0; i < (int) num_el_blk; i++)
2035  {
2036  delete [] elem_blk[i];
2037  }
2038 
2039  delete [] elem_blk;
2040  delete [] start_of_block;
2041 
2042  for (int i = 0; i < (int) num_side_sets; i++)
2043  {
2044  delete [] ss_node_id[i];
2045  }
2046  delete [] ss_node_id;
2047  delete [] ebprop;
2048  delete [] ssprop;
2049 
2050 }
2051 #endif // #ifdef MFEM_USE_NETCDF
2052 
2053 } // namespace mfem
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:289
int Size() const
Logical size of the array.
Definition: array.hpp:124
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
Data type for vertex.
Definition: vertex.hpp:21
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
Data type quadrilateral element.
Data type Wedge element.
Definition: wedge.hpp:22
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:27
void DeleteAll()
Delete whole array.
Definition: array.hpp:802
Data type hexahedron element.
Definition: hexahedron.hpp:22
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:1540
Data type triangle element.
Definition: triangle.hpp:23
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
Definition: ncmesh.hpp:100
Data type tetrahedron element.
Definition: tetrahedron.hpp:22
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
double a
Definition: lissajous.cpp:41
void Destroy()
Destroy a vector.
Definition: vector.hpp:478
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:368
void filter_dos(std::string &line)
Definition: text.hpp:41
Vector data type.
Definition: vector.hpp:48
Data type point element.
Definition: point.hpp:22
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Data type line segment element.
Definition: segment.hpp:22
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:107