MFEM  v4.5.2
Finite element discretization library
mesh_readers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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/binaryio.hpp"
15 #include "../general/text.hpp"
16 #include "../general/tinyxml2.h"
17 #include "gmsh.hpp"
18 
19 #include <iostream>
20 #include <cstdio>
21 #include <vector>
22 #include <algorithm>
23 
24 #ifdef MFEM_USE_NETCDF
25 #include "netcdf.h"
26 #endif
27 
28 #ifdef MFEM_USE_ZLIB
29 #include <zlib.h>
30 #endif
31 
32 using namespace std;
33 
34 namespace mfem
35 {
36 
37 bool Mesh::remove_unused_vertices = true;
38 
39 void Mesh::ReadMFEMMesh(std::istream &input, int version, int &curved)
40 {
41  // Read MFEM mesh v1.0 or v1.2 format
42  MFEM_VERIFY(version == 10 || version == 12,
43  "unknown MFEM mesh version");
44 
45  string ident;
46 
47  // read lines beginning with '#' (comments)
48  skip_comment_lines(input, '#');
49  input >> ident; // 'dimension'
50 
51  MFEM_VERIFY(ident == "dimension", "invalid mesh file");
52  input >> Dim;
53 
54  skip_comment_lines(input, '#');
55  input >> ident; // 'elements'
56 
57  MFEM_VERIFY(ident == "elements", "invalid mesh file");
58  input >> NumOfElements;
59  elements.SetSize(NumOfElements);
60  for (int j = 0; j < NumOfElements; j++)
61  {
62  elements[j] = ReadElement(input);
63  }
64 
65  skip_comment_lines(input, '#');
66  input >> ident; // 'boundary'
67 
68  MFEM_VERIFY(ident == "boundary", "invalid mesh file");
69  input >> NumOfBdrElements;
70  boundary.SetSize(NumOfBdrElements);
71  for (int j = 0; j < NumOfBdrElements; j++)
72  {
73  boundary[j] = ReadElement(input);
74  }
75 
76  skip_comment_lines(input, '#');
77  input >> ident; // 'vertices'
78 
79  MFEM_VERIFY(ident == "vertices", "invalid mesh file");
80  input >> NumOfVertices;
81  vertices.SetSize(NumOfVertices);
82 
83  input >> ws >> ident;
84  if (ident != "nodes")
85  {
86  // read the vertices
87  spaceDim = atoi(ident.c_str());
88  for (int j = 0; j < NumOfVertices; j++)
89  {
90  for (int i = 0; i < spaceDim; i++)
91  {
92  input >> vertices[j](i);
93  }
94  }
95  }
96  else
97  {
98  // prepare to read the nodes
99  input >> ws;
100  curved = 1;
101  }
102 
103  // When visualizing solutions on non-conforming grids, PETSc
104  // may dump additional vertices
105  if (remove_unused_vertices) { RemoveUnusedVertices(); }
106 }
107 
108 void Mesh::ReadLineMesh(std::istream &input)
109 {
110  int j,p1,p2,a;
111 
112  Dim = 1;
113 
114  input >> NumOfVertices;
115  vertices.SetSize(NumOfVertices);
116  // Sets vertices and the corresponding coordinates
117  for (j = 0; j < NumOfVertices; j++)
118  {
119  input >> vertices[j](0);
120  }
121 
122  input >> NumOfElements;
123  elements.SetSize(NumOfElements);
124  // Sets elements and the corresponding indices of vertices
125  for (j = 0; j < NumOfElements; j++)
126  {
127  input >> a >> p1 >> p2;
128  elements[j] = new Segment(p1-1, p2-1, a);
129  }
130 
131  int ind[1];
132  input >> NumOfBdrElements;
133  boundary.SetSize(NumOfBdrElements);
134  for (j = 0; j < NumOfBdrElements; j++)
135  {
136  input >> a >> ind[0];
137  ind[0]--;
138  boundary[j] = new Point(ind,a);
139  }
140 }
141 
142 void Mesh::ReadNetgen2DMesh(std::istream &input, int &curved)
143 {
144  int ints[32], attr, n;
145 
146  // Read planar mesh in Netgen format.
147  Dim = 2;
148 
149  // Read the boundary elements.
150  input >> NumOfBdrElements;
151  boundary.SetSize(NumOfBdrElements);
152  for (int i = 0; i < NumOfBdrElements; i++)
153  {
154  input >> attr
155  >> ints[0] >> ints[1];
156  ints[0]--; ints[1]--;
157  boundary[i] = new Segment(ints, attr);
158  }
159 
160  // Read the elements.
161  input >> NumOfElements;
162  elements.SetSize(NumOfElements);
163  for (int i = 0; i < NumOfElements; i++)
164  {
165  input >> attr >> n;
166  for (int j = 0; j < n; j++)
167  {
168  input >> ints[j];
169  ints[j]--;
170  }
171  switch (n)
172  {
173  case 2:
174  elements[i] = new Segment(ints, attr);
175  break;
176  case 3:
177  elements[i] = new Triangle(ints, attr);
178  break;
179  case 4:
180  elements[i] = new Quadrilateral(ints, attr);
181  break;
182  }
183  }
184 
185  if (!curved)
186  {
187  // Read the vertices.
188  input >> NumOfVertices;
189  vertices.SetSize(NumOfVertices);
190  for (int i = 0; i < NumOfVertices; i++)
191  for (int j = 0; j < Dim; j++)
192  {
193  input >> vertices[i](j);
194  }
195  }
196  else
197  {
198  input >> NumOfVertices;
199  vertices.SetSize(NumOfVertices);
200  input >> ws;
201  }
202 }
203 
204 void Mesh::ReadNetgen3DMesh(std::istream &input)
205 {
206  int ints[32], attr;
207 
208  // Read a Netgen format mesh of tetrahedra.
209  Dim = 3;
210 
211  // Read the vertices
212  input >> NumOfVertices;
213 
214  vertices.SetSize(NumOfVertices);
215  for (int i = 0; i < NumOfVertices; i++)
216  for (int j = 0; j < Dim; j++)
217  {
218  input >> vertices[i](j);
219  }
220 
221  // Read the elements
222  input >> NumOfElements;
223  elements.SetSize(NumOfElements);
224  for (int i = 0; i < NumOfElements; i++)
225  {
226  input >> attr;
227  for (int j = 0; j < 4; j++)
228  {
229  input >> ints[j];
230  ints[j]--;
231  }
232 #ifdef MFEM_USE_MEMALLOC
233  Tetrahedron *tet;
234  tet = TetMemory.Alloc();
235  tet->SetVertices(ints);
236  tet->SetAttribute(attr);
237  elements[i] = tet;
238 #else
239  elements[i] = new Tetrahedron(ints, attr);
240 #endif
241  }
242 
243  // Read the boundary information.
244  input >> NumOfBdrElements;
245  boundary.SetSize(NumOfBdrElements);
246  for (int i = 0; i < NumOfBdrElements; i++)
247  {
248  input >> attr;
249  for (int j = 0; j < 3; j++)
250  {
251  input >> ints[j];
252  ints[j]--;
253  }
254  boundary[i] = new Triangle(ints, attr);
255  }
256 }
257 
258 void Mesh::ReadTrueGridMesh(std::istream &input)
259 {
260  int i, j, ints[32], attr;
261  const int buflen = 1024;
262  char buf[buflen];
263 
264  // TODO: find the actual dimension
265  Dim = 3;
266 
267  if (Dim == 2)
268  {
269  int vari;
270  double varf;
271 
272  input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
273  input.getline(buf, buflen);
274  input.getline(buf, buflen);
275  input >> vari;
276  input.getline(buf, buflen);
277  input.getline(buf, buflen);
278  input.getline(buf, buflen);
279 
280  // Read the vertices.
281  vertices.SetSize(NumOfVertices);
282  for (i = 0; i < NumOfVertices; i++)
283  {
284  input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
285  input.getline(buf, buflen);
286  }
287 
288  // Read the elements.
289  elements.SetSize(NumOfElements);
290  for (i = 0; i < NumOfElements; i++)
291  {
292  input >> vari >> attr;
293  for (j = 0; j < 4; j++)
294  {
295  input >> ints[j];
296  ints[j]--;
297  }
298  input.getline(buf, buflen);
299  input.getline(buf, buflen);
300  elements[i] = new Quadrilateral(ints, attr);
301  }
302  }
303  else if (Dim == 3)
304  {
305  int vari;
306  double varf;
307  input >> vari >> NumOfVertices >> NumOfElements;
308  input.getline(buf, buflen);
309  input.getline(buf, buflen);
310  input >> vari >> vari >> NumOfBdrElements;
311  input.getline(buf, buflen);
312  input.getline(buf, buflen);
313  input.getline(buf, buflen);
314  // Read the vertices.
315  vertices.SetSize(NumOfVertices);
316  for (i = 0; i < NumOfVertices; i++)
317  {
318  input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
319  >> vertices[i](2);
320  input.getline(buf, buflen);
321  }
322  // Read the elements.
323  elements.SetSize(NumOfElements);
324  for (i = 0; i < NumOfElements; i++)
325  {
326  input >> vari >> attr;
327  for (j = 0; j < 8; j++)
328  {
329  input >> ints[j];
330  ints[j]--;
331  }
332  input.getline(buf, buflen);
333  elements[i] = new Hexahedron(ints, attr);
334  }
335  // Read the boundary elements.
336  boundary.SetSize(NumOfBdrElements);
337  for (i = 0; i < NumOfBdrElements; i++)
338  {
339  input >> attr;
340  for (j = 0; j < 4; j++)
341  {
342  input >> ints[j];
343  ints[j]--;
344  }
345  input.getline(buf, buflen);
346  boundary[i] = new Quadrilateral(ints, attr);
347  }
348  }
349 }
350 
351 // see Tetrahedron::edges
352 const int Mesh::vtk_quadratic_tet[10] =
353 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
354 
355 // see Pyramid::edges & Mesh::GenerateFaces
356 // https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
357 const int Mesh::vtk_quadratic_pyramid[13] =
358 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
359 
360 // see Wedge::edges & Mesh::GenerateFaces
361 // https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
362 const int Mesh::vtk_quadratic_wedge[18] =
363 { 0, 2, 1, 3, 5, 4, 8, 7, 6, 11, 10, 9, 12, 14, 13, 17, 16, 15};
364 
365 // see Hexahedron::edges & Mesh::GenerateFaces
366 const int Mesh::vtk_quadratic_hex[27] =
367 {
368  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
369  24, 22, 21, 23, 20, 25, 26
370 };
371 
372 void Mesh::CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
373  const Array<int> &cell_offsets,
374  const Array<int> &cell_types,
375  const Array<int> &cell_attributes,
376  int &curved, int &read_gf, bool &finalize_topo)
377 {
378  int np = points.Size()/3;
379  Dim = -1;
380  NumOfElements = cell_types.Size();
381  elements.SetSize(NumOfElements);
382 
383  int order = -1;
384  bool legacy_elem = false, lagrange_elem = false;
385 
386  for (int i = 0; i < NumOfElements; i++)
387  {
388  int j = (i > 0) ? cell_offsets[i-1] : 0;
389  int ct = cell_types[i];
390  Geometry::Type geom = VTKGeometry::GetMFEMGeometry(ct);
391  elements[i] = NewElement(geom);
392  if (cell_attributes.Size() > 0)
393  {
394  elements[i]->SetAttribute(cell_attributes[i]);
395  }
396  // VTK ordering of vertices is the same as MFEM ordering of vertices
397  // for all element types *except* prisms, which require a permutation
398  if (geom == Geometry::PRISM && ct != VTKGeometry::LAGRANGE_PRISM)
399  {
400  int prism_vertices[6];
401  for (int k=0; k<6; ++k)
402  {
403  prism_vertices[k] = cell_data[j+VTKGeometry::PrismMap[k]];
404  }
405  elements[i]->SetVertices(prism_vertices);
406  }
407  else
408  {
409  elements[i]->SetVertices(&cell_data[j]);
410  }
411 
412  int elem_dim = Geometry::Dimension[geom];
413  int elem_order = VTKGeometry::GetOrder(ct, cell_offsets[i] - j);
414 
415  if (VTKGeometry::IsLagrange(ct)) { lagrange_elem = true; }
416  else { legacy_elem = true; }
417 
418  MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
419  "Elements with different dimensions are not supported");
420  MFEM_VERIFY(order == -1 || order == elem_order,
421  "Elements with different orders are not supported");
422  MFEM_VERIFY(legacy_elem != lagrange_elem,
423  "Mixing of legacy and Lagrange cell types is not supported");
424  Dim = elem_dim;
425  order = elem_order;
426  }
427 
428  // determine spaceDim based on min/max differences detected each dimension
429  spaceDim = 0;
430  if (np > 0)
431  {
432  double min_value, max_value;
433  for (int d = 3; d > 0; --d)
434  {
435  min_value = max_value = points(3*0 + d-1);
436  for (int i = 1; i < np; i++)
437  {
438  min_value = std::min(min_value, points(3*i + d-1));
439  max_value = std::max(max_value, points(3*i + d-1));
440  if (min_value != max_value)
441  {
442  spaceDim = d;
443  break;
444  }
445  }
446  if (spaceDim > 0) { break; }
447  }
448  }
449 
450  if (order == 1 && !lagrange_elem)
451  {
452  NumOfVertices = np;
453  vertices.SetSize(np);
454  for (int i = 0; i < np; i++)
455  {
456  vertices[i](0) = points(3*i+0);
457  vertices[i](1) = points(3*i+1);
458  vertices[i](2) = points(3*i+2);
459  }
460  // No boundary is defined in a VTK mesh
461  NumOfBdrElements = 0;
462  FinalizeTopology();
463  CheckElementOrientation(true);
464  }
465  else
466  {
467  // The following section of code is shared for legacy quadratic and the
468  // Lagrange high order elements
469  curved = 1;
470 
471  // generate new enumeration for the vertices
472  Array<int> pts_dof(np);
473  pts_dof = -1;
474  // mark vertex points
475  for (int i = 0; i < NumOfElements; i++)
476  {
477  int *v = elements[i]->GetVertices();
478  int nv = elements[i]->GetNVertices();
479  for (int j = 0; j < nv; j++)
480  {
481  if (pts_dof[v[j]] == -1) { pts_dof[v[j]] = 0; }
482  }
483  }
484 
485  // The following loop reorders pts_dofs so vertices are visited in
486  // canonical order
487 
488  // Keep the original ordering of the vertices
489  NumOfVertices = 0;
490  for (int i = 0; i < np; i++)
491  {
492  if (pts_dof[i] != -1)
493  {
494  pts_dof[i] = NumOfVertices++;
495  }
496  }
497  // update the element vertices
498  for (int i = 0; i < NumOfElements; i++)
499  {
500  int *v = elements[i]->GetVertices();
501  int nv = elements[i]->GetNVertices();
502  for (int j = 0; j < nv; j++)
503  {
504  v[j] = pts_dof[v[j]];
505  }
506  }
507  // Define the 'vertices' from the 'points' through the 'pts_dof' map
508  vertices.SetSize(NumOfVertices);
509  for (int i = 0; i < np; i++)
510  {
511  int j = pts_dof[i];
512  if (j != -1)
513  {
514  vertices[j](0) = points(3*i+0);
515  vertices[j](1) = points(3*i+1);
516  vertices[j](2) = points(3*i+2);
517  }
518  }
519 
520  // No boundary is defined in a VTK mesh
521  NumOfBdrElements = 0;
522 
523  // Generate faces and edges so that we can define FE space on the mesh
524  FinalizeTopology();
525 
527  FiniteElementSpace *fes;
528  if (legacy_elem)
529  {
530  // Define quadratic FE space
531  fec = new QuadraticFECollection;
532  fes = new FiniteElementSpace(this, fec, spaceDim);
533  Nodes = new GridFunction(fes);
534  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
535  own_nodes = 1;
536 
537  // Map vtk points to edge/face/element dofs
538  Array<int> dofs;
539  for (int i = 0; i < NumOfElements; i++)
540  {
541  fes->GetElementDofs(i, dofs);
542  const int *vtk_mfem;
543  switch (elements[i]->GetGeometryType())
544  {
545  case Geometry::TRIANGLE:
546  case Geometry::SQUARE:
547  vtk_mfem = vtk_quadratic_hex; break; // identity map
548  case Geometry::TETRAHEDRON:
549  vtk_mfem = vtk_quadratic_tet; break;
550  case Geometry::CUBE:
551  vtk_mfem = vtk_quadratic_hex; break;
552  case Geometry::PRISM:
553  vtk_mfem = vtk_quadratic_wedge; break;
554  case Geometry::PYRAMID:
555  vtk_mfem = vtk_quadratic_pyramid; break;
556  default:
557  vtk_mfem = NULL; // suppress a warning
558  break;
559  }
560 
561  int offset = (i == 0) ? 0 : cell_offsets[i-1];
562  for (int j = 0; j < dofs.Size(); j++)
563  {
564  if (pts_dof[cell_data[offset+j]] == -1)
565  {
566  pts_dof[cell_data[offset+j]] = dofs[vtk_mfem[j]];
567  }
568  else
569  {
570  if (pts_dof[cell_data[offset+j]] != dofs[vtk_mfem[j]])
571  {
572  MFEM_ABORT("VTK mesh: inconsistent quadratic mesh!");
573  }
574  }
575  }
576  }
577  }
578  else
579  {
580  // Define H1 FE space
581  fec = new H1_FECollection(order,Dim,BasisType::ClosedUniform);
582  fes = new FiniteElementSpace(this, fec, spaceDim);
583  Nodes = new GridFunction(fes);
584  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
585  own_nodes = 1;
586  Array<int> dofs;
587 
588  std::map<Geometry::Type,Array<int>> vtk_inv_maps;
589  std::map<Geometry::Type,const Array<int>*> lex_orderings;
590 
591  int i, n;
592  for (n = i = 0; i < NumOfElements; i++)
593  {
594  Geometry::Type geom = GetElementBaseGeometry(i);
595  fes->GetElementDofs(i, dofs);
596 
597  Array<int> &vtk_inv_map = vtk_inv_maps[geom];
598  if (vtk_inv_map.Size() == 0)
599  {
600  Array<int> vtk_map;
601  CreateVTKElementConnectivity(vtk_map, geom, order);
602  vtk_inv_map.SetSize(vtk_map.Size());
603  for (int j=0; j<vtk_map.Size(); ++j)
604  {
605  vtk_inv_map[vtk_map[j]] = j;
606  }
607  }
608  const Array<int> *&lex_ordering = lex_orderings[geom];
609  if (!lex_ordering)
610  {
611  const FiniteElement *fe = fes->GetFE(i);
612  const NodalFiniteElement *nodal_fe =
613  dynamic_cast<const NodalFiniteElement*>(fe);
614  MFEM_ASSERT(nodal_fe != NULL, "Unsupported element type");
615  lex_ordering = &nodal_fe->GetLexicographicOrdering();
616  }
617 
618  for (int lex_idx = 0; lex_idx < dofs.Size(); lex_idx++)
619  {
620  int mfem_idx = (*lex_ordering)[lex_idx];
621  int vtk_idx = vtk_inv_map[lex_idx];
622  int pt_idx = cell_data[n + vtk_idx];
623  if (pts_dof[pt_idx] == -1)
624  {
625  pts_dof[pt_idx] = dofs[mfem_idx];
626  }
627  else
628  {
629  if (pts_dof[pt_idx] != dofs[mfem_idx])
630  {
631  MFEM_ABORT("VTK mesh: inconsistent Lagrange mesh!");
632  }
633  }
634  }
635  n += dofs.Size();
636  }
637  }
638  // Define the 'Nodes' from the 'points' through the 'pts_dof' map
639  Array<int> dofs;
640  for (int i = 0; i < np; i++)
641  {
642  dofs.SetSize(1);
643  if (pts_dof[i] != -1)
644  {
645  dofs[0] = pts_dof[i];
646  fes->DofsToVDofs(dofs);
647  for (int d = 0; d < dofs.Size(); d++)
648  {
649  (*Nodes)(dofs[d]) = points(3*i+d);
650  }
651  }
652  }
653  read_gf = 0;
654  }
655 }
656 
657 namespace vtk_xml
658 {
659 
660 using namespace tinyxml2;
661 
662 /// Return false if either string is NULL or if the strings differ, return true
663 /// if the strings are the same.
664 bool StringCompare(const char *s1, const char *s2)
665 {
666  if (s1 == NULL || s2 == NULL) { return false; }
667  return strcmp(s1, s2) == 0;
668 }
669 
670 /// Abstract base class for reading contiguous arrays of (potentially
671 /// compressed, potentially base-64 encoded) binary data from a buffer into a
672 /// destination array. The types of the source and destination arrays may be
673 /// different (e.g. read data of type uint8_t into destination array of
674 /// uint32_t), which is handled by the templated derived class @a BufferReader.
675 struct BufferReaderBase
676 {
677  enum HeaderType { UINT32_HEADER, UINT64_HEADER };
678  virtual void ReadBinary(const char *buf, void *dest, int n) const = 0;
679  virtual void ReadBase64(const char *txt, void *dest, int n) const = 0;
680  virtual ~BufferReaderBase() { }
681 };
682 
683 /// Read an array of source data stored as (potentially compressed, potentially
684 /// base-64 encoded) into a destination array. The types of the elements in the
685 /// source array are given by template parameter @a F ("from") and the types of
686 /// the elements of the destination array are given by @a T ("to"). The binary
687 /// data has a header, which is one integer if the data is uncompressed, and is
688 /// four integers if the data is compressed. The integers may either by uint32_t
689 /// or uint64_t, according to the @a header_type. If the data is compressed and
690 /// base-64 encoded, then the header is encoded separately from the data. If the
691 /// data is uncompressed and base-64 encoded, then the header and data are
692 /// encoded together.
693 template <typename T, typename F>
694 struct BufferReader : BufferReaderBase
695 {
696  bool compressed;
697  HeaderType header_type;
698  BufferReader(bool compressed_, HeaderType header_type_)
699  : compressed(compressed_), header_type(header_type_) { }
700 
701  /// Return the number of bytes of each header entry.
702  size_t HeaderEntrySize() const
703  {
704  return header_type == UINT64_HEADER ? sizeof(uint64_t) : sizeof(uint32_t);
705  }
706 
707  /// Return the value of the header entry pointer to by @a header_buf. The
708  /// value is stored as either uint32_t or uint64_t, according to the @a
709  /// header_type, and is returned as uint64_t.
710  uint64_t ReadHeaderEntry(const char *header_buf) const
711  {
712  return (header_type == UINT64_HEADER) ? bin_io::read<uint64_t>(header_buf)
713  : bin_io::read<uint32_t>(header_buf);
714  }
715 
716  /// Return the number of bytes in the header. The header consists of one
717  /// integer if the data is uncompressed, and @a N + 3 integers if the data is
718  /// compressed, where @a N is the number of blocks. The integers are either
719  /// 32 or 64 bytes depending on the value of @a header_type. The number of
720  /// blocks is determined by reading the first integer (of type @a
721  /// header_type) pointed to by @a header_buf.
722  int NumHeaderBytes(const char *header_buf) const
723  {
724  if (!compressed) { return HeaderEntrySize(); }
725  return (3 + ReadHeaderEntry(header_buf))*HeaderEntrySize();
726  }
727 
728  /// Read @a n elements of type @a F from the source buffer @a buf into the
729  /// (pre-allocated) destination array of elements of type @a T stored in
730  /// the buffer @a dest_void. The header is stored @b separately from the
731  /// rest of the data, in the buffer @a header_buf. The data buffer @a buf
732  /// does @b not contain a header.
733  void ReadBinaryWithHeader(const char *header_buf, const char *buf,
734  void *dest_void, int n) const
735  {
736  std::vector<char> uncompressed_data;
737  T *dest = static_cast<T*>(dest_void);
738 
739  if (compressed)
740  {
741 #ifdef MFEM_USE_ZLIB
742  // The header has format (where header_t is uint32_t or uint64_t):
743  // header_t number_of_blocks;
744  // header_t uncompressed_block_size;
745  // header_t uncompressed_last_block_size;
746  // header_t compressed_size[number_of_blocks];
747  int header_entry_size = HeaderEntrySize();
748  int nblocks = ReadHeaderEntry(header_buf);
749  header_buf += header_entry_size;
750  std::vector<int> header(nblocks + 2);
751  for (int i=0; i<nblocks+2; ++i)
752  {
753  header[i] = ReadHeaderEntry(header_buf);
754  header_buf += header_entry_size;
755  }
756  uncompressed_data.resize((nblocks-1)*header[0] + header[1]);
757  Bytef *dest_ptr = (Bytef *)uncompressed_data.data();
758  Bytef *dest_start = dest_ptr;
759  const Bytef *source_ptr = (const Bytef *)buf;
760  for (int i=0; i<nblocks; ++i)
761  {
762  uLongf source_len = header[i+2];
763  uLong dest_len = (i == nblocks-1) ? header[1] : header[0];
764  int res = uncompress(dest_ptr, &dest_len, source_ptr, source_len);
765  MFEM_VERIFY(res == Z_OK, "Error uncompressing");
766  dest_ptr += dest_len;
767  source_ptr += source_len;
768  }
769  MFEM_VERIFY(int(sizeof(F)*n) == (dest_ptr - dest_start),
770  "AppendedData: wrong data size");
771  buf = uncompressed_data.data();
772 #else
773  MFEM_ABORT("MFEM must be compiled with zlib enabled to uncompress.")
774 #endif
775  }
776  else
777  {
778  // Each "data block" is preceded by a header that is either UInt32 or
779  // UInt64. The rest of the data follows.
780  uint64_t data_size;
781  if (header_type == UINT32_HEADER)
782  {
783  uint32_t *data_size_32 = (uint32_t *)header_buf;
784  data_size = *data_size_32;
785  }
786  else
787  {
788  uint64_t *data_size_64 = (uint64_t *)header_buf;
789  data_size = *data_size_64;
790  }
791  MFEM_VERIFY(sizeof(F)*n == data_size, "AppendedData: wrong data size");
792  }
793 
794  if (std::is_same<T, F>::value)
795  {
796  // Special case: no type conversions necessary, so can just memcpy
797  memcpy(dest, buf, sizeof(T)*n);
798  }
799  else
800  {
801  for (int i=0; i<n; ++i)
802  {
803  // Read binary data as type F, place in array as type T
804  dest[i] = bin_io::read<F>(buf + i*sizeof(F));
805  }
806  }
807  }
808 
809  /// Read @a n elements of type @a F from source buffer @a buf into
810  /// (pre-allocated) array of elements of type @a T stored in destination
811  /// buffer @a dest. The input buffer contains both the header and the data.
812  void ReadBinary(const char *buf, void *dest, int n) const override
813  {
814  ReadBinaryWithHeader(buf, buf + NumHeaderBytes(buf), dest, n);
815  }
816 
817  /// Read @a n elements of type @a F from base-64 encoded source buffer into
818  /// (pre-allocated) array of elements of type @a T stored in destination
819  /// buffer @a dest. The base-64-encoded data is given by the null-terminated
820  /// string @a txt, which contains both the header and the data.
821  void ReadBase64(const char *txt, void *dest, int n) const override
822  {
823  // Skip whitespace
824  while (*txt)
825  {
826  if (*txt != ' ' && *txt != '\n') { break; }
827  ++txt;
828  }
829  if (compressed)
830  {
831  // Decode the first entry of the header, which we need to determine
832  // how long the rest of the header is.
833  std::vector<char> nblocks_buf;
834  int nblocks_b64 = bin_io::NumBase64Chars(HeaderEntrySize());
835  bin_io::DecodeBase64(txt, nblocks_b64, nblocks_buf);
836  std::vector<char> data, header;
837  // Compute number of characters needed to encode header in base 64,
838  // then round to nearest multiple of 4 to take padding into account.
839  int header_b64 = bin_io::NumBase64Chars(NumHeaderBytes(nblocks_buf.data()));
840  // If data is compressed, header is encoded separately
841  bin_io::DecodeBase64(txt, header_b64, header);
842  bin_io::DecodeBase64(txt + header_b64, strlen(txt)-header_b64, data);
843  ReadBinaryWithHeader(header.data(), data.data(), dest, n);
844  }
845  else
846  {
847  std::vector<char> data;
848  bin_io::DecodeBase64(txt, strlen(txt), data);
849  ReadBinary(data.data(), dest, n);
850  }
851  }
852 };
853 
854 /// Class to read data from VTK's @a DataArary elements. Each @a DataArray can
855 /// contain inline ASCII data, inline base-64-encoded data (potentially
856 /// compressed), or reference "appended data", which may be raw or base-64, and
857 /// may be compressed or uncompressed.
858 struct XMLDataReader
859 {
860  const char *appended_data, *byte_order, *compressor;
861  enum AppendedDataEncoding { RAW, BASE64 };
862  map<string,BufferReaderBase*> type_map;
863  AppendedDataEncoding encoding;
864 
865  /// Create the data reader, where @a vtk is the @a VTKFile XML element, and
866  /// @a vtu is the child @a UnstructuredGrid XML element. This will determine
867  /// the header type (32 or 64 bit integers) and whether compression is
868  /// enabled or not. The appended data will be loaded.
869  XMLDataReader(const XMLElement *vtk, const XMLElement *vtu)
870  {
871  // Determine whether binary data header is 32 or 64 bit integer
872  BufferReaderBase::HeaderType htype;
873  if (StringCompare(vtk->Attribute("header_type"), "UInt64"))
874  {
875  htype = BufferReaderBase::UINT64_HEADER;
876  }
877  else
878  {
879  htype = BufferReaderBase::UINT32_HEADER;
880  }
881 
882  // Get the byte order of the file (will check if we encounter binary data)
883  byte_order = vtk->Attribute("byte_order");
884 
885  // Get the compressor. We will check that MFEM can handle the compression
886  // if we encounter binary data.
887  compressor = vtk->Attribute("compressor");
888  bool compressed = (compressor != NULL);
889 
890  // Find the appended data.
891  appended_data = NULL;
892  for (const XMLElement *xml_elem = vtu->NextSiblingElement();
893  xml_elem != NULL;
894  xml_elem = xml_elem->NextSiblingElement())
895  {
896  if (StringCompare(xml_elem->Name(), "AppendedData"))
897  {
898  const char *encoding_str = xml_elem->Attribute("encoding");
899  if (StringCompare(encoding_str, "raw"))
900  {
901  appended_data = xml_elem->GetAppendedData();
902  encoding = RAW;
903  }
904  else if (StringCompare(encoding_str, "base64"))
905  {
906  appended_data = xml_elem->GetText();
907  encoding = BASE64;
908  }
909  MFEM_VERIFY(appended_data != NULL, "Invalid AppendedData");
910  // Appended data follows first underscore
911  bool found_leading_underscore = false;
912  while (*appended_data)
913  {
914  ++appended_data;
915  if (*appended_data == '_')
916  {
917  found_leading_underscore = true;
918  ++appended_data;
919  break;
920  }
921  }
922  MFEM_VERIFY(found_leading_underscore, "Invalid AppendedData");
923  break;
924  }
925  }
926 
927  type_map["Int8"] = new BufferReader<int, int8_t>(compressed, htype);
928  type_map["Int16"] = new BufferReader<int, int16_t>(compressed, htype);
929  type_map["Int32"] = new BufferReader<int, int32_t>(compressed, htype);
930  type_map["Int64"] = new BufferReader<int, int64_t>(compressed, htype);
931  type_map["UInt8"] = new BufferReader<int, uint8_t>(compressed, htype);
932  type_map["UInt16"] = new BufferReader<int, uint16_t>(compressed, htype);
933  type_map["UInt32"] = new BufferReader<int, uint32_t>(compressed, htype);
934  type_map["UInt64"] = new BufferReader<int, uint64_t>(compressed, htype);
935  type_map["Float32"] = new BufferReader<double, float>(compressed, htype);
936  type_map["Float64"] = new BufferReader<double, double>(compressed, htype);
937  }
938 
939  /// Read the @a DataArray XML element given by @a xml_elem into
940  /// (pre-allocated) destination array @a dest, where @a dest stores @a n
941  /// elements of type @a T.
942  template <typename T>
943  void Read(const XMLElement *xml_elem, T *dest, int n)
944  {
945  static const char *erstr = "Error reading XML DataArray";
946  MFEM_VERIFY(StringCompare(xml_elem->Name(), "DataArray"), erstr);
947  const char *format = xml_elem->Attribute("format");
948  if (StringCompare(format, "ascii"))
949  {
950  const char *txt = xml_elem->GetText();
951  MFEM_VERIFY(txt != NULL, erstr);
952  std::istringstream data_stream(txt);
953  for (int i=0; i<n; ++i) { data_stream >> dest[i]; }
954  }
955  else if (StringCompare(format, "appended"))
956  {
957  VerifyBinaryOptions();
958  int offset = xml_elem->IntAttribute("offset");
959  const char *type = xml_elem->Attribute("type");
960  MFEM_VERIFY(type != NULL, erstr);
961  BufferReaderBase *reader = type_map[type];
962  MFEM_VERIFY(reader != NULL, erstr);
963  MFEM_VERIFY(appended_data != NULL, "No AppendedData found");
964  if (encoding == RAW)
965  {
966  reader->ReadBinary(appended_data + offset, dest, n);
967  }
968  else
969  {
970  reader->ReadBase64(appended_data + offset, dest, n);
971  }
972  }
973  else if (StringCompare(format, "binary"))
974  {
975  VerifyBinaryOptions();
976  const char *txt = xml_elem->GetText();
977  MFEM_VERIFY(txt != NULL, erstr);
978  const char *type = xml_elem->Attribute("type");
979  if (type == NULL) { MFEM_ABORT(erstr); }
980  BufferReaderBase *reader = type_map[type];
981  if (reader == NULL) { MFEM_ABORT(erstr); }
982  reader->ReadBase64(txt, dest, n);
983  }
984  else
985  {
986  MFEM_ABORT("Invalid XML VTK DataArray format");
987  }
988  }
989 
990  /// Check that the byte order of the file is the same as the native byte
991  /// order that we're running with. We don't currently support converting
992  /// between byte orders. The byte order is only verified if we encounter
993  /// binary data.
994  void VerifyByteOrder() const
995  {
996  // Can't handle reading big endian from little endian or vice versa
997  if (byte_order && !StringCompare(byte_order, VTKByteOrder()))
998  {
999  MFEM_ABORT("Converting between different byte orders is unsupported.");
1000  }
1001  }
1002 
1003  /// Check that the compressor is compatible (MFEM currently only supports
1004  /// zlib compression). If MFEM is not compiled with zlib, then we cannot
1005  /// read binary data with compression.
1006  void VerifyCompressor() const
1007  {
1008  if (compressor && !StringCompare(compressor, "vtkZLibDataCompressor"))
1009  {
1010  MFEM_ABORT("Unsupported compressor. Only zlib is supported.")
1011  }
1012 #ifndef MFEM_USE_ZLIB
1013  MFEM_VERIFY(compressor == NULL, "MFEM must be compiled with zlib enabled "
1014  "to support reading compressed data.");
1015 #endif
1016  }
1017 
1018  /// Verify that the binary data is stored with compatible options (i.e.
1019  /// native byte order and compatible compression).
1020  void VerifyBinaryOptions() const
1021  {
1022  VerifyByteOrder();
1023  VerifyCompressor();
1024  }
1025 
1026  ~XMLDataReader()
1027  {
1028  for (auto &x : type_map) { delete x.second; }
1029  }
1030 };
1031 
1032 } // namespace vtk_xml
1033 
1034 void Mesh::ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
1035  bool &finalize_topo, const std::string &xml_prefix)
1036 {
1037  using namespace vtk_xml;
1038 
1039  static const char *erstr = "XML parsing error";
1040 
1041  // Create buffer beginning with xml_prefix, then read the rest of the stream
1042  std::vector<char> buf(xml_prefix.begin(), xml_prefix.end());
1043  std::istreambuf_iterator<char> eos;
1044  buf.insert(buf.end(), std::istreambuf_iterator<char>(input), eos);
1045  buf.push_back('\0'); // null-terminate buffer
1046 
1047  XMLDocument xml;
1048  xml.Parse(buf.data(), buf.size());
1049  if (xml.ErrorID() != XML_SUCCESS)
1050  {
1051  MFEM_ABORT("Error parsing XML VTK file.\n" << xml.ErrorStr());
1052  }
1053 
1054  const XMLElement *vtkfile = xml.FirstChildElement();
1055  MFEM_VERIFY(vtkfile, erstr);
1056  MFEM_VERIFY(StringCompare(vtkfile->Name(), "VTKFile"), erstr);
1057  const XMLElement *vtu = vtkfile->FirstChildElement();
1058  MFEM_VERIFY(vtu, erstr);
1059  MFEM_VERIFY(StringCompare(vtu->Name(), "UnstructuredGrid"), erstr);
1060 
1061  XMLDataReader data_reader(vtkfile, vtu);
1062 
1063  // Count the number of points and cells
1064  const XMLElement *piece = vtu->FirstChildElement();
1065  MFEM_VERIFY(StringCompare(piece->Name(), "Piece"), erstr);
1066  MFEM_VERIFY(piece->NextSiblingElement() == NULL,
1067  "XML VTK meshes with more than one Piece are not supported");
1068  int npts = piece->IntAttribute("NumberOfPoints");
1069  int ncells = piece->IntAttribute("NumberOfCells");
1070 
1071  // Read the points
1072  Vector points(3*npts);
1073  const XMLElement *pts_xml;
1074  for (pts_xml = piece->FirstChildElement();
1075  pts_xml != NULL;
1076  pts_xml = pts_xml->NextSiblingElement())
1077  {
1078  if (StringCompare(pts_xml->Name(), "Points"))
1079  {
1080  const XMLElement *pts_data = pts_xml->FirstChildElement();
1081  MFEM_VERIFY(pts_data->IntAttribute("NumberOfComponents") == 3,
1082  "XML VTK Points DataArray must have 3 components");
1083  data_reader.Read(pts_data, points.GetData(), points.Size());
1084  break;
1085  }
1086  }
1087  if (pts_xml == NULL) { MFEM_ABORT(erstr); }
1088 
1089  // Read the cells
1090  Array<int> cell_data, cell_offsets(ncells), cell_types(ncells);
1091  const XMLElement *cells_xml;
1092  for (cells_xml = piece->FirstChildElement();
1093  cells_xml != NULL;
1094  cells_xml = cells_xml->NextSiblingElement())
1095  {
1096  if (StringCompare(cells_xml->Name(), "Cells"))
1097  {
1098  const XMLElement *cell_data_xml = NULL;
1099  for (const XMLElement *data_xml = cells_xml->FirstChildElement();
1100  data_xml != NULL;
1101  data_xml = data_xml->NextSiblingElement())
1102  {
1103  const char *data_name = data_xml->Attribute("Name");
1104  if (StringCompare(data_name, "offsets"))
1105  {
1106  data_reader.Read(data_xml, cell_offsets.GetData(), ncells);
1107  }
1108  else if (StringCompare(data_name, "types"))
1109  {
1110  data_reader.Read(data_xml, cell_types.GetData(), ncells);
1111  }
1112  else if (StringCompare(data_name, "connectivity"))
1113  {
1114  // Have to read the connectivity after the offsets, because we
1115  // don't know how many points to read until we have the offsets
1116  // (size of connectivity array is equal to the last offset), so
1117  // store the XML element pointer and read this data later.
1118  cell_data_xml = data_xml;
1119  }
1120  }
1121  MFEM_VERIFY(cell_data_xml != NULL, erstr);
1122  int cell_data_size = cell_offsets.Last();
1123  cell_data.SetSize(cell_data_size);
1124  data_reader.Read(cell_data_xml, cell_data.GetData(), cell_data_size);
1125  break;
1126  }
1127  }
1128  if (cells_xml == NULL) { MFEM_ABORT(erstr); }
1129 
1130  // Read the element attributes, which are stored as CellData named "material"
1131  Array<int> cell_attributes;
1132  for (const XMLElement *cell_data_xml = piece->FirstChildElement();
1133  cell_data_xml != NULL;
1134  cell_data_xml = cell_data_xml->NextSiblingElement())
1135  {
1136  if (StringCompare(cell_data_xml->Name(), "CellData")
1137  && StringCompare(cell_data_xml->Attribute("Scalars"), "material"))
1138  {
1139  const XMLElement *data_xml = cell_data_xml->FirstChildElement();
1140  if (data_xml != NULL && StringCompare(data_xml->Name(), "DataArray"))
1141  {
1142  cell_attributes.SetSize(ncells);
1143  data_reader.Read(data_xml, cell_attributes.GetData(), ncells);
1144  }
1145  }
1146  }
1147 
1148  CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
1149  curved, read_gf, finalize_topo);
1150 }
1151 
1152 void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
1153  bool &finalize_topo)
1154 {
1155  // VTK resources:
1156  // * https://www.vtk.org/doc/nightly/html/vtkCellType_8h_source.html
1157  // * https://www.vtk.org/doc/nightly/html/classvtkCell.html
1158  // * https://lorensen.github.io/VTKExamples/site/VTKFileFormats
1159  // * https://www.kitware.com/products/books/VTKUsersGuide.pdf
1160 
1161  string buff;
1162  getline(input, buff); // comment line
1163  getline(input, buff);
1164  filter_dos(buff);
1165  if (buff != "ASCII")
1166  {
1167  MFEM_ABORT("VTK mesh is not in ASCII format!");
1168  return;
1169  }
1170  do
1171  {
1172  getline(input, buff);
1173  filter_dos(buff);
1174  if (!input.good()) { MFEM_ABORT("VTK mesh is not UNSTRUCTURED_GRID!"); }
1175  }
1176  while (buff != "DATASET UNSTRUCTURED_GRID");
1177 
1178  // Read the points, skipping optional sections such as the FIELD data from
1179  // VisIt's VTK export (or from Mesh::PrintVTK with field_data==1).
1180  do
1181  {
1182  input >> buff;
1183  if (!input.good())
1184  {
1185  MFEM_ABORT("VTK mesh does not have POINTS data!");
1186  }
1187  }
1188  while (buff != "POINTS");
1189 
1190  Vector points;
1191  int np;
1192  input >> np >> ws;
1193  getline(input, buff); // "double"
1194  points.Load(input, 3*np);
1195 
1196  //skip metadata
1197  // Looks like:
1198  // METADATA
1199  //INFORMATION 2
1200  //NAME L2_NORM_RANGE LOCATION vtkDataArray
1201  //DATA 2 0 5.19615
1202  //NAME L2_NORM_FINITE_RANGE LOCATION vtkDataArray
1203  //DATA 2 0 5.19615
1204  do
1205  {
1206  input >> buff;
1207  if (!input.good())
1208  {
1209  MFEM_ABORT("VTK mesh does not have CELLS data!");
1210  }
1211  }
1212  while (buff != "CELLS");
1213 
1214  // Read the cells
1215  Array<int> cell_data, cell_offsets;
1216  if (buff == "CELLS")
1217  {
1218  int ncells, n;
1219  input >> ncells >> n >> ws;
1220  cell_offsets.SetSize(ncells);
1221  cell_data.SetSize(n - ncells);
1222  int offset = 0;
1223  for (int i=0; i<ncells; ++i)
1224  {
1225  int nv;
1226  input >> nv;
1227  cell_offsets[i] = offset + nv;
1228  for (int j=0; j<nv; ++j)
1229  {
1230  input >> cell_data[offset + j];
1231  }
1232  offset += nv;
1233  }
1234  }
1235 
1236  // Read the cell types
1237  input >> ws >> buff;
1238  Array<int> cell_types;
1239  int ncells;
1240  MFEM_VERIFY(buff == "CELL_TYPES", "CELL_TYPES not provided in VTK mesh.")
1241  input >> ncells;
1242  cell_types.Load(ncells, input);
1243 
1244  while ((input.good()) && (buff != "CELL_DATA"))
1245  {
1246  input >> buff;
1247  }
1248  getline(input, buff); // finish the line
1249 
1250  // Read the cell materials
1251  // bool found_material = false;
1252  Array<int> cell_attributes;
1253  while ((input.good()))
1254  {
1255  getline(input, buff);
1256  if (buff.rfind("POINT_DATA") == 0)
1257  {
1258  break; // We have entered the POINT_DATA block. Quit.
1259  }
1260  else if (buff.rfind("SCALARS material") == 0)
1261  {
1262  getline(input, buff); // LOOKUP_TABLE default
1263  if (buff.rfind("LOOKUP_TABLE default") != 0)
1264  {
1265  MFEM_ABORT("Invalid LOOKUP_TABLE for material array in VTK file.");
1266  }
1267  cell_attributes.Load(ncells, input);
1268  // found_material = true;
1269  break;
1270  }
1271  }
1272 
1273  // if (!found_material)
1274  // {
1275  // MFEM_WARNING("Material array not found in VTK file. "
1276  // "Assuming uniform material composition.");
1277  // }
1278 
1279  CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
1280  curved, read_gf, finalize_topo);
1281 } // end ReadVTKMesh
1282 
1283 void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
1284 {
1285  NURBSext = new NURBSExtension(input);
1286 
1287  Dim = NURBSext->Dimension();
1288  NumOfVertices = NURBSext->GetNV();
1289  NumOfElements = NURBSext->GetNE();
1290  NumOfBdrElements = NURBSext->GetNBE();
1291 
1292  NURBSext->GetElementTopo(elements);
1293  NURBSext->GetBdrElementTopo(boundary);
1294 
1295  vertices.SetSize(NumOfVertices);
1296  curved = 1;
1297  if (NURBSext->HavePatches())
1298  {
1299  NURBSFECollection *fec = new NURBSFECollection(NURBSext->GetOrder());
1300  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
1301  Ordering::byVDIM);
1302  Nodes = new GridFunction(fes);
1303  Nodes->MakeOwner(fec);
1304  NURBSext->SetCoordsFromPatches(*Nodes);
1305  own_nodes = 1;
1306  read_gf = 0;
1307  spaceDim = Nodes->VectorDim();
1308  for (int i = 0; i < spaceDim; i++)
1309  {
1310  Vector vert_val;
1311  Nodes->GetNodalValues(vert_val, i+1);
1312  for (int j = 0; j < NumOfVertices; j++)
1313  {
1314  vertices[j](i) = vert_val(j);
1315  }
1316  }
1317  }
1318  else
1319  {
1320  read_gf = 1;
1321  }
1322 }
1323 
1324 void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
1325 {
1326  // Initialize to negative numbers so that we know if they've been set. We're
1327  // using Element::POINT as our flag, since we're not going to make a 0D mesh,
1328  // ever.
1329  int nx = -1;
1330  int ny = -1;
1331  int nz = -1;
1332  double sx = -1.0;
1333  double sy = -1.0;
1334  double sz = -1.0;
1335  Element::Type type = Element::POINT;
1336 
1337  while (true)
1338  {
1339  skip_comment_lines(input, '#');
1340  // Break out if we reached the end of the file after gobbling up the
1341  // whitespace and comments after the last keyword.
1342  if (!input.good())
1343  {
1344  break;
1345  }
1346 
1347  // Read the next keyword
1348  std::string name;
1349  input >> name;
1350  input >> std::ws;
1351  // Make sure there's an equal sign
1352  MFEM_VERIFY(input.get() == '=',
1353  "Inline mesh expected '=' after keyword " << name);
1354  input >> std::ws;
1355 
1356  if (name == "nx")
1357  {
1358  input >> nx;
1359  }
1360  else if (name == "ny")
1361  {
1362  input >> ny;
1363  }
1364  else if (name == "nz")
1365  {
1366  input >> nz;
1367  }
1368  else if (name == "sx")
1369  {
1370  input >> sx;
1371  }
1372  else if (name == "sy")
1373  {
1374  input >> sy;
1375  }
1376  else if (name == "sz")
1377  {
1378  input >> sz;
1379  }
1380  else if (name == "type")
1381  {
1382  std::string eltype;
1383  input >> eltype;
1384  if (eltype == "segment")
1385  {
1386  type = Element::SEGMENT;
1387  }
1388  else if (eltype == "quad")
1389  {
1390  type = Element::QUADRILATERAL;
1391  }
1392  else if (eltype == "tri")
1393  {
1394  type = Element::TRIANGLE;
1395  }
1396  else if (eltype == "hex")
1397  {
1398  type = Element::HEXAHEDRON;
1399  }
1400  else if (eltype == "wedge")
1401  {
1402  type = Element::WEDGE;
1403  }
1404  else if (eltype == "pyramid")
1405  {
1406  type = Element::PYRAMID;
1407  }
1408  else if (eltype == "tet")
1409  {
1410  type = Element::TETRAHEDRON;
1411  }
1412  else
1413  {
1414  MFEM_ABORT("unrecognized element type (read '" << eltype
1415  << "') in inline mesh format. "
1416  "Allowed: segment, tri, quad, tet, hex, wedge");
1417  }
1418  }
1419  else
1420  {
1421  MFEM_ABORT("unrecognized keyword (" << name
1422  << ") in inline mesh format. "
1423  "Allowed: nx, ny, nz, type, sx, sy, sz");
1424  }
1425 
1426  input >> std::ws;
1427  // Allow an optional semi-colon at the end of each line.
1428  if (input.peek() == ';')
1429  {
1430  input.get();
1431  }
1432 
1433  // Done reading file
1434  if (!input)
1435  {
1436  break;
1437  }
1438  }
1439 
1440  // Now make the mesh.
1441  if (type == Element::SEGMENT)
1442  {
1443  MFEM_VERIFY(nx > 0 && sx > 0.0,
1444  "invalid 1D inline mesh format, all values must be "
1445  "positive\n"
1446  << " nx = " << nx << "\n"
1447  << " sx = " << sx << "\n");
1448  Make1D(nx, sx);
1449  }
1450  else if (type == Element::TRIANGLE || type == Element::QUADRILATERAL)
1451  {
1452  MFEM_VERIFY(nx > 0 && ny > 0 && sx > 0.0 && sy > 0.0,
1453  "invalid 2D inline mesh format, all values must be "
1454  "positive\n"
1455  << " nx = " << nx << "\n"
1456  << " ny = " << ny << "\n"
1457  << " sx = " << sx << "\n"
1458  << " sy = " << sy << "\n");
1459  Make2D(nx, ny, type, sx, sy, generate_edges, true);
1460  }
1461  else if (type == Element::TETRAHEDRON || type == Element::WEDGE ||
1462  type == Element::HEXAHEDRON || type == Element::PYRAMID)
1463  {
1464  MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
1465  sx > 0.0 && sy > 0.0 && sz > 0.0,
1466  "invalid 3D inline mesh format, all values must be "
1467  "positive\n"
1468  << " nx = " << nx << "\n"
1469  << " ny = " << ny << "\n"
1470  << " nz = " << nz << "\n"
1471  << " sx = " << sx << "\n"
1472  << " sy = " << sy << "\n"
1473  << " sz = " << sz << "\n");
1474  Make3D(nx, ny, nz, type, sx, sy, sz, true);
1475  // TODO: maybe have an option in the file to control ordering?
1476  }
1477  else
1478  {
1479  MFEM_ABORT("For inline mesh, must specify an element type ="
1480  " [segment, tri, quad, tet, hex, wedge]");
1481  }
1482 }
1483 
1484 void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
1485 {
1486  string buff;
1487  double version;
1488  int binary, dsize;
1489  input >> version >> binary >> dsize;
1490  if (version < 2.2)
1491  {
1492  MFEM_ABORT("Gmsh file version < 2.2");
1493  }
1494  if (dsize != sizeof(double))
1495  {
1496  MFEM_ABORT("Gmsh file : dsize != sizeof(double)");
1497  }
1498  getline(input, buff);
1499  // There is a number 1 in binary format
1500  if (binary)
1501  {
1502  int one;
1503  input.read(reinterpret_cast<char*>(&one), sizeof(one));
1504  if (one != 1)
1505  {
1506  MFEM_ABORT("Gmsh file : wrong binary format");
1507  }
1508  }
1509 
1510  // A map between a serial number of the vertex and its number in the file
1511  // (there may be gaps in the numbering, and also Gmsh enumerates vertices
1512  // starting from 1, not 0)
1513  map<int, int> vertices_map;
1514 
1515  // Gmsh always outputs coordinates in 3D, but MFEM distinguishes between the
1516  // mesh element dimension (Dim) and the dimension of the space in which the
1517  // mesh is embedded (spaceDim). For example, a 2D MFEM mesh has Dim = 2 and
1518  // spaceDim = 2, while a 2D surface mesh in 3D has Dim = 2 but spaceDim = 3.
1519  // Below we set spaceDim by measuring the mesh bounding box and checking for
1520  // a lower dimensional subspace. The assumption is that the mesh is at least
1521  // 2D if the y-dimension of the box is non-trivial and 3D if the z-dimension
1522  // is non-trivial. Note that with these assumptions a 2D mesh parallel to the
1523  // yz plane will be considered a surface mesh embedded in 3D whereas the same
1524  // 2D mesh parallel to the xy plane will be considered a 2D mesh.
1525  double bb_tol = 1e-14;
1526  double bb_min[3];
1527  double bb_max[3];
1528 
1529  // Mesh order
1530  int mesh_order = 1;
1531 
1532  // Mesh type
1533  bool periodic = false;
1534 
1535  // Vector field to store uniformly spaced Gmsh high order mesh coords
1536  GridFunction Nodes_gf;
1537 
1538  // Read the lines of the mesh file. If we face specific keyword, we'll treat
1539  // the section.
1540  while (input >> buff)
1541  {
1542  if (buff == "$Nodes") // reading mesh vertices
1543  {
1544  input >> NumOfVertices;
1545  getline(input, buff);
1546  vertices.SetSize(NumOfVertices);
1547  int serial_number;
1548  const int gmsh_dim = 3; // Gmsh always outputs 3 coordinates
1549  double coord[gmsh_dim];
1550  for (int ver = 0; ver < NumOfVertices; ++ver)
1551  {
1552  if (binary)
1553  {
1554  input.read(reinterpret_cast<char*>(&serial_number), sizeof(int));
1555  input.read(reinterpret_cast<char*>(coord), gmsh_dim*sizeof(double));
1556  }
1557  else // ASCII
1558  {
1559  input >> serial_number;
1560  for (int ci = 0; ci < gmsh_dim; ++ci)
1561  {
1562  input >> coord[ci];
1563  }
1564  }
1565  vertices[ver] = Vertex(coord, gmsh_dim);
1566  vertices_map[serial_number] = ver;
1567 
1568  for (int ci = 0; ci < gmsh_dim; ++ci)
1569  {
1570  bb_min[ci] = (ver == 0) ? coord[ci] :
1571  std::min(bb_min[ci], coord[ci]);
1572  bb_max[ci] = (ver == 0) ? coord[ci] :
1573  std::max(bb_max[ci], coord[ci]);
1574  }
1575  }
1576  double bb_size = std::max(bb_max[0] - bb_min[0],
1577  std::max(bb_max[1] - bb_min[1],
1578  bb_max[2] - bb_min[2]));
1579  spaceDim = 1;
1580  if (bb_max[1] - bb_min[1] > bb_size * bb_tol)
1581  {
1582  spaceDim++;
1583  }
1584  if (bb_max[2] - bb_min[2] > bb_size * bb_tol)
1585  {
1586  spaceDim++;
1587  }
1588 
1589  if (static_cast<int>(vertices_map.size()) != NumOfVertices)
1590  {
1591  MFEM_ABORT("Gmsh file : vertices indices are not unique");
1592  }
1593  } // section '$Nodes'
1594  else if (buff == "$Elements") // reading mesh elements
1595  {
1596  int num_of_all_elements;
1597  input >> num_of_all_elements;
1598  // = NumOfElements + NumOfBdrElements + (maybe, PhysicalPoints)
1599  getline(input, buff);
1600 
1601  int serial_number; // serial number of an element
1602  int type_of_element; // ID describing a type of a mesh element
1603  int n_tags; // number of different tags describing an element
1604  int phys_domain; // element's attribute
1605  int elem_domain; // another element's attribute (rarely used)
1606  int n_partitions; // number of partitions where an element takes place
1607 
1608  // number of nodes for each type of Gmsh elements, type is the index of
1609  // the array + 1
1610  int nodes_of_gmsh_element[] =
1611  {
1612  2, // 2-node line.
1613  3, // 3-node triangle.
1614  4, // 4-node quadrangle.
1615  4, // 4-node tetrahedron.
1616  8, // 8-node hexahedron.
1617  6, // 6-node prism.
1618  5, // 5-node pyramid.
1619  3, /* 3-node second order line (2 nodes associated with the
1620  vertices and 1 with the edge). */
1621  6, /* 6-node second order triangle (3 nodes associated with the
1622  vertices and 3 with the edges). */
1623  9, /* 9-node second order quadrangle (4 nodes associated with the
1624  vertices, 4 with the edges and 1 with the face). */
1625  10, /* 10-node second order tetrahedron (4 nodes associated with
1626  the vertices and 6 with the edges). */
1627  27, /* 27-node second order hexahedron (8 nodes associated with the
1628  vertices, 12 with the edges, 6 with the faces and 1 with the
1629  volume). */
1630  18, /* 18-node second order prism (6 nodes associated with the
1631  vertices, 9 with the edges and 3 with the quadrangular
1632  faces). */
1633  14, /* 14-node second order pyramid (5 nodes associated with the
1634  vertices, 8 with the edges and 1 with the quadrangular
1635  face). */
1636  1, // 1-node point.
1637  8, /* 8-node second order quadrangle (4 nodes associated with the
1638  vertices and 4 with the edges). */
1639  20, /* 20-node second order hexahedron (8 nodes associated with the
1640  vertices and 12 with the edges). */
1641  15, /* 15-node second order prism (6 nodes associated with the
1642  vertices and 9 with the edges). */
1643  13, /* 13-node second order pyramid (5 nodes associated with the
1644  vertices and 8 with the edges). */
1645  9, /* 9-node third order incomplete triangle (3 nodes associated
1646  with the vertices, 6 with the edges) */
1647  10, /* 10-node third order triangle (3 nodes associated with the
1648  vertices, 6 with the edges, 1 with the face) */
1649  12, /* 12-node fourth order incomplete triangle (3 nodes associated
1650  with the vertices, 9 with the edges) */
1651  15, /* 15-node fourth order triangle (3 nodes associated with the
1652  vertices, 9 with the edges, 3 with the face) */
1653  15, /* 15-node fifth order incomplete triangle (3 nodes associated
1654  with the vertices, 12 with the edges) */
1655  21, /* 21-node fifth order complete triangle (3 nodes associated
1656  with the vertices, 12 with the edges, 6 with the face) */
1657  4, /* 4-node third order edge (2 nodes associated with the
1658  vertices, 2 internal to the edge) */
1659  5, /* 5-node fourth order edge (2 nodes associated with the
1660  vertices, 3 internal to the edge) */
1661  6, /* 6-node fifth order edge (2 nodes associated with the
1662  vertices, 4 internal to the edge) */
1663  20, /* 20-node third order tetrahedron (4 nodes associated with the
1664  vertices, 12 with the edges, 4 with the faces) */
1665  35, /* 35-node fourth order tetrahedron (4 nodes associated with
1666  the vertices, 18 with the edges, 12 with the faces, and 1
1667  with the volume) */
1668  56, /* 56-node fifth order tetrahedron (4 nodes associated with the
1669  vertices, 24 with the edges, 24 with the faces, and 4 with
1670  the volume) */
1671  -1,-1, /* unsupported tetrahedral types */
1672  -1,-1, /* unsupported polygonal and polyhedral types */
1673  16, /* 16-node third order quadrilateral (4 nodes associated with
1674  the vertices, 8 with the edges, 4 with the face) */
1675  25, /* 25-node fourth order quadrilateral (4 nodes associated with
1676  the vertices, 12 with the edges, 9 with the face) */
1677  36, /* 36-node fifth order quadrilateral (4 nodes associated with
1678  the vertices, 16 with the edges, 16 with the face) */
1679  -1,-1,-1, /* unsupported quadrilateral types */
1680  28, /* 28-node sixth order complete triangle (3 nodes associated
1681  with the vertices, 15 with the edges, 10 with the face) */
1682  36, /* 36-node seventh order complete triangle (3 nodes associated
1683  with the vertices, 18 with the edges, 15 with the face) */
1684  45, /* 45-node eighth order complete triangle (3 nodes associated
1685  with the vertices, 21 with the edges, 21 with the face) */
1686  55, /* 55-node ninth order complete triangle (3 nodes associated
1687  with the vertices, 24 with the edges, 28 with the face) */
1688  66, /* 66-node tenth order complete triangle (3 nodes associated
1689  with the vertices, 27 with the edges, 36 with the face) */
1690  49, /* 49-node sixth order quadrilateral (4 nodes associated with
1691  the vertices, 20 with the edges, 25 with the face) */
1692  64, /* 64-node seventh order quadrilateral (4 nodes associated with
1693  the vertices, 24 with the edges, 36 with the face) */
1694  81, /* 81-node eighth order quadrilateral (4 nodes associated with
1695  the vertices, 28 with the edges, 49 with the face) */
1696  100, /* 100-node ninth order quadrilateral (4 nodes associated with
1697  the vertices, 32 with the edges, 64 with the face) */
1698  121, /* 121-node tenth order quadrilateral (4 nodes associated with
1699  the vertices, 36 with the edges, 81 with the face) */
1700  -1,-1,-1,-1,-1, /* unsupported triangular types */
1701  -1,-1,-1,-1,-1, /* unsupported quadrilateral types */
1702  7, /* 7-node sixth order edge (2 nodes associated with the
1703  vertices, 5 internal to the edge) */
1704  8, /* 8-node seventh order edge (2 nodes associated with the
1705  vertices, 6 internal to the edge) */
1706  9, /* 9-node eighth order edge (2 nodes associated with the
1707  vertices, 7 internal to the edge) */
1708  10, /* 10-node ninth order edge (2 nodes associated with the
1709  vertices, 8 internal to the edge) */
1710  11, /* 11-node tenth order edge (2 nodes associated with the
1711  vertices, 9 internal to the edge) */
1712  -1, /* unsupported linear types */
1713  -1,-1,-1, /* unsupported types */
1714  84, /* 84-node sixth order tetrahedron (4 nodes associated with the
1715  vertices, 30 with the edges, 40 with the faces, and 10 with
1716  the volume) */
1717  120, /* 120-node seventh order tetrahedron (4 nodes associated with
1718  the vertices, 36 with the edges, 60 with the faces, and 20
1719  with the volume) */
1720  165, /* 165-node eighth order tetrahedron (4 nodes associated with
1721  the vertices, 42 with the edges, 84 with the faces, and 35
1722  with the volume) */
1723  220, /* 220-node ninth order tetrahedron (4 nodes associated with
1724  the vertices, 48 with the edges, 112 with the faces, and 56
1725  with the volume) */
1726  286, /* 286-node tenth order tetrahedron (4 nodes associated with
1727  the vertices, 54 with the edges, 144 with the faces, and 84
1728  with the volume) */
1729  -1,-1,-1, /* undefined types */
1730  -1,-1,-1,-1,-1, /* unsupported tetrahedral types */
1731  -1,-1,-1,-1,-1,-1, /* unsupported types */
1732  40, /* 40-node third order prism (6 nodes associated with the
1733  vertices, 18 with the edges, 14 with the faces, and 2 with
1734  the volume) */
1735  75, /* 75-node fourth order prism (6 nodes associated with the
1736  vertices, 27 with the edges, 33 with the faces, and 9 with
1737  the volume) */
1738  64, /* 64-node third order hexahedron (8 nodes associated with the
1739  vertices, 24 with the edges, 24 with the faces and 8 with
1740  the volume).*/
1741  125, /* 125-node fourth order hexahedron (8 nodes associated with
1742  the vertices, 36 with the edges, 54 with the faces and 27
1743  with the volume).*/
1744  216, /* 216-node fifth order hexahedron (8 nodes associated with the
1745  vertices, 48 with the edges, 96 with the faces and 64 with
1746  the volume).*/
1747  343, /* 343-node sixth order hexahedron (8 nodes associated with the
1748  vertices, 60 with the edges, 150 with the faces and 125 with
1749  the volume).*/
1750  512, /* 512-node seventh order hexahedron (8 nodes associated with
1751  the vertices, 72 with the edges, 216 with the faces and 216
1752  with the volume).*/
1753  729, /* 729-node eighth order hexahedron (8 nodes associated with
1754  the vertices, 84 with the edges, 294 with the faces and 343
1755  with the volume).*/
1756  1000,/* 1000-node ninth order hexahedron (8 nodes associated with
1757  the vertices, 96 with the edges, 384 with the faces and 512
1758  with the volume).*/
1759  -1,-1,-1,-1,-1,-1,-1, /* unsupported hexahedron types */
1760  126, /* 126-node fifth order prism (6 nodes associated with the
1761  vertices, 36 with the edges, 60 with the faces, and 24 with
1762  the volume) */
1763  196, /* 196-node sixth order prism (6 nodes associated with the
1764  vertices, 45 with the edges, 95 with the faces, and 50 with
1765  the volume) */
1766  288, /* 288-node seventh order prism (6 nodes associated with the
1767  vertices, 54 with the edges, 138 with the faces, and 90 with
1768  the volume) */
1769  405, /* 405-node eighth order prism (6 nodes associated with the
1770  vertices, 63 with the edges, 189 with the faces, and 147
1771  with the volume) */
1772  550, /* 550-node ninth order prism (6 nodes associated with the
1773  vertices, 72 with the edges, 248 with the faces, and 224
1774  with the volume) */
1775  -1,-1,-1,-1,-1,-1,-1, /* unsupported prism types */
1776  30, /* 30-node third order pyramid (5 nodes associated with the
1777  vertices, 16 with the edges and 8 with the faces, and 1 with
1778  the volume). */
1779  55, /* 55-node fourth order pyramid (5 nodes associated with the
1780  vertices, 24 with the edges and 21 with the faces, and 5
1781  with the volume). */
1782  91, /* 91-node fifth order pyramid (5 nodes associated with the
1783  vertices, 32 with the edges and 40 with the faces, and 14
1784  with the volume). */
1785  140, /* 140-node sixth order pyramid (5 nodes associated with the
1786  vertices, 40 with the edges and 65 with the faces, and 30
1787  with the volume). */
1788  204, /* 204-node seventh order pyramid (5 nodes associated with the
1789  vertices, 48 with the edges and 96 with the faces, and 55
1790  with the volume). */
1791  285, /* 285-node eighth order pyramid (5 nodes associated with the
1792  vertices, 56 with the edges and 133 with the faces, and 91
1793  with the volume). */
1794  385 /* 385-node ninth order pyramid (5 nodes associated with the
1795  vertices, 64 with the edges and 176 with the faces, and 140
1796  with the volume). */
1797  };
1798 
1799  /** The following mappings convert the Gmsh node orderings for high
1800  order elements to MFEM's L2 degree of freedom ordering. To support
1801  more options examine Gmsh's ordering and read off the indices in
1802  MFEM's order. For example 2nd order Gmsh quadrilaterals use the
1803  following ordering:
1804 
1805  3--6--2
1806  | | |
1807  7 8 5
1808  | | |
1809  0--4--1
1810 
1811  (from https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering)
1812 
1813  Whereas MFEM uses a tensor product ordering with the horizontal
1814  axis cycling fastest so we would read off:
1815 
1816  0 4 1 7 8 5 3 6 2
1817 
1818  This corresponds to the quad9 mapping below.
1819  */
1820  int lin3[] = {0,2,1}; // 2nd order segment
1821  int lin4[] = {0,2,3,1}; // 3rd order segment
1822  int tri6[] = {0,3,1,5,4,2}; // 2nd order triangle
1823  int tri10[] = {0,3,4,1,8,9,5,7,6,2}; // 3rd order triangle
1824  int quad9[] = {0,4,1,7,8,5,3,6,2}; // 2nd order quadrilateral
1825  int quad16[] = {0,4,5,1,11,12,13,6, // 3rd order quadrilateral
1826  10,15,14,7,3,9,8,2
1827  };
1828  int tet10[] {0,4,1,6,5,2,7,9,8,3}; // 2nd order tetrahedron
1829  int tet20[] = {0,4,5,1,9,16,6,8,7,2, // 3rd order tetrahedron
1830  11,17,15,18,19,13,10,14,12,3
1831  };
1832  int hex27[] {0,8,1,9,20,11,3,13,2, // 2nd order hexahedron
1833  10,21,12,22,26,23,15,24,14,
1834  4,16,5,17,25,18,7,19,6
1835  };
1836  int hex64[] {0,8,9,1,10,32,35,14, // 3rd order hexahedron
1837  11,33,34,15,3,19,18,2,
1838  12,36,37,16,40,56,57,44,
1839  43,59,58,45,22,49,48,20,
1840  13,39,38,17,41,60,61,47,
1841  42,63,62,46,23,50,51,21,
1842  4,24,25,5,26,52,53,28,
1843  27,55,54,29,7,31,30,6
1844  };
1845 
1846  int wdg18[] = {0,6,1,7,9,2,8,15,10, // 2nd order wedge/prism
1847  16,17,11,3,12,4,13,14,5
1848  };
1849  int wdg40[] = {0,6,7,1,8,24,12,9,13,2, // 3rd order wedge/prism
1850  10,26,27,14,30,38,34,33,35,16,
1851  11,29,28,15,31,39,37,32,36,17,
1852  3,18,19,4,20,25,22,21,23,5
1853  };
1854 
1855  int pyr14[] = {0,5,1,6,13,8,3, // 2nd order pyramid
1856  10,2,7,9,12,11,4
1857  };
1858  int pyr30[] = {0,5,6,1,7,25,28,11,8,26, // 3rd order pyramid
1859  27,12,3,16,15,2,9,21,13,22,
1860  29,23,19,24,17,10,14,20,18,4
1861  };
1862 
1863  vector<Element*> elements_0D, elements_1D, elements_2D, elements_3D;
1864  elements_0D.reserve(num_of_all_elements);
1865  elements_1D.reserve(num_of_all_elements);
1866  elements_2D.reserve(num_of_all_elements);
1867  elements_3D.reserve(num_of_all_elements);
1868 
1869  // Temporary storage for high order vertices, if present
1870  vector<Array<int>*> ho_verts_1D, ho_verts_2D, ho_verts_3D;
1871  ho_verts_1D.reserve(num_of_all_elements);
1872  ho_verts_2D.reserve(num_of_all_elements);
1873  ho_verts_3D.reserve(num_of_all_elements);
1874 
1875  // Temporary storage for order of elements
1876  vector<int> ho_el_order_1D, ho_el_order_2D, ho_el_order_3D;
1877  ho_el_order_1D.reserve(num_of_all_elements);
1878  ho_el_order_2D.reserve(num_of_all_elements);
1879  ho_el_order_3D.reserve(num_of_all_elements);
1880 
1881  // Vertex order mappings
1882  Array<int*> ho_lin(11); ho_lin = NULL;
1883  Array<int*> ho_tri(11); ho_tri = NULL;
1884  Array<int*> ho_sqr(11); ho_sqr = NULL;
1885  Array<int*> ho_tet(11); ho_tet = NULL;
1886  Array<int*> ho_hex(10); ho_hex = NULL;
1887  Array<int*> ho_wdg(10); ho_wdg = NULL;
1888  Array<int*> ho_pyr(10); ho_pyr = NULL;
1889 
1890  // Use predefined arrays at lowest orders (for efficiency)
1891  ho_lin[2] = lin3; ho_lin[3] = lin4;
1892  ho_tri[2] = tri6; ho_tri[3] = tri10;
1893  ho_sqr[2] = quad9; ho_sqr[3] = quad16;
1894  ho_tet[2] = tet10; ho_tet[3] = tet20;
1895  ho_hex[2] = hex27; ho_hex[3] = hex64;
1896  ho_wdg[2] = wdg18; ho_wdg[3] = wdg40;
1897  ho_pyr[2] = pyr14; ho_pyr[3] = pyr30;
1898 
1899  bool has_nonpositive_phys_domain = false;
1900  bool has_positive_phys_domain = false;
1901 
1902  if (binary)
1903  {
1904  int n_elem_part = 0; // partial sum of elements that are read
1905  const int header_size = 3;
1906  // header consists of 3 numbers: type of the element, number of
1907  // elements of this type, and number of tags
1908  int header[header_size];
1909  int n_elem_one_type; // number of elements of a specific type
1910 
1911  while (n_elem_part < num_of_all_elements)
1912  {
1913  input.read(reinterpret_cast<char*>(header),
1914  header_size*sizeof(int));
1915  type_of_element = header[0];
1916  n_elem_one_type = header[1];
1917  n_tags = header[2];
1918 
1919  n_elem_part += n_elem_one_type;
1920 
1921  const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
1922  vector<int> data(1+n_tags+n_elem_nodes);
1923  for (int el = 0; el < n_elem_one_type; ++el)
1924  {
1925  input.read(reinterpret_cast<char*>(&data[0]),
1926  data.size()*sizeof(int));
1927  int dd = 0; // index for data array
1928  serial_number = data[dd++];
1929  // physical domain - the most important value (to distinguish
1930  // materials with different properties)
1931  phys_domain = (n_tags > 0) ? data[dd++] : 1;
1932  // elementary domain - to distinguish different geometrical
1933  // domains (typically, it's used rarely)
1934  elem_domain = (n_tags > 1) ? data[dd++] : 0;
1935  // the number of tags is bigger than 2 if there are some
1936  // partitions (domain decompositions)
1937  n_partitions = (n_tags > 2) ? data[dd++] : 0;
1938  // we currently just skip the partitions if they exist, and go
1939  // directly to vertices describing the mesh element
1940  vector<int> vert_indices(n_elem_nodes);
1941  for (int vi = 0; vi < n_elem_nodes; ++vi)
1942  {
1943  map<int, int>::const_iterator it =
1944  vertices_map.find(data[1+n_tags+vi]);
1945  if (it == vertices_map.end())
1946  {
1947  MFEM_ABORT("Gmsh file : vertex index doesn't exist");
1948  }
1949  vert_indices[vi] = it->second;
1950  }
1951 
1952  // Non-positive attributes are not allowed in MFEM. However,
1953  // by default, Gmsh sets the physical domain of all elements
1954  // to zero. In the case that all elements have physical domain
1955  // zero, we will given them attribute 1. If only some elements
1956  // have physical domain zero, we will throw an error.
1957  if (phys_domain <= 0)
1958  {
1959  has_nonpositive_phys_domain = true;
1960  phys_domain = 1;
1961  }
1962  else
1963  {
1964  has_positive_phys_domain = true;
1965  }
1966 
1967  // initialize the mesh element
1968  int el_order = 11;
1969  switch (type_of_element)
1970  {
1971  case 1: // 2-node line
1972  case 8: // 3-node line (2nd order)
1973  case 26: // 4-node line (3rd order)
1974  case 27: // 5-node line (4th order)
1975  case 28: // 6-node line (5th order)
1976  case 62: // 7-node line (6th order)
1977  case 63: // 8-node line (7th order)
1978  case 64: // 9-node line (8th order)
1979  case 65: // 10-node line (9th order)
1980  case 66: // 11-node line (10th order)
1981  {
1982  elements_1D.push_back(
1983  new Segment(&vert_indices[0], phys_domain));
1984  if (type_of_element != 1)
1985  {
1986  el_order = n_elem_nodes - 1;
1987  Array<int> * hov = new Array<int>;
1988  hov->Append(&vert_indices[0], n_elem_nodes);
1989  ho_verts_1D.push_back(hov);
1990  ho_el_order_1D.push_back(el_order);
1991  }
1992  break;
1993  }
1994  case 2: el_order--; // 3-node triangle
1995  case 9: el_order--; // 6-node triangle (2nd order)
1996  case 21: el_order--; // 10-node triangle (3rd order)
1997  case 23: el_order--; // 15-node triangle (4th order)
1998  case 25: el_order--; // 21-node triangle (5th order)
1999  case 42: el_order--; // 28-node triangle (6th order)
2000  case 43: el_order--; // 36-node triangle (7th order)
2001  case 44: el_order--; // 45-node triangle (8th order)
2002  case 45: el_order--; // 55-node triangle (9th order)
2003  case 46: el_order--; // 66-node triangle (10th order)
2004  {
2005  elements_2D.push_back(
2006  new Triangle(&vert_indices[0], phys_domain));
2007  if (el_order > 1)
2008  {
2009  Array<int> * hov = new Array<int>;
2010  hov->Append(&vert_indices[0], n_elem_nodes);
2011  ho_verts_2D.push_back(hov);
2012  ho_el_order_2D.push_back(el_order);
2013  }
2014  break;
2015  }
2016  case 3: el_order--; // 4-node quadrangle
2017  case 10: el_order--; // 9-node quadrangle (2nd order)
2018  case 36: el_order--; // 16-node quadrangle (3rd order)
2019  case 37: el_order--; // 25-node quadrangle (4th order)
2020  case 38: el_order--; // 36-node quadrangle (5th order)
2021  case 47: el_order--; // 49-node quadrangle (6th order)
2022  case 48: el_order--; // 64-node quadrangle (7th order)
2023  case 49: el_order--; // 81-node quadrangle (8th order)
2024  case 50: el_order--; // 100-node quadrangle (9th order)
2025  case 51: el_order--; // 121-node quadrangle (10th order)
2026  {
2027  elements_2D.push_back(
2028  new Quadrilateral(&vert_indices[0], phys_domain));
2029  if (el_order > 1)
2030  {
2031  Array<int> * hov = new Array<int>;
2032  hov->Append(&vert_indices[0], n_elem_nodes);
2033  ho_verts_2D.push_back(hov);
2034  ho_el_order_2D.push_back(el_order);
2035  }
2036  break;
2037  }
2038  case 4: el_order--; // 4-node tetrahedron
2039  case 11: el_order--; // 10-node tetrahedron (2nd order)
2040  case 29: el_order--; // 20-node tetrahedron (3rd order)
2041  case 30: el_order--; // 35-node tetrahedron (4th order)
2042  case 31: el_order--; // 56-node tetrahedron (5th order)
2043  case 71: el_order--; // 84-node tetrahedron (6th order)
2044  case 72: el_order--; // 120-node tetrahedron (7th order)
2045  case 73: el_order--; // 165-node tetrahedron (8th order)
2046  case 74: el_order--; // 220-node tetrahedron (9th order)
2047  case 75: el_order--; // 286-node tetrahedron (10th order)
2048  {
2049 #ifdef MFEM_USE_MEMALLOC
2050  elements_3D.push_back(TetMemory.Alloc());
2051  elements_3D.back()->SetVertices(&vert_indices[0]);
2052  elements_3D.back()->SetAttribute(phys_domain);
2053 #else
2054  elements_3D.push_back(
2055  new Tetrahedron(&vert_indices[0], phys_domain));
2056 #endif
2057  if (el_order > 1)
2058  {
2059  Array<int> * hov = new Array<int>;
2060  hov->Append(&vert_indices[0], n_elem_nodes);
2061  ho_verts_3D.push_back(hov);
2062  ho_el_order_3D.push_back(el_order);
2063  }
2064  break;
2065  }
2066  case 5: el_order--; // 8-node hexahedron
2067  case 12: el_order--; // 27-node hexahedron (2nd order)
2068  case 92: el_order--; // 64-node hexahedron (3rd order)
2069  case 93: el_order--; // 125-node hexahedron (4th order)
2070  case 94: el_order--; // 216-node hexahedron (5th order)
2071  case 95: el_order--; // 343-node hexahedron (6th order)
2072  case 96: el_order--; // 512-node hexahedron (7th order)
2073  case 97: el_order--; // 729-node hexahedron (8th order)
2074  case 98: el_order--; // 1000-node hexahedron (9th order)
2075  {
2076  el_order--; // Gmsh does not define an order 10 hex
2077  elements_3D.push_back(
2078  new Hexahedron(&vert_indices[0], phys_domain));
2079  if (el_order > 1)
2080  {
2081  Array<int> * hov = new Array<int>;
2082  hov->Append(&vert_indices[0], n_elem_nodes);
2083  ho_verts_3D.push_back(hov);
2084  ho_el_order_3D.push_back(el_order);
2085  }
2086  break;
2087  }
2088  case 6: el_order--; // 6-node wedge
2089  case 13: el_order--; // 18-node wedge (2nd order)
2090  case 90: el_order--; // 40-node wedge (3rd order)
2091  case 91: el_order--; // 75-node wedge (4th order)
2092  case 106: el_order--; // 126-node wedge (5th order)
2093  case 107: el_order--; // 196-node wedge (6th order)
2094  case 108: el_order--; // 288-node wedge (7th order)
2095  case 109: el_order--; // 405-node wedge (8th order)
2096  case 110: el_order--; // 550-node wedge (9th order)
2097  {
2098  el_order--; // Gmsh does not define an order 10 wedge
2099  elements_3D.push_back(
2100  new Wedge(&vert_indices[0], phys_domain));
2101  if (el_order > 1)
2102  {
2103  Array<int> * hov = new Array<int>;
2104  hov->Append(&vert_indices[0], n_elem_nodes);
2105  ho_verts_3D.push_back(hov);
2106  ho_el_order_3D.push_back(el_order);
2107  }
2108  break;
2109  }
2110  case 7: el_order--; // 5-node pyramid
2111  case 14: el_order--; // 14-node pyramid (2nd order)
2112  case 118: el_order--; // 30-node pyramid (3rd order)
2113  case 119: el_order--; // 55-node pyramid (4th order)
2114  case 120: el_order--; // 91-node pyramid (5th order)
2115  case 121: el_order--; // 140-node pyramid (6th order)
2116  case 122: el_order--; // 204-node pyramid (7th order)
2117  case 123: el_order--; // 285-node pyramid (8th order)
2118  case 124: el_order--; // 385-node pyramid (9th order)
2119  {
2120  el_order--; // Gmsh does not define an order 10 pyr
2121  elements_3D.push_back(
2122  new Pyramid(&vert_indices[0], phys_domain));
2123  if (el_order > 1)
2124  {
2125  Array<int> * hov = new Array<int>;
2126  hov->Append(&vert_indices[0], n_elem_nodes);
2127  ho_verts_3D.push_back(hov);
2128  ho_el_order_3D.push_back(el_order);
2129  }
2130  break;
2131  }
2132  case 15: // 1-node point
2133  {
2134  elements_0D.push_back(
2135  new Point(&vert_indices[0], phys_domain));
2136  break;
2137  }
2138  default: // any other element
2139  MFEM_WARNING("Unsupported Gmsh element type.");
2140  break;
2141 
2142  } // switch (type_of_element)
2143  } // el (elements of one type)
2144  } // all elements
2145  } // if binary
2146  else // ASCII
2147  {
2148  for (int el = 0; el < num_of_all_elements; ++el)
2149  {
2150  input >> serial_number >> type_of_element >> n_tags;
2151  vector<int> data(n_tags);
2152  for (int i = 0; i < n_tags; ++i) { input >> data[i]; }
2153  // physical domain - the most important value (to distinguish
2154  // materials with different properties)
2155  phys_domain = (n_tags > 0) ? data[0] : 1;
2156  // elementary domain - to distinguish different geometrical
2157  // domains (typically, it's used rarely)
2158  elem_domain = (n_tags > 1) ? data[1] : 0;
2159  // the number of tags is bigger than 2 if there are some
2160  // partitions (domain decompositions)
2161  n_partitions = (n_tags > 2) ? data[2] : 0;
2162  // we currently just skip the partitions if they exist, and go
2163  // directly to vertices describing the mesh element
2164  const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
2165  vector<int> vert_indices(n_elem_nodes);
2166  int index;
2167  for (int vi = 0; vi < n_elem_nodes; ++vi)
2168  {
2169  input >> index;
2170  map<int, int>::const_iterator it = vertices_map.find(index);
2171  if (it == vertices_map.end())
2172  {
2173  MFEM_ABORT("Gmsh file : vertex index doesn't exist");
2174  }
2175  vert_indices[vi] = it->second;
2176  }
2177 
2178  // Non-positive attributes are not allowed in MFEM. However,
2179  // by default, Gmsh sets the physical domain of all elements
2180  // to zero. In the case that all elements have physical domain
2181  // zero, we will given them attribute 1. If only some elements
2182  // have physical domain zero, we will throw an error.
2183  if (phys_domain <= 0)
2184  {
2185  has_nonpositive_phys_domain = true;
2186  phys_domain = 1;
2187  }
2188  else
2189  {
2190  has_positive_phys_domain = true;
2191  }
2192 
2193  // initialize the mesh element
2194  int el_order = 11;
2195  switch (type_of_element)
2196  {
2197  case 1: // 2-node line
2198  case 8: // 3-node line (2nd order)
2199  case 26: // 4-node line (3rd order)
2200  case 27: // 5-node line (4th order)
2201  case 28: // 6-node line (5th order)
2202  case 62: // 7-node line (6th order)
2203  case 63: // 8-node line (7th order)
2204  case 64: // 9-node line (8th order)
2205  case 65: // 10-node line (9th order)
2206  case 66: // 11-node line (10th order)
2207  {
2208  elements_1D.push_back(
2209  new Segment(&vert_indices[0], phys_domain));
2210  if (type_of_element != 1)
2211  {
2212  Array<int> * hov = new Array<int>;
2213  hov->Append(&vert_indices[0], n_elem_nodes);
2214  ho_verts_1D.push_back(hov);
2215  el_order = n_elem_nodes - 1;
2216  ho_el_order_1D.push_back(el_order);
2217  }
2218  break;
2219  }
2220  case 2: el_order--; // 3-node triangle
2221  case 9: el_order--; // 6-node triangle (2nd order)
2222  case 21: el_order--; // 10-node triangle (3rd order)
2223  case 23: el_order--; // 15-node triangle (4th order)
2224  case 25: el_order--; // 21-node triangle (5th order)
2225  case 42: el_order--; // 28-node triangle (6th order)
2226  case 43: el_order--; // 36-node triangle (7th order)
2227  case 44: el_order--; // 45-node triangle (8th order)
2228  case 45: el_order--; // 55-node triangle (9th order)
2229  case 46: el_order--; // 66-node triangle (10th order)
2230  {
2231  elements_2D.push_back(
2232  new Triangle(&vert_indices[0], phys_domain));
2233  if (el_order > 1)
2234  {
2235  Array<int> * hov = new Array<int>;
2236  hov->Append(&vert_indices[0], n_elem_nodes);
2237  ho_verts_2D.push_back(hov);
2238  ho_el_order_2D.push_back(el_order);
2239  }
2240  break;
2241  }
2242  case 3: el_order--; // 4-node quadrangle
2243  case 10: el_order--; // 9-node quadrangle (2nd order)
2244  case 36: el_order--; // 16-node quadrangle (3rd order)
2245  case 37: el_order--; // 25-node quadrangle (4th order)
2246  case 38: el_order--; // 36-node quadrangle (5th order)
2247  case 47: el_order--; // 49-node quadrangle (6th order)
2248  case 48: el_order--; // 64-node quadrangle (7th order)
2249  case 49: el_order--; // 81-node quadrangle (8th order)
2250  case 50: el_order--; // 100-node quadrangle (9th order)
2251  case 51: el_order--; // 121-node quadrangle (10th order)
2252  {
2253  elements_2D.push_back(
2254  new Quadrilateral(&vert_indices[0], phys_domain));
2255  if (el_order > 1)
2256  {
2257  Array<int> * hov = new Array<int>;
2258  hov->Append(&vert_indices[0], n_elem_nodes);
2259  ho_verts_2D.push_back(hov);
2260  ho_el_order_2D.push_back(el_order);
2261  }
2262  break;
2263  }
2264  case 4: el_order--; // 4-node tetrahedron
2265  case 11: el_order--; // 10-node tetrahedron (2nd order)
2266  case 29: el_order--; // 20-node tetrahedron (3rd order)
2267  case 30: el_order--; // 35-node tetrahedron (4th order)
2268  case 31: el_order--; // 56-node tetrahedron (5th order)
2269  case 71: el_order--; // 84-node tetrahedron (6th order)
2270  case 72: el_order--; // 120-node tetrahedron (7th order)
2271  case 73: el_order--; // 165-node tetrahedron (8th order)
2272  case 74: el_order--; // 220-node tetrahedron (9th order)
2273  case 75: el_order--; // 286-node tetrahedron (10th order)
2274  {
2275 #ifdef MFEM_USE_MEMALLOC
2276  elements_3D.push_back(TetMemory.Alloc());
2277  elements_3D.back()->SetVertices(&vert_indices[0]);
2278  elements_3D.back()->SetAttribute(phys_domain);
2279 #else
2280  elements_3D.push_back(
2281  new Tetrahedron(&vert_indices[0], phys_domain));
2282 #endif
2283  if (el_order > 1)
2284  {
2285  Array<int> * hov = new Array<int>;
2286  hov->Append(&vert_indices[0], n_elem_nodes);
2287  ho_verts_3D.push_back(hov);
2288  ho_el_order_3D.push_back(el_order);
2289  }
2290  break;
2291  }
2292  case 5: el_order--; // 8-node hexahedron
2293  case 12: el_order--; // 27-node hexahedron (2nd order)
2294  case 92: el_order--; // 64-node hexahedron (3rd order)
2295  case 93: el_order--; // 125-node hexahedron (4th order)
2296  case 94: el_order--; // 216-node hexahedron (5th order)
2297  case 95: el_order--; // 343-node hexahedron (6th order)
2298  case 96: el_order--; // 512-node hexahedron (7th order)
2299  case 97: el_order--; // 729-node hexahedron (8th order)
2300  case 98: el_order--; // 1000-node hexahedron (9th order)
2301  {
2302  el_order--;
2303  elements_3D.push_back(
2304  new Hexahedron(&vert_indices[0], phys_domain));
2305  if (el_order > 1)
2306  {
2307  Array<int> * hov = new Array<int>;
2308  hov->Append(&vert_indices[0], n_elem_nodes);
2309  ho_verts_3D.push_back(hov);
2310  ho_el_order_3D.push_back(el_order);
2311  }
2312  break;
2313  }
2314  case 6: el_order--; // 6-node wedge
2315  case 13: el_order--; // 18-node wedge (2nd order)
2316  case 90: el_order--; // 40-node wedge (3rd order)
2317  case 91: el_order--; // 75-node wedge (4th order)
2318  case 106: el_order--; // 126-node wedge (5th order)
2319  case 107: el_order--; // 196-node wedge (6th order)
2320  case 108: el_order--; // 288-node wedge (7th order)
2321  case 109: el_order--; // 405-node wedge (8th order)
2322  case 110: el_order--; // 550-node wedge (9th order)
2323  {
2324  el_order--;
2325  elements_3D.push_back(
2326  new Wedge(&vert_indices[0], phys_domain));
2327  if (el_order > 1)
2328  {
2329  Array<int> * hov = new Array<int>;
2330  hov->Append(&vert_indices[0], n_elem_nodes);
2331  ho_verts_3D.push_back(hov);
2332  ho_el_order_3D.push_back(el_order);
2333  }
2334  break;
2335  }
2336  case 7: el_order--; // 5-node pyramid
2337  case 14: el_order--; // 14-node pyramid (2nd order)
2338  case 118: el_order--; // 30-node pyramid (3rd order)
2339  case 119: el_order--; // 55-node pyramid (4th order)
2340  case 120: el_order--; // 91-node pyramid (5th order)
2341  case 121: el_order--; // 140-node pyramid (6th order)
2342  case 122: el_order--; // 204-node pyramid (7th order)
2343  case 123: el_order--; // 285-node pyramid (8th order)
2344  case 124: el_order--; // 385-node pyramid (9th order)
2345  {
2346  el_order--;
2347  elements_3D.push_back(
2348  new Pyramid(&vert_indices[0], phys_domain));
2349  if (el_order > 1)
2350  {
2351  Array<int> * hov = new Array<int>;
2352  hov->Append(&vert_indices[0], n_elem_nodes);
2353  ho_verts_3D.push_back(hov);
2354  ho_el_order_3D.push_back(el_order);
2355  }
2356  break;
2357  }
2358  case 15: // 1-node point
2359  {
2360  elements_0D.push_back(
2361  new Point(&vert_indices[0], phys_domain));
2362  break;
2363  }
2364  default: // any other element
2365  MFEM_WARNING("Unsupported Gmsh element type.");
2366  break;
2367 
2368  } // switch (type_of_element)
2369  } // el (all elements)
2370  } // if ASCII
2371 
2372  if (has_positive_phys_domain && has_nonpositive_phys_domain)
2373  {
2374  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
2375  "By default Gmsh sets element tags (attributes)"
2376  " to '0' but MFEM requires that they be"
2377  " positive integers.\n"
2378  "Use \"Physical Curve\", \"Physical Surface\","
2379  " or \"Physical Volume\" to set tags/attributes"
2380  " for all curves, surfaces, or volumes in your"
2381  " Gmsh geometry to values which are >= 1.");
2382  }
2383  else if (has_nonpositive_phys_domain)
2384  {
2385  mfem::out << "\nGmsh reader: all element attributes were zero.\n"
2386  << "MFEM only supports positive element attributes.\n"
2387  << "Setting element attributes to 1.\n\n";
2388  }
2389 
2390  if (!elements_3D.empty())
2391  {
2392  Dim = 3;
2393  NumOfElements = elements_3D.size();
2394  elements.SetSize(NumOfElements);
2395  for (int el = 0; el < NumOfElements; ++el)
2396  {
2397  elements[el] = elements_3D[el];
2398  }
2399  NumOfBdrElements = elements_2D.size();
2400  boundary.SetSize(NumOfBdrElements);
2401  for (int el = 0; el < NumOfBdrElements; ++el)
2402  {
2403  boundary[el] = elements_2D[el];
2404  }
2405  for (size_t el = 0; el < ho_el_order_3D.size(); el++)
2406  {
2407  mesh_order = max(mesh_order, ho_el_order_3D[el]);
2408  }
2409  // discard other elements
2410  for (size_t el = 0; el < elements_1D.size(); ++el)
2411  {
2412  delete elements_1D[el];
2413  }
2414  for (size_t el = 0; el < elements_0D.size(); ++el)
2415  {
2416  delete elements_0D[el];
2417  }
2418  }
2419  else if (!elements_2D.empty())
2420  {
2421  Dim = 2;
2422  NumOfElements = elements_2D.size();
2423  elements.SetSize(NumOfElements);
2424  for (int el = 0; el < NumOfElements; ++el)
2425  {
2426  elements[el] = elements_2D[el];
2427  }
2428  NumOfBdrElements = elements_1D.size();
2429  boundary.SetSize(NumOfBdrElements);
2430  for (int el = 0; el < NumOfBdrElements; ++el)
2431  {
2432  boundary[el] = elements_1D[el];
2433  }
2434  for (size_t el = 0; el < ho_el_order_2D.size(); el++)
2435  {
2436  mesh_order = max(mesh_order, ho_el_order_2D[el]);
2437  }
2438  // discard other elements
2439  for (size_t el = 0; el < elements_0D.size(); ++el)
2440  {
2441  delete elements_0D[el];
2442  }
2443  }
2444  else if (!elements_1D.empty())
2445  {
2446  Dim = 1;
2447  NumOfElements = elements_1D.size();
2448  elements.SetSize(NumOfElements);
2449  for (int el = 0; el < NumOfElements; ++el)
2450  {
2451  elements[el] = elements_1D[el];
2452  }
2453  NumOfBdrElements = elements_0D.size();
2454  boundary.SetSize(NumOfBdrElements);
2455  for (int el = 0; el < NumOfBdrElements; ++el)
2456  {
2457  boundary[el] = elements_0D[el];
2458  }
2459  for (size_t el = 0; el < ho_el_order_1D.size(); el++)
2460  {
2461  mesh_order = max(mesh_order, ho_el_order_1D[el]);
2462  }
2463  }
2464  else
2465  {
2466  MFEM_ABORT("Gmsh file : no elements found");
2467  return;
2468  }
2469 
2470  if (mesh_order > 1)
2471  {
2472  curved = 1;
2473  read_gf = 0;
2474 
2475  // initialize mesh_geoms so we can create Nodes FE space below
2476  this->SetMeshGen();
2477 
2478  // Generate faces and edges so that we can define
2479  // FE space on the mesh
2480  this->FinalizeTopology();
2481 
2482  // Construct GridFunction for uniformly spaced high order coords
2484  FiniteElementSpace* nfes;
2485  nfec = new L2_FECollection(mesh_order, Dim,
2486  BasisType::ClosedUniform);
2487  nfes = new FiniteElementSpace(this, nfec, spaceDim,
2488  Ordering::byVDIM);
2489  Nodes_gf.SetSpace(nfes);
2490  Nodes_gf.MakeOwner(nfec);
2491 
2492  int o = 0;
2493  int el_order = 1;
2494  for (int el = 0; el < NumOfElements; el++)
2495  {
2496  const int * vm = NULL;
2497  Array<int> * ho_verts = NULL;
2498  switch (GetElementType(el))
2499  {
2500  case Element::SEGMENT:
2501  ho_verts = ho_verts_1D[el];
2502  el_order = ho_el_order_1D[el];
2503  if (!ho_lin[el_order])
2504  {
2505  ho_lin[el_order] = new int[ho_verts->Size()];
2506  GmshHOSegmentMapping(el_order, ho_lin[el_order]);
2507  }
2508  vm = ho_lin[el_order];
2509  break;
2510  case Element::TRIANGLE:
2511  ho_verts = ho_verts_2D[el];
2512  el_order = ho_el_order_2D[el];
2513  if (!ho_tri[el_order])
2514  {
2515  ho_tri[el_order] = new int[ho_verts->Size()];
2516  GmshHOTriangleMapping(el_order, ho_tri[el_order]);
2517  }
2518  vm = ho_tri[el_order];
2519  break;
2520  case Element::QUADRILATERAL:
2521  ho_verts = ho_verts_2D[el];
2522  el_order = ho_el_order_2D[el];
2523  if (!ho_sqr[el_order])
2524  {
2525  ho_sqr[el_order] = new int[ho_verts->Size()];
2526  GmshHOQuadrilateralMapping(el_order, ho_sqr[el_order]);
2527  }
2528  vm = ho_sqr[el_order];
2529  break;
2530  case Element::TETRAHEDRON:
2531  ho_verts = ho_verts_3D[el];
2532  el_order = ho_el_order_3D[el];
2533  if (!ho_tet[el_order])
2534  {
2535  ho_tet[el_order] = new int[ho_verts->Size()];
2536  GmshHOTetrahedronMapping(el_order, ho_tet[el_order]);
2537  }
2538  vm = ho_tet[el_order];
2539  break;
2540  case Element::HEXAHEDRON:
2541  ho_verts = ho_verts_3D[el];
2542  el_order = ho_el_order_3D[el];
2543  if (!ho_hex[el_order])
2544  {
2545  ho_hex[el_order] = new int[ho_verts->Size()];
2546  GmshHOHexahedronMapping(el_order, ho_hex[el_order]);
2547  }
2548  vm = ho_hex[el_order];
2549  break;
2550  case Element::WEDGE:
2551  ho_verts = ho_verts_3D[el];
2552  el_order = ho_el_order_3D[el];
2553  if (!ho_wdg[el_order])
2554  {
2555  ho_wdg[el_order] = new int[ho_verts->Size()];
2556  GmshHOWedgeMapping(el_order, ho_wdg[el_order]);
2557  }
2558  vm = ho_wdg[el_order];
2559  break;
2560  case Element::PYRAMID:
2561  ho_verts = ho_verts_3D[el];
2562  el_order = ho_el_order_3D[el];
2563  if (!ho_pyr[el_order])
2564  {
2565  ho_pyr[el_order] = new int[ho_verts->Size()];
2566  GmshHOPyramidMapping(el_order, ho_pyr[el_order]);
2567  }
2568  vm = ho_pyr[el_order];
2569  break;
2570  default: // Any other element type
2571  MFEM_WARNING("Unsupported Gmsh element type.");
2572  break;
2573  }
2574  int nv = (ho_verts) ? ho_verts->Size() : 0;
2575 
2576  for (int v = 0; v<nv; v++)
2577  {
2578  double * c = GetVertex((*ho_verts)[vm[v]]);
2579  for (int d=0; d<spaceDim; d++)
2580  {
2581  Nodes_gf(spaceDim * (o + v) + d) = c[d];
2582  }
2583  }
2584  o += nv;
2585  }
2586  }
2587 
2588  // Delete any high order element to vertex connectivity
2589  for (size_t el=0; el<ho_verts_1D.size(); el++)
2590  {
2591  delete ho_verts_1D[el];
2592  }
2593  for (size_t el=0; el<ho_verts_2D.size(); el++)
2594  {
2595  delete ho_verts_2D[el];
2596  }
2597  for (size_t el=0; el<ho_verts_3D.size(); el++)
2598  {
2599  delete ho_verts_3D[el];
2600  }
2601 
2602  // Delete dynamically allocated high vertex order mappings
2603  for (int ord=4; ord<ho_lin.Size(); ord++)
2604  {
2605  if (ho_lin[ord] != NULL) { delete [] ho_lin[ord]; }
2606  }
2607  for (int ord=4; ord<ho_tri.Size(); ord++)
2608  {
2609  if (ho_tri[ord] != NULL) { delete [] ho_tri[ord]; }
2610  }
2611  for (int ord=4; ord<ho_sqr.Size(); ord++)
2612  {
2613  if (ho_sqr[ord] != NULL) { delete [] ho_sqr[ord]; }
2614  }
2615  for (int ord=4; ord<ho_tet.Size(); ord++)
2616  {
2617  if (ho_tet[ord] != NULL) { delete [] ho_tet[ord]; }
2618  }
2619  for (int ord=4; ord<ho_hex.Size(); ord++)
2620  {
2621  if (ho_hex[ord] != NULL) { delete [] ho_hex[ord]; }
2622  }
2623  for (int ord=4; ord<ho_wdg.Size(); ord++)
2624  {
2625  if (ho_wdg[ord] != NULL) { delete [] ho_wdg[ord]; }
2626  }
2627  for (int ord=4; ord<ho_pyr.Size(); ord++)
2628  {
2629  if (ho_pyr[ord] != NULL) { delete [] ho_pyr[ord]; }
2630  }
2631 
2632  // Suppress warnings (MFEM_CONTRACT_VAR does not work here with nvcc):
2633  ++n_partitions;
2634  ++elem_domain;
2635 
2636  } // section '$Elements'
2637  else if (buff == "$Periodic") // Reading master/slave node pairs
2638  {
2639  curved = 1;
2640  read_gf = 0;
2641  periodic = true;
2642 
2643  Array<int> v2v(NumOfVertices);
2644  for (int i = 0; i < v2v.Size(); i++)
2645  {
2646  v2v[i] = i;
2647  }
2648  int num_per_ent;
2649  int num_nodes;
2650  input >> num_per_ent;
2651  getline(input, buff); // Read end-of-line
2652  for (int i = 0; i < num_per_ent; i++)
2653  {
2654  getline(input, buff); // Read and ignore entity dimension and tags
2655  getline(input, buff); // If affine mapping exist, read and ignore
2656  if (!strncmp(buff.c_str(), "Affine", 6))
2657  {
2658  input >> num_nodes;
2659  }
2660  else
2661  {
2662  num_nodes = atoi(buff.c_str());
2663  }
2664  for (int j=0; j<num_nodes; j++)
2665  {
2666  int slave, master;
2667  input >> slave >> master;
2668  v2v[slave - 1] = master - 1;
2669  }
2670  getline(input, buff); // Read end-of-line
2671  }
2672 
2673  // Follow existing long chains of slave->master in v2v array.
2674  // Upon completion of this loop, each v2v[slave] will point to a true
2675  // master vertex. This algorithm is useful for periodicity defined in
2676  // multiple directions.
2677  for (int slave = 0; slave < v2v.Size(); slave++)
2678  {
2679  int master = v2v[slave];
2680  if (master != slave)
2681  {
2682  // This loop will end if it finds a circular dependency.
2683  while (v2v[master] != master && master != slave)
2684  {
2685  master = v2v[master];
2686  }
2687  if (master == slave)
2688  {
2689  // if master and slave are the same vertex, circular dependency
2690  // exists. We need to fix the problem, we choose slave.
2691  v2v[slave] = slave;
2692  }
2693  else
2694  {
2695  // the long chain has ended on the true master vertex.
2696  v2v[slave] = master;
2697  }
2698  }
2699  }
2700 
2701  // Convert nodes to discontinuous GridFunction (if they aren't already)
2702  if (mesh_order == 1)
2703  {
2704  this->FinalizeTopology();
2705  this->SetMeshGen();
2706  this->SetCurvature(1, true, spaceDim, Ordering::byVDIM);
2707  }
2708 
2709  // Replace "slave" vertex indices in the element connectivity
2710  // with their corresponding "master" vertex indices.
2711  for (int i = 0; i < this->GetNE(); i++)
2712  {
2713  Element *el = this->GetElement(i);
2714  int *v = el->GetVertices();
2715  int nv = el->GetNVertices();
2716  for (int j = 0; j < nv; j++)
2717  {
2718  v[j] = v2v[v[j]];
2719  }
2720  }
2721  // Replace "slave" vertex indices in the boundary element connectivity
2722  // with their corresponding "master" vertex indices.
2723  for (int i = 0; i < this->GetNBE(); i++)
2724  {
2725  Element *el = this->GetBdrElement(i);
2726  int *v = el->GetVertices();
2727  int nv = el->GetNVertices();
2728  for (int j = 0; j < nv; j++)
2729  {
2730  v[j] = v2v[v[j]];
2731  }
2732  }
2733  }
2734  } // we reach the end of the file
2735 
2736  this->RemoveUnusedVertices();
2737  if (periodic)
2738  {
2739  this->RemoveInternalBoundaries();
2740  }
2741  this->FinalizeTopology();
2742 
2743  // If a high order coordinate field was created project it onto the mesh
2744  if (mesh_order > 1)
2745  {
2746  SetCurvature(mesh_order, periodic, spaceDim, Ordering::byVDIM);
2747 
2748  VectorGridFunctionCoefficient NodesCoef(&Nodes_gf);
2749  Nodes->ProjectCoefficient(NodesCoef);
2750  }
2751 }
2752 
2753 
2754 #ifdef MFEM_USE_NETCDF
2755 void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
2756 {
2757  read_gf = 0;
2758 
2759  // curved set to zero will change if mesh is indeed curved
2760  curved = 0;
2761 
2762  const int sideMapTri3[3][2] =
2763  {
2764  {1,2},
2765  {2,3},
2766  {3,1},
2767  };
2768 
2769  const int sideMapQuad4[4][2] =
2770  {
2771  {1,2},
2772  {2,3},
2773  {3,4},
2774  {4,1},
2775  };
2776 
2777  const int sideMapTri6[3][3] =
2778  {
2779  {1,2,4},
2780  {2,3,5},
2781  {3,1,6},
2782  };
2783 
2784  const int sideMapQuad9[4][3] =
2785  {
2786  {1,2,5},
2787  {2,3,6},
2788  {3,4,7},
2789  {4,1,8},
2790  };
2791 
2792  const int sideMapTet4[4][3] =
2793  {
2794  {1,2,4},
2795  {2,3,4},
2796  {1,4,3},
2797  {1,3,2}
2798  };
2799 
2800  const int sideMapTet10[4][6] =
2801  {
2802  {1,2,4,5,9,8},
2803  {2,3,4,6,10,9},
2804  {1,4,3,8,10,7},
2805  {1,3,2,7,6,5}
2806  };
2807 
2808  const int sideMapHex8[6][4] =
2809  {
2810  {1,2,6,5},
2811  {2,3,7,6},
2812  {4,3,7,8},
2813  {1,4,8,5},
2814  {1,4,3,2},
2815  {5,8,7,6}
2816  };
2817 
2818  const int sideMapHex27[6][9] =
2819  {
2820  {1,2,6,5,9,14,17,13,26},
2821  {2,3,7,6,10,15,18,14,25},
2822  {4,3,7,8,11,15,19,16,27},
2823  {1,4,8,5,12,16,20,13,24},
2824  {1,4,3,2,12,11,10,9,22},
2825  {5,8,7,6,20,19,18,17,23}
2826  };
2827 
2828 
2829  // 1,2,3,4,5,6,7,8,9,10
2830  const int mfemToGenesisTet10[10] = {1,2,3,4,5,7,8,6,9,10};
2831 
2832  // 1,2,3,4,5,6,7,8,9,10,11,
2833  const int mfemToGenesisHex27[27] = {1,2,3,4,5,6,7,8,9,10,11,
2834  // 12,13,14,15,16,17,18,19
2835  12,17,18,19,20,13,14,15,
2836  // 20,21,22,23,24,25,26,27
2837  16,22,26,25,27,24,23,21
2838  };
2839 
2840  const int mfemToGenesisTri6[6] = {1,2,3,4,5,6};
2841  const int mfemToGenesisQuad9[9] = {1,2,3,4,5,6,7,8,9};
2842 
2843 
2844  // error handling.
2845  int retval;
2846 
2847  // dummy string
2848  constexpr size_t buf_size = 256;
2849  char str_dummy[buf_size];
2850 
2851  char temp_str[buf_size];
2852  int temp_id;
2853 
2854  // open the file.
2855  int ncid;
2856  if ((retval = nc_open(filename, NC_NOWRITE, &ncid)))
2857  {
2858  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
2859  }
2860 
2861  // read important dimensions
2862 
2863  int id;
2864  size_t num_dim=0, num_nodes=0, num_elem=0, num_el_blk=0, num_side_sets=0;
2865 
2866  if ((retval = nc_inq_dimid(ncid, "num_dim", &id)) ||
2867  (retval = nc_inq_dim(ncid, id, str_dummy, &num_dim)) ||
2868 
2869  (retval = nc_inq_dimid(ncid, "num_nodes", &id)) ||
2870  (retval = nc_inq_dim(ncid, id, str_dummy, &num_nodes)) ||
2871 
2872  (retval = nc_inq_dimid(ncid, "num_elem", &id)) ||
2873  (retval = nc_inq_dim(ncid, id, str_dummy, &num_elem)) ||
2874 
2875  (retval = nc_inq_dimid(ncid, "num_el_blk", &id)) ||
2876  (retval = nc_inq_dim(ncid, id, str_dummy, &num_el_blk)))
2877  {
2878  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
2879  }
2880  if ((retval = nc_inq_dimid(ncid, "num_side_sets", &id)) ||
2881  (retval = nc_inq_dim(ncid, id, str_dummy, &num_side_sets)))
2882  {
2883  num_side_sets = 0;
2884  }
2885 
2886  Dim = num_dim;
2887 
2888  // create arrays for element blocks
2889  size_t *num_el_in_blk = new size_t[num_el_blk];
2890  size_t num_node_per_el;
2891 
2892  int previous_num_node_per_el = 0;
2893  for (int i = 0; i < (int) num_el_blk; i++)
2894  {
2895  snprintf(temp_str, buf_size, "num_el_in_blk%d", i+1);
2896  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
2897  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_el_in_blk[i])))
2898  {
2899  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
2900  }
2901 
2902  snprintf(temp_str, buf_size, "num_nod_per_el%d", i+1);
2903  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
2904  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_node_per_el)))
2905  {
2906  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
2907  }
2908 
2909  // check for different element types in each block
2910  // which is not currently supported
2911  if (i != 0)
2912  {
2913  if ((int) num_node_per_el != previous_num_node_per_el)
2914  {
2915  MFEM_ABORT("Element blocks of different element types not supported");
2916  }
2917  }
2918  previous_num_node_per_el = num_node_per_el;
2919  }
2920 
2921  // Determine CUBIT element and face type
2922  enum CubitElementType
2923  {
2924  ELEMENT_TRI3,
2925  ELEMENT_TRI6,
2926  ELEMENT_QUAD4,
2927  ELEMENT_QUAD9,
2928  ELEMENT_TET4,
2929  ELEMENT_TET10,
2930  ELEMENT_HEX8,
2931  ELEMENT_HEX27
2932  };
2933 
2934  enum CubitFaceType
2935  {
2936  FACE_EDGE2,
2937  FACE_EDGE3,
2938  FACE_TRI3,
2939  FACE_TRI6,
2940  FACE_QUAD4,
2941  FACE_QUAD9
2942  };
2943 
2944  CubitElementType cubit_element_type = ELEMENT_TRI3; // suppress a warning
2945  CubitFaceType cubit_face_type = FACE_EDGE2; // suppress a warning
2946  int num_element_linear_nodes = 0; // initialize to suppress a warning
2947 
2948  if (num_dim == 2)
2949  {
2950  switch (num_node_per_el)
2951  {
2952  case (3) :
2953  {
2954  cubit_element_type = ELEMENT_TRI3;
2955  cubit_face_type = FACE_EDGE2;
2956  num_element_linear_nodes = 3;
2957  break;
2958  }
2959  case (6) :
2960  {
2961  cubit_element_type = ELEMENT_TRI6;
2962  cubit_face_type = FACE_EDGE3;
2963  num_element_linear_nodes = 3;
2964  break;
2965  }
2966  case (4) :
2967  {
2968  cubit_element_type = ELEMENT_QUAD4;
2969  cubit_face_type = FACE_EDGE2;
2970  num_element_linear_nodes = 4;
2971  break;
2972  }
2973  case (9) :
2974  {
2975  cubit_element_type = ELEMENT_QUAD9;
2976  cubit_face_type = FACE_EDGE3;
2977  num_element_linear_nodes = 4;
2978  break;
2979  }
2980  default :
2981  {
2982  MFEM_ABORT("Don't know what to do with a " << num_node_per_el <<
2983  " node 2D element\n");
2984  }
2985  }
2986  }
2987  else if (num_dim == 3)
2988  {
2989  switch (num_node_per_el)
2990  {
2991  case (4) :
2992  {
2993  cubit_element_type = ELEMENT_TET4;
2994  cubit_face_type = FACE_TRI3;
2995  num_element_linear_nodes = 4;
2996  break;
2997  }
2998  case (10) :
2999  {
3000  cubit_element_type = ELEMENT_TET10;
3001  cubit_face_type = FACE_TRI6;
3002  num_element_linear_nodes = 4;
3003  break;
3004  }
3005  case (8) :
3006  {
3007  cubit_element_type = ELEMENT_HEX8;
3008  cubit_face_type = FACE_QUAD4;
3009  num_element_linear_nodes = 8;
3010  break;
3011  }
3012  case (27) :
3013  {
3014  cubit_element_type = ELEMENT_HEX27;
3015  cubit_face_type = FACE_QUAD9;
3016  num_element_linear_nodes = 8;
3017  break;
3018  }
3019  default :
3020  {
3021  MFEM_ABORT("Don't know what to do with a " << num_node_per_el <<
3022  " node 3D element\n");
3023  }
3024  }
3025  }
3026  else
3027  {
3028  MFEM_ABORT("Invalid dimension: num_dim = " << num_dim);
3029  }
3030 
3031  // Determine order of elements
3032  int order = 0;
3033  if (cubit_element_type == ELEMENT_TRI3 || cubit_element_type == ELEMENT_QUAD4 ||
3034  cubit_element_type == ELEMENT_TET4 || cubit_element_type == ELEMENT_HEX8)
3035  {
3036  order = 1;
3037  }
3038  else if (cubit_element_type == ELEMENT_TRI6 ||
3039  cubit_element_type == ELEMENT_QUAD9 ||
3040  cubit_element_type == ELEMENT_TET10 || cubit_element_type == ELEMENT_HEX27)
3041  {
3042  order = 2;
3043  }
3044 
3045  // create array for number of sides in side sets
3046  size_t *num_side_in_ss = new size_t[num_side_sets];
3047  for (int i = 0; i < (int) num_side_sets; i++)
3048  {
3049  snprintf(temp_str, buf_size, "num_side_ss%d", i+1);
3050  if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
3051  (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_side_in_ss[i])))
3052  {
3053  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3054  }
3055  }
3056 
3057  // read the coordinates
3058  double *coordx = new double[num_nodes];
3059  double *coordy = new double[num_nodes];
3060  double *coordz = new double[num_nodes];
3061 
3062  if ((retval = nc_inq_varid(ncid, "coordx", &id)) ||
3063  (retval = nc_get_var_double(ncid, id, coordx)) ||
3064  (retval = nc_inq_varid(ncid, "coordy", &id)) ||
3065  (retval = nc_get_var_double(ncid, id, coordy)))
3066  {
3067  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3068  }
3069 
3070  if (num_dim == 3)
3071  {
3072  if ((retval = nc_inq_varid(ncid, "coordz", &id)) ||
3073  (retval = nc_get_var_double(ncid, id, coordz)))
3074  {
3075  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3076  }
3077  }
3078 
3079  // read the element blocks
3080  int **elem_blk = new int*[num_el_blk];
3081  for (int i = 0; i < (int) num_el_blk; i++)
3082  {
3083  elem_blk[i] = new int[num_el_in_blk[i] * num_node_per_el];
3084  snprintf(temp_str, buf_size, "connect%d", i+1);
3085  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
3086  (retval = nc_get_var_int(ncid, temp_id, elem_blk[i])))
3087  {
3088  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3089  }
3090  }
3091  int *ebprop = new int[num_el_blk];
3092  if ((retval = nc_inq_varid(ncid, "eb_prop1", &id)) ||
3093  (retval = nc_get_var_int(ncid, id, ebprop)))
3094  {
3095  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3096  }
3097 
3098  // read the side sets, a side is is given by (element, face) pairs
3099 
3100  int **elem_ss = new int*[num_side_sets];
3101  int **side_ss = new int*[num_side_sets];
3102 
3103  for (int i = 0; i < (int) num_side_sets; i++)
3104  {
3105  elem_ss[i] = new int[num_side_in_ss[i]];
3106  side_ss[i] = new int[num_side_in_ss[i]];
3107 
3108  snprintf(temp_str, buf_size, "elem_ss%d", i+1);
3109  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
3110  (retval = nc_get_var_int(ncid, temp_id, elem_ss[i])))
3111  {
3112  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3113  }
3114 
3115  snprintf(temp_str, buf_size,"side_ss%d",i+1);
3116  if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
3117  (retval = nc_get_var_int(ncid, temp_id, side_ss[i])))
3118  {
3119  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3120  }
3121  }
3122 
3123  int *ssprop = new int[num_side_sets];
3124  if ((num_side_sets > 0) &&
3125  ((retval = nc_inq_varid(ncid, "ss_prop1", &id)) ||
3126  (retval = nc_get_var_int(ncid, id, ssprop))))
3127  {
3128  MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
3129  }
3130 
3131  // convert (elem,side) pairs to 2D elements
3132 
3133 
3134  int num_face_nodes = 0;
3135  int num_face_linear_nodes = 0;
3136 
3137  switch (cubit_face_type)
3138  {
3139  case (FACE_EDGE2):
3140  {
3141  num_face_nodes = 2;
3142  num_face_linear_nodes = 2;
3143  break;
3144  }
3145  case (FACE_EDGE3):
3146  {
3147  num_face_nodes = 3;
3148  num_face_linear_nodes = 2;
3149  break;
3150  }
3151  case (FACE_TRI3):
3152  {
3153  num_face_nodes = 3;
3154  num_face_linear_nodes = 3;
3155  break;
3156  }
3157  case (FACE_TRI6):
3158  {
3159  num_face_nodes = 6;
3160  num_face_linear_nodes = 3;
3161  break;
3162  }
3163  case (FACE_QUAD4):
3164  {
3165  num_face_nodes = 4;
3166  num_face_linear_nodes = 4;
3167  break;
3168  }
3169  case (FACE_QUAD9):
3170  {
3171  num_face_nodes = 9;
3172  num_face_linear_nodes = 4;
3173  break;
3174  }
3175  }
3176 
3177  // given a global element number, determine the element block and local
3178  // element number
3179  int *start_of_block = new int[num_el_blk+1];
3180  start_of_block[0] = 0;
3181  for (int i = 1; i < (int) num_el_blk+1; i++)
3182  {
3183  start_of_block[i] = start_of_block[i-1] + num_el_in_blk[i-1];
3184  }
3185 
3186  int **ss_node_id = new int*[num_side_sets];
3187 
3188  for (int i = 0; i < (int) num_side_sets; i++)
3189  {
3190  ss_node_id[i] = new int[num_side_in_ss[i]*num_face_nodes];
3191  for (int j = 0; j < (int) num_side_in_ss[i]; j++)
3192  {
3193  int glob_ind = elem_ss[i][j]-1;
3194  int iblk = 0;
3195  int loc_ind;
3196  while (iblk < (int) num_el_blk && glob_ind >= start_of_block[iblk+1])
3197  {
3198  iblk++;
3199  }
3200  if (iblk >= (int) num_el_blk)
3201  {
3202  MFEM_ABORT("Sideset element does not exist");
3203  }
3204  loc_ind = glob_ind - start_of_block[iblk];
3205  int this_side = side_ss[i][j];
3206  int ielem = loc_ind*num_node_per_el;
3207 
3208  for (int k = 0; k < num_face_nodes; k++)
3209  {
3210  int inode;
3211  switch (cubit_element_type)
3212  {
3213  case (ELEMENT_TRI3):
3214  {
3215  inode = sideMapTri3[this_side-1][k];
3216  break;
3217  }
3218  case (ELEMENT_TRI6):
3219  {
3220  inode = sideMapTri6[this_side-1][k];
3221  break;
3222  }
3223  case (ELEMENT_QUAD4):
3224  {
3225  inode = sideMapQuad4[this_side-1][k];
3226  break;
3227  }
3228  case (ELEMENT_QUAD9):
3229  {
3230  inode = sideMapQuad9[this_side-1][k];
3231  break;
3232  }
3233  case (ELEMENT_TET4):
3234  {
3235  inode = sideMapTet4[this_side-1][k];
3236  break;
3237  }
3238  case (ELEMENT_TET10):
3239  {
3240  inode = sideMapTet10[this_side-1][k];
3241  break;
3242  }
3243  case (ELEMENT_HEX8):
3244  {
3245  inode = sideMapHex8[this_side-1][k];
3246  break;
3247  }
3248  case (ELEMENT_HEX27):
3249  {
3250  inode = sideMapHex27[this_side-1][k];
3251  break;
3252  }
3253  }
3254  ss_node_id[i][j*num_face_nodes+k] =
3255  elem_blk[iblk][ielem + inode - 1];
3256  }
3257  }
3258  }
3259 
3260  // we need another node ID mapping since MFEM needs contiguous vertex IDs
3261  std::vector<int> uniqueVertexID;
3262 
3263  for (int iblk = 0; iblk < (int) num_el_blk; iblk++)
3264  {
3265  for (int i = 0; i < (int) num_el_in_blk[iblk]; i++)
3266  {
3267  for (int j = 0; j < num_element_linear_nodes; j++)
3268  {
3269  uniqueVertexID.push_back(elem_blk[iblk][i*num_node_per_el + j]);
3270  }
3271  }
3272  }
3273  std::sort(uniqueVertexID.begin(), uniqueVertexID.end());
3274  std::vector<int>::iterator newEnd;
3275  newEnd = std::unique(uniqueVertexID.begin(), uniqueVertexID.end());
3276  uniqueVertexID.resize(std::distance(uniqueVertexID.begin(), newEnd));
3277 
3278  // OK at this point uniqueVertexID contains a list of all the nodes that are
3279  // actually used by the mesh, 1-based, and sorted. We need to invert this
3280  // list, the inverse is a map
3281 
3282  std::map<int,int> cubitToMFEMVertMap;
3283  for (int i = 0; i < (int) uniqueVertexID.size(); i++)
3284  {
3285  cubitToMFEMVertMap[uniqueVertexID[i]] = i+1;
3286  }
3287  MFEM_ASSERT(cubitToMFEMVertMap.size() == uniqueVertexID.size(),
3288  "This should never happen\n");
3289 
3290  // OK now load up the MFEM mesh structures
3291 
3292  // load up the vertices
3293 
3294  NumOfVertices = uniqueVertexID.size();
3295  vertices.SetSize(NumOfVertices);
3296  for (int i = 0; i < (int) uniqueVertexID.size(); i++)
3297  {
3298  vertices[i](0) = coordx[uniqueVertexID[i] - 1];
3299  vertices[i](1) = coordy[uniqueVertexID[i] - 1];
3300  if (Dim == 3)
3301  {
3302  vertices[i](2) = coordz[uniqueVertexID[i] - 1];
3303  }
3304  }
3305 
3306  NumOfElements = num_elem;
3307  elements.SetSize(num_elem);
3308  int elcount = 0;
3309  int renumberedVertID[8];
3310  for (int iblk = 0; iblk < (int) num_el_blk; iblk++)
3311  {
3312  int NumNodePerEl = num_node_per_el;
3313  for (int i = 0; i < (int) num_el_in_blk[iblk]; i++)
3314  {
3315  for (int j = 0; j < num_element_linear_nodes; j++)
3316  {
3317  renumberedVertID[j] =
3318  cubitToMFEMVertMap[elem_blk[iblk][i*NumNodePerEl+j]]-1;
3319  }
3320 
3321  switch (cubit_element_type)
3322  {
3323  case (ELEMENT_TRI3):
3324  case (ELEMENT_TRI6):
3325  {
3326  elements[elcount] = new Triangle(renumberedVertID,ebprop[iblk]);
3327  break;
3328  }
3329  case (ELEMENT_QUAD4):
3330  case (ELEMENT_QUAD9):
3331  {
3332  elements[elcount] = new Quadrilateral(renumberedVertID,ebprop[iblk]);
3333  break;
3334  }
3335  case (ELEMENT_TET4):
3336  case (ELEMENT_TET10):
3337  {
3338 #ifdef MFEM_USE_MEMALLOC
3339  elements[elcount] = TetMemory.Alloc();
3340  elements[elcount]->SetVertices(renumberedVertID);
3341  elements[elcount]->SetAttribute(ebprop[iblk]);
3342 #else
3343  elements[elcount] = new Tetrahedron(renumberedVertID,
3344  ebprop[iblk]);
3345 #endif
3346  break;
3347  }
3348  case (ELEMENT_HEX8):
3349  case (ELEMENT_HEX27):
3350  {
3351  elements[elcount] = new Hexahedron(renumberedVertID,ebprop[iblk]);
3352  break;
3353  }
3354  }
3355  elcount++;
3356  }
3357  }
3358 
3359  // load up the boundary elements
3360 
3361  NumOfBdrElements = 0;
3362  for (int iss = 0; iss < (int) num_side_sets; iss++)
3363  {
3364  NumOfBdrElements += num_side_in_ss[iss];
3365  }
3366  boundary.SetSize(NumOfBdrElements);
3367  int sidecount = 0;
3368  for (int iss = 0; iss < (int) num_side_sets; iss++)
3369  {
3370  for (int i = 0; i < (int) num_side_in_ss[iss]; i++)
3371  {
3372  for (int j = 0; j < num_face_linear_nodes; j++)
3373  {
3374  renumberedVertID[j] =
3375  cubitToMFEMVertMap[ss_node_id[iss][i*num_face_nodes+j]] - 1;
3376  }
3377  switch (cubit_face_type)
3378  {
3379  case (FACE_EDGE2):
3380  case (FACE_EDGE3):
3381  {
3382  boundary[sidecount] = new Segment(renumberedVertID,ssprop[iss]);
3383  break;
3384  }
3385  case (FACE_TRI3):
3386  case (FACE_TRI6):
3387  {
3388  boundary[sidecount] = new Triangle(renumberedVertID,ssprop[iss]);
3389  break;
3390  }
3391  case (FACE_QUAD4):
3392  case (FACE_QUAD9):
3393  {
3394  boundary[sidecount] = new Quadrilateral(renumberedVertID,ssprop[iss]);
3395  break;
3396  }
3397  }
3398  sidecount++;
3399  }
3400  }
3401 
3402  if (order == 2)
3403  {
3404  curved = 1;
3405  int *mymap = NULL;
3406 
3407  switch (cubit_element_type)
3408  {
3409  case (ELEMENT_TRI6):
3410  {
3411  mymap = (int *) mfemToGenesisTri6;
3412  break;
3413  }
3414  case (ELEMENT_QUAD9):
3415  {
3416  mymap = (int *) mfemToGenesisQuad9;
3417  break;
3418  }
3419  case (ELEMENT_TET10):
3420  {
3421  mymap = (int *) mfemToGenesisTet10;
3422  break;
3423  }
3424  case (ELEMENT_HEX27):
3425  {
3426  mymap = (int *) mfemToGenesisHex27;
3427  break;
3428  }
3429  case (ELEMENT_TRI3):
3430  case (ELEMENT_QUAD4):
3431  case (ELEMENT_TET4):
3432  case (ELEMENT_HEX8):
3433  {
3434  MFEM_ABORT("Something went wrong. Linear elements detected when order is 2.");
3435  break;
3436  }
3437  }
3438 
3439  FinalizeTopology();
3440 
3441  // Define quadratic FE space
3442  FiniteElementCollection *fec = new H1_FECollection(2,3);
3443  FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
3444  Ordering::byVDIM);
3445  Nodes = new GridFunction(fes);
3446  Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
3447  own_nodes = 1;
3448 
3449  // int nTotDofs = fes->GetNDofs();
3450  // int nTotVDofs = fes->GetVSize();
3451  // mfem::out << endl << "nTotDofs = " << nTotDofs << " nTotVDofs "
3452  // << nTotVDofs << endl << endl;
3453 
3454  for (int i = 0; i < NumOfElements; i++)
3455  {
3456  Array<int> dofs;
3457 
3458  fes->GetElementDofs(i, dofs);
3459  Array<int> vdofs;
3460  vdofs.SetSize(dofs.Size());
3461  for (int l = 0; l < dofs.Size(); l++) { vdofs[l] = dofs[l]; }
3462  fes->DofsToVDofs(vdofs);
3463  int iblk = 0;
3464  int loc_ind;
3465  while (iblk < (int) num_el_blk && i >= start_of_block[iblk+1]) { iblk++; }
3466  loc_ind = i - start_of_block[iblk];
3467  for (int j = 0; j < dofs.Size(); j++)
3468  {
3469  int point_id = elem_blk[iblk][loc_ind*num_node_per_el + mymap[j] - 1] - 1;
3470  (*Nodes)(vdofs[j]) = coordx[point_id];
3471  (*Nodes)(vdofs[j]+1) = coordy[point_id];
3472  if (Dim == 3)
3473  {
3474  (*Nodes)(vdofs[j]+2) = coordz[point_id];
3475  }
3476  }
3477  }
3478  }
3479 
3480  // clean up all netcdf stuff
3481 
3482  nc_close(ncid);
3483 
3484  for (int i = 0; i < (int) num_side_sets; i++)
3485  {
3486  delete [] elem_ss[i];
3487  delete [] side_ss[i];
3488  }
3489 
3490  delete [] elem_ss;
3491  delete [] side_ss;
3492  delete [] num_el_in_blk;
3493  delete [] num_side_in_ss;
3494  delete [] coordx;
3495  delete [] coordy;
3496  delete [] coordz;
3497 
3498  for (int i = 0; i < (int) num_el_blk; i++)
3499  {
3500  delete [] elem_blk[i];
3501  }
3502 
3503  delete [] elem_blk;
3504  delete [] start_of_block;
3505 
3506  for (int i = 0; i < (int) num_side_sets; i++)
3507  {
3508  delete [] ss_node_id[i];
3509  }
3510  delete [] ss_node_id;
3511  delete [] ebprop;
3512  delete [] ssprop;
3513 
3514 }
3515 #endif // #ifdef MFEM_USE_NETCDF
3516 
3517 } // namespace mfem
Abstract class for all finite elements.
Definition: fe_base.hpp:232
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:638
Vector bb_max
Definition: ex9.cpp:67
Mesh * Make3D(int nsteps, double rstep, double aspect, int order, bool sfc)
Definition: polar-nc.cpp:370
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
Definition: array.cpp:53
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void GmshHOSegmentMapping(int order, int *map)
Generate Gmsh vertex mapping for a Segment.
Definition: gmsh.cpp:378
T * GetData()
Returns the data.
Definition: array.hpp:115
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:122
Data type for vertex.
Definition: vertex.hpp:22
STL namespace.
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)
Check if the stream starts with comment_char. If so skip it.
Definition: text.hpp:31
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
void DecodeBase64(const char *src, size_t len, std::vector< char > &buf)
Decode len base-64 encoded characters in the buffer src, and store the resulting decoded data in buf...
Definition: binaryio.cpp:74
Mesh * Make2D(int nsteps, double rstep, double phi, double aspect, int order, bool sfc)
Definition: polar-nc.cpp:49
Class for standard nodal finite elements.
Definition: fe_base.hpp:686
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:57
Data type hexahedron element.
Definition: hexahedron.hpp:22
Data type Pyramid element.
Definition: pyramid.hpp:22
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
Data type triangle element.
Definition: triangle.hpp:23
size_t NumBase64Chars(size_t nbytes)
Return the number of characters needed to encode nbytes in base-64.
Definition: binaryio.cpp:103
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
void GmshHOQuadrilateralMapping(int order, int *map)
Generate Gmsh vertex mapping for a Quadrilateral.
Definition: gmsh.cpp:403
const Array< int > & GetLexicographicOrdering() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:761
Data type tetrahedron element.
Definition: tetrahedron.hpp:22
Vector bb_min
Definition: ex9.cpp:67
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
void GmshHOTriangleMapping(int order, int *map)
Generate Gmsh vertex mapping for a Triangle.
Definition: gmsh.cpp:388
void GmshHOTetrahedronMapping(int order, int *map)
Generate Gmsh vertex mapping for a Tetrahedron.
Definition: gmsh.cpp:417
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
double a
Definition: lissajous.cpp:41
bool StringCompare(const char *s1, const char *s2)
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
Definition: vtk.cpp:602
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:723
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:196
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void filter_dos(std::string &line)
Check for, and remove, a trailing &#39;\r&#39; from and std::string.
Definition: text.hpp:45
Vector data type.
Definition: vector.hpp:60
Data type point element.
Definition: point.hpp:22
Vector coefficient defined by a vector GridFunction.
void GmshHOHexahedronMapping(int order, int *map)
Generate Gmsh vertex mapping for a Hexahedron.
Definition: gmsh.cpp:436
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Create the VTK element connectivity array for a given element geometry and refinement level...
Definition: vtk.cpp:497
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
void GmshHOWedgeMapping(int order, int *map)
Generate Gmsh vertex mapping for a Wedge.
Definition: gmsh.cpp:453
Abstract data type element.
Definition: element.hpp:28
Data type line segment element.
Definition: segment.hpp:22
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:215
void GmshHOPyramidMapping(int order, int *map)
Generate Gmsh vertex mapping for a Pyramid.
Definition: gmsh.cpp:470
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320