MFEM  v4.5.2
Finite element discretization library
mesh.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 // Implementation of data type mesh
13 
14 #include "mesh_headers.hpp"
15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 #include "../general/binaryio.hpp"
18 #include "../general/text.hpp"
19 #include "../general/device.hpp"
20 #include "../general/tic_toc.hpp"
21 #include "../general/gecko.hpp"
22 #include "../fem/quadinterpolator.hpp"
23 
24 #include <iostream>
25 #include <sstream>
26 #include <fstream>
27 #include <limits>
28 #include <cmath>
29 #include <cstring>
30 #include <ctime>
31 #include <functional>
32 #include <map>
33 #include <set>
34 
35 // Include the METIS header, if using version 5. If using METIS 4, the needed
36 // declarations are inlined below, i.e. no header is needed.
37 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
38 #include "metis.h"
39 #endif
40 
41 // METIS 4 prototypes
42 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
43 typedef int idx_t;
44 typedef int idxtype;
45 extern "C" {
47  int*, int*, int*, int*, int*, idxtype*);
49  int*, int*, int*, int*, int*, idxtype*);
51  int*, int*, int*, int*, int*, idxtype*);
52 }
53 #endif
54 
55 using namespace std;
56 
57 namespace mfem
58 {
59 
60 void Mesh::GetElementJacobian(int i, DenseMatrix &J)
61 {
62  Geometry::Type geom = GetElementBaseGeometry(i);
63  ElementTransformation *eltransf = GetElementTransformation(i);
64  eltransf->SetIntPoint(&Geometries.GetCenter(geom));
65  Geometries.JacToPerfJac(geom, eltransf->Jacobian(), J);
66 }
67 
68 void Mesh::GetElementCenter(int i, Vector &center)
69 {
70  center.SetSize(spaceDim);
71  int geom = GetElementBaseGeometry(i);
72  ElementTransformation *eltransf = GetElementTransformation(i);
73  eltransf->Transform(Geometries.GetCenter(geom), center);
74 }
75 
76 double Mesh::GetElementSize(ElementTransformation *T, int type)
77 {
78  DenseMatrix J(spaceDim, Dim);
79 
80  Geometry::Type geom = T->GetGeometryType();
81  T->SetIntPoint(&Geometries.GetCenter(geom));
82  Geometries.JacToPerfJac(geom, T->Jacobian(), J);
83 
84  if (type == 0)
85  {
86  return pow(fabs(J.Weight()), 1./Dim);
87  }
88  else if (type == 1)
89  {
90  return J.CalcSingularvalue(Dim-1); // h_min
91  }
92  else
93  {
94  return J.CalcSingularvalue(0); // h_max
95  }
96 }
97 
98 double Mesh::GetElementSize(int i, int type)
99 {
100  return GetElementSize(GetElementTransformation(i), type);
101 }
102 
103 double Mesh::GetElementSize(int i, const Vector &dir)
104 {
105  DenseMatrix J(spaceDim, Dim);
106  Vector d_hat(Dim);
107  GetElementJacobian(i, J);
108  J.MultTranspose(dir, d_hat);
109  return sqrt((d_hat * d_hat) / (dir * dir));
110 }
111 
112 double Mesh::GetElementVolume(int i)
113 {
114  ElementTransformation *et = GetElementTransformation(i);
115  const IntegrationRule &ir = IntRules.Get(GetElementBaseGeometry(i),
116  et->OrderJ());
117  double volume = 0.0;
118  for (int j = 0; j < ir.GetNPoints(); j++)
119  {
120  const IntegrationPoint &ip = ir.IntPoint(j);
121  et->SetIntPoint(&ip);
122  volume += ip.weight * et->Weight();
123  }
124 
125  return volume;
126 }
127 
128 // Similar to VisualizationSceneSolution3d::FindNewBox in GLVis
129 void Mesh::GetBoundingBox(Vector &min, Vector &max, int ref)
130 {
131  min.SetSize(spaceDim);
132  max.SetSize(spaceDim);
133 
134  for (int d = 0; d < spaceDim; d++)
135  {
136  min(d) = infinity();
137  max(d) = -infinity();
138  }
139 
140  if (Nodes == NULL)
141  {
142  double *coord;
143  for (int i = 0; i < NumOfVertices; i++)
144  {
145  coord = GetVertex(i);
146  for (int d = 0; d < spaceDim; d++)
147  {
148  if (coord[d] < min(d)) { min(d) = coord[d]; }
149  if (coord[d] > max(d)) { max(d) = coord[d]; }
150  }
151  }
152  }
153  else
154  {
155  const bool use_boundary = false; // make this a parameter?
156  int ne = use_boundary ? GetNBE() : GetNE();
157  int fn, fo;
158  DenseMatrix pointmat;
159  RefinedGeometry *RefG;
160  IntegrationRule eir;
163 
164  for (int i = 0; i < ne; i++)
165  {
166  if (use_boundary)
167  {
168  GetBdrElementFace(i, &fn, &fo);
169  RefG = GlobGeometryRefiner.Refine(GetFaceGeometry(fn), ref);
170  Tr = GetFaceElementTransformations(fn, 5);
171  eir.SetSize(RefG->RefPts.GetNPoints());
172  Tr->Loc1.Transform(RefG->RefPts, eir);
173  Tr->Elem1->Transform(eir, pointmat);
174  }
175  else
176  {
177  T = GetElementTransformation(i);
178  RefG = GlobGeometryRefiner.Refine(GetElementBaseGeometry(i), ref);
179  T->Transform(RefG->RefPts, pointmat);
180  }
181  for (int j = 0; j < pointmat.Width(); j++)
182  {
183  for (int d = 0; d < pointmat.Height(); d++)
184  {
185  if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
186  if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
187  }
188  }
189  }
190  }
191 }
192 
193 void Mesh::GetCharacteristics(double &h_min, double &h_max,
194  double &kappa_min, double &kappa_max,
195  Vector *Vh, Vector *Vk)
196 {
197  int i, dim, sdim;
198  DenseMatrix J;
199  double h, kappa;
200 
201  dim = Dimension();
202  sdim = SpaceDimension();
203 
204  if (Vh) { Vh->SetSize(NumOfElements); }
205  if (Vk) { Vk->SetSize(NumOfElements); }
206 
207  h_min = kappa_min = infinity();
208  h_max = kappa_max = -h_min;
209  if (dim == 0) { if (Vh) { *Vh = 1.0; } if (Vk) {*Vk = 1.0; } return; }
210  J.SetSize(sdim, dim);
211  for (i = 0; i < NumOfElements; i++)
212  {
213  GetElementJacobian(i, J);
214  h = pow(fabs(J.Weight()), 1.0/double(dim));
215  kappa = (dim == sdim) ?
216  J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1) : -1.0;
217  if (Vh) { (*Vh)(i) = h; }
218  if (Vk) { (*Vk)(i) = kappa; }
219 
220  if (h < h_min) { h_min = h; }
221  if (h > h_max) { h_max = h; }
222  if (kappa < kappa_min) { kappa_min = kappa; }
223  if (kappa > kappa_max) { kappa_max = kappa; }
224  }
225 }
226 
227 // static method
228 void Mesh::PrintElementsByGeometry(int dim,
229  const Array<int> &num_elems_by_geom,
230  std::ostream &os)
231 {
232  for (int g = Geometry::DimStart[dim], first = 1;
233  g < Geometry::DimStart[dim+1]; g++)
234  {
235  if (!num_elems_by_geom[g]) { continue; }
236  if (!first) { os << " + "; }
237  else { first = 0; }
238  os << num_elems_by_geom[g] << ' ' << Geometry::Name[g] << "(s)";
239  }
240 }
241 
242 void Mesh::PrintCharacteristics(Vector *Vh, Vector *Vk, std::ostream &os)
243 {
244  double h_min, h_max, kappa_min, kappa_max;
245 
246  os << "Mesh Characteristics:";
247 
248  this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
249 
250  Array<int> num_elems_by_geom(Geometry::NumGeom);
251  num_elems_by_geom = 0;
252  for (int i = 0; i < GetNE(); i++)
253  {
254  num_elems_by_geom[GetElementBaseGeometry(i)]++;
255  }
256 
257  os << '\n'
258  << "Dimension : " << Dimension() << '\n'
259  << "Space dimension : " << SpaceDimension();
260  if (Dim == 0)
261  {
262  os << '\n'
263  << "Number of vertices : " << GetNV() << '\n'
264  << "Number of elements : " << GetNE() << '\n'
265  << "Number of bdr elem : " << GetNBE() << '\n';
266  }
267  else if (Dim == 1)
268  {
269  os << '\n'
270  << "Number of vertices : " << GetNV() << '\n'
271  << "Number of elements : " << GetNE() << '\n'
272  << "Number of bdr elem : " << GetNBE() << '\n'
273  << "h_min : " << h_min << '\n'
274  << "h_max : " << h_max << '\n';
275  }
276  else if (Dim == 2)
277  {
278  os << '\n'
279  << "Number of vertices : " << GetNV() << '\n'
280  << "Number of edges : " << GetNEdges() << '\n'
281  << "Number of elements : " << GetNE() << " -- ";
282  PrintElementsByGeometry(2, num_elems_by_geom, os);
283  os << '\n'
284  << "Number of bdr elem : " << GetNBE() << '\n'
285  << "Euler Number : " << EulerNumber2D() << '\n'
286  << "h_min : " << h_min << '\n'
287  << "h_max : " << h_max << '\n'
288  << "kappa_min : " << kappa_min << '\n'
289  << "kappa_max : " << kappa_max << '\n';
290  }
291  else
292  {
293  Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
294  num_bdr_elems_by_geom = 0;
295  for (int i = 0; i < GetNBE(); i++)
296  {
297  num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
298  }
299  Array<int> num_faces_by_geom(Geometry::NumGeom);
300  num_faces_by_geom = 0;
301  for (int i = 0; i < GetNFaces(); i++)
302  {
303  num_faces_by_geom[GetFaceGeometry(i)]++;
304  }
305 
306  os << '\n'
307  << "Number of vertices : " << GetNV() << '\n'
308  << "Number of edges : " << GetNEdges() << '\n'
309  << "Number of faces : " << GetNFaces() << " -- ";
310  PrintElementsByGeometry(Dim-1, num_faces_by_geom, os);
311  os << '\n'
312  << "Number of elements : " << GetNE() << " -- ";
313  PrintElementsByGeometry(Dim, num_elems_by_geom, os);
314  os << '\n'
315  << "Number of bdr elem : " << GetNBE() << " -- ";
316  PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, os);
317  os << '\n'
318  << "Euler Number : " << EulerNumber() << '\n'
319  << "h_min : " << h_min << '\n'
320  << "h_max : " << h_max << '\n'
321  << "kappa_min : " << kappa_min << '\n'
322  << "kappa_max : " << kappa_max << '\n';
323  }
324  os << '\n' << std::flush;
325 }
326 
327 FiniteElement *Mesh::GetTransformationFEforElementType(Element::Type ElemType)
328 {
329  switch (ElemType)
330  {
331  case Element::POINT : return &PointFE;
332  case Element::SEGMENT : return &SegmentFE;
333  case Element::TRIANGLE : return &TriangleFE;
334  case Element::QUADRILATERAL : return &QuadrilateralFE;
335  case Element::TETRAHEDRON : return &TetrahedronFE;
336  case Element::HEXAHEDRON : return &HexahedronFE;
337  case Element::WEDGE : return &WedgeFE;
338  case Element::PYRAMID : return &PyramidFE;
339  default:
340  MFEM_ABORT("Unknown element type \"" << ElemType << "\"");
341  break;
342  }
343  MFEM_ABORT("Unknown element type");
344  return NULL;
345 }
346 
347 
348 void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
349 {
350  ElTr->Attribute = GetAttribute(i);
351  ElTr->ElementNo = i;
352  ElTr->ElementType = ElementTransformation::ELEMENT;
353  ElTr->mesh = this;
354  ElTr->Reset();
355  if (Nodes == NULL)
356  {
357  GetPointMatrix(i, ElTr->GetPointMat());
358  ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
359  }
360  else
361  {
362  DenseMatrix &pm = ElTr->GetPointMat();
363  Array<int> vdofs;
364  Nodes->FESpace()->GetElementVDofs(i, vdofs);
365  Nodes->HostRead();
366  const GridFunction &nodes = *Nodes;
367  int n = vdofs.Size()/spaceDim;
368  pm.SetSize(spaceDim, n);
369  for (int k = 0; k < spaceDim; k++)
370  {
371  for (int j = 0; j < n; j++)
372  {
373  pm(k,j) = nodes(vdofs[n*k+j]);
374  }
375  }
376  ElTr->SetFE(Nodes->FESpace()->GetFE(i));
377  }
378 }
379 
380 void Mesh::GetElementTransformation(int i, const Vector &nodes,
382 {
383  ElTr->Attribute = GetAttribute(i);
384  ElTr->ElementNo = i;
385  ElTr->ElementType = ElementTransformation::ELEMENT;
386  ElTr->mesh = this;
387  DenseMatrix &pm = ElTr->GetPointMat();
388  ElTr->Reset();
389  nodes.HostRead();
390  if (Nodes == NULL)
391  {
392  MFEM_ASSERT(nodes.Size() == spaceDim*GetNV(), "");
393  int nv = elements[i]->GetNVertices();
394  const int *v = elements[i]->GetVertices();
395  int n = vertices.Size();
396  pm.SetSize(spaceDim, nv);
397  for (int k = 0; k < spaceDim; k++)
398  {
399  for (int j = 0; j < nv; j++)
400  {
401  pm(k, j) = nodes(k*n+v[j]);
402  }
403  }
404  ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
405  }
406  else
407  {
408  MFEM_ASSERT(nodes.Size() == Nodes->Size(), "");
409  Array<int> vdofs;
410  Nodes->FESpace()->GetElementVDofs(i, vdofs);
411  int n = vdofs.Size()/spaceDim;
412  pm.SetSize(spaceDim, n);
413  for (int k = 0; k < spaceDim; k++)
414  {
415  for (int j = 0; j < n; j++)
416  {
417  pm(k,j) = nodes(vdofs[n*k+j]);
418  }
419  }
420  ElTr->SetFE(Nodes->FESpace()->GetFE(i));
421  }
422 }
423 
424 ElementTransformation *Mesh::GetElementTransformation(int i)
425 {
426  GetElementTransformation(i, &Transformation);
427 
428  return &Transformation;
429 }
430 
431 ElementTransformation *Mesh::GetBdrElementTransformation(int i)
432 {
433  GetBdrElementTransformation(i, &BdrTransformation);
434  return &BdrTransformation;
435 }
436 
437 void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
438 {
439  ElTr->Attribute = GetBdrAttribute(i);
440  ElTr->ElementNo = i; // boundary element number
441  ElTr->ElementType = ElementTransformation::BDR_ELEMENT;
442  ElTr->mesh = this;
443  DenseMatrix &pm = ElTr->GetPointMat();
444  ElTr->Reset();
445  if (Nodes == NULL)
446  {
447  GetBdrPointMatrix(i, pm);
448  ElTr->SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
449  }
450  else
451  {
452  const FiniteElement *bdr_el = Nodes->FESpace()->GetBE(i);
453  Nodes->HostRead();
454  const GridFunction &nodes = *Nodes;
455  if (bdr_el)
456  {
457  Array<int> vdofs;
458  Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
459  int n = vdofs.Size()/spaceDim;
460  pm.SetSize(spaceDim, n);
461  for (int k = 0; k < spaceDim; k++)
462  {
463  for (int j = 0; j < n; j++)
464  {
465  pm(k,j) = nodes(vdofs[n*k+j]);
466  }
467  }
468  ElTr->SetFE(bdr_el);
469  }
470  else // L2 Nodes (e.g., periodic mesh)
471  {
472  int elem_id, face_info;
473  GetBdrElementAdjacentElement2(i, elem_id, face_info);
474 
475  GetLocalFaceTransformation(GetBdrElementType(i),
476  GetElementType(elem_id),
477  FaceElemTr.Loc1.Transf, face_info);
478  // NOTE: FaceElemTr.Loc1 is overwritten here -- used as a temporary
479 
480  Geometry::Type face_geom = GetBdrElementBaseGeometry(i);
481  const FiniteElement *face_el =
482  Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
483  MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
484  "Mesh requires nodal Finite Element.");
485  IntegrationRule eir(face_el->GetDof());
486  FaceElemTr.Loc1.Transf.ElementNo = elem_id;
487  FaceElemTr.Loc1.Transf.mesh = this;
488  FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
489  FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
490  Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
491 
492  ElTr->SetFE(face_el);
493  }
494  }
495 }
496 
497 void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
498 {
499  FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
500  FTr->ElementNo = FaceNo;
501  FTr->ElementType = ElementTransformation::FACE;
502  FTr->mesh = this;
503  DenseMatrix &pm = FTr->GetPointMat();
504  FTr->Reset();
505  if (Nodes == NULL)
506  {
507  const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
508  const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
509  pm.SetSize(spaceDim, nv);
510  for (int i = 0; i < spaceDim; i++)
511  {
512  for (int j = 0; j < nv; j++)
513  {
514  pm(i, j) = vertices[v[j]](i);
515  }
516  }
517  FTr->SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
518  }
519  else // curved mesh
520  {
521  const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
522  Nodes->HostRead();
523  const GridFunction &nodes = *Nodes;
524  if (face_el)
525  {
526  Array<int> vdofs;
527  Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
528  int n = vdofs.Size()/spaceDim;
529  pm.SetSize(spaceDim, n);
530  for (int i = 0; i < spaceDim; i++)
531  {
532  for (int j = 0; j < n; j++)
533  {
534  pm(i, j) = nodes(vdofs[n*i+j]);
535  }
536  }
537  FTr->SetFE(face_el);
538  }
539  else // L2 Nodes (e.g., periodic mesh), go through the volume of Elem1
540  {
541  FaceInfo &face_info = faces_info[FaceNo];
542 
543  Geometry::Type face_geom = GetFaceGeometry(FaceNo);
544  Element::Type face_type = GetFaceElementType(FaceNo);
545 
546  GetLocalFaceTransformation(face_type,
547  GetElementType(face_info.Elem1No),
548  FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
549  // NOTE: FaceElemTr.Loc1 is overwritten here -- used as a temporary
550 
551  face_el = Nodes->FESpace()->GetTraceElement(face_info.Elem1No,
552  face_geom);
553  MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
554  "Mesh requires nodal Finite Element.");
555 
556  IntegrationRule eir(face_el->GetDof());
557  FaceElemTr.Loc1.Transf.ElementNo = face_info.Elem1No;
558  FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
559  FaceElemTr.Loc1.Transf.mesh = this;
560  FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
561  Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
562 
563  FTr->SetFE(face_el);
564  }
565  }
566 }
567 
568 ElementTransformation *Mesh::GetFaceTransformation(int FaceNo)
569 {
570  GetFaceTransformation(FaceNo, &FaceTransformation);
571  return &FaceTransformation;
572 }
573 
574 void Mesh::GetEdgeTransformation(int EdgeNo, IsoparametricTransformation *EdTr)
575 {
576  if (Dim == 2)
577  {
578  GetFaceTransformation(EdgeNo, EdTr);
579  return;
580  }
581  if (Dim == 1)
582  {
583  mfem_error("Mesh::GetEdgeTransformation not defined in 1D \n");
584  }
585 
586  EdTr->Attribute = 1;
587  EdTr->ElementNo = EdgeNo;
588  EdTr->ElementType = ElementTransformation::EDGE;
589  EdTr->mesh = this;
590  DenseMatrix &pm = EdTr->GetPointMat();
591  EdTr->Reset();
592  if (Nodes == NULL)
593  {
594  Array<int> v;
595  GetEdgeVertices(EdgeNo, v);
596  const int nv = 2;
597  pm.SetSize(spaceDim, nv);
598  for (int i = 0; i < spaceDim; i++)
599  {
600  for (int j = 0; j < nv; j++)
601  {
602  pm(i, j) = vertices[v[j]](i);
603  }
604  }
605  EdTr->SetFE(GetTransformationFEforElementType(Element::SEGMENT));
606  }
607  else
608  {
609  const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
610  Nodes->HostRead();
611  const GridFunction &nodes = *Nodes;
612  if (edge_el)
613  {
614  Array<int> vdofs;
615  Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
616  int n = vdofs.Size()/spaceDim;
617  pm.SetSize(spaceDim, n);
618  for (int i = 0; i < spaceDim; i++)
619  {
620  for (int j = 0; j < n; j++)
621  {
622  pm(i, j) = nodes(vdofs[n*i+j]);
623  }
624  }
625  EdTr->SetFE(edge_el);
626  }
627  else
628  {
629  MFEM_ABORT("Not implemented.");
630  }
631  }
632 }
633 
634 ElementTransformation *Mesh::GetEdgeTransformation(int EdgeNo)
635 {
636  GetEdgeTransformation(EdgeNo, &EdgeTransformation);
637  return &EdgeTransformation;
638 }
639 
640 
641 void Mesh::GetLocalPtToSegTransformation(
642  IsoparametricTransformation &Transf, int i)
643 {
644  const IntegrationRule *SegVert;
645  DenseMatrix &locpm = Transf.GetPointMat();
646  Transf.Reset();
647 
648  Transf.SetFE(&PointFE);
649  SegVert = Geometries.GetVertices(Geometry::SEGMENT);
650  locpm.SetSize(1, 1);
651  locpm(0, 0) = SegVert->IntPoint(i/64).x;
652  // (i/64) is the local face no. in the segment
653  // (i%64) is the orientation of the point (not used)
654 }
655 
656 void Mesh::GetLocalSegToTriTransformation(
657  IsoparametricTransformation &Transf, int i)
658 {
659  const int *tv, *so;
660  const IntegrationRule *TriVert;
661  DenseMatrix &locpm = Transf.GetPointMat();
662  Transf.Reset();
663 
664  Transf.SetFE(&SegmentFE);
665  tv = tri_t::Edges[i/64]; // (i/64) is the local face no. in the triangle
666  so = seg_t::Orient[i%64]; // (i%64) is the orientation of the segment
667  TriVert = Geometries.GetVertices(Geometry::TRIANGLE);
668  locpm.SetSize(2, 2);
669  for (int j = 0; j < 2; j++)
670  {
671  locpm(0, so[j]) = TriVert->IntPoint(tv[j]).x;
672  locpm(1, so[j]) = TriVert->IntPoint(tv[j]).y;
673  }
674 }
675 
676 void Mesh::GetLocalSegToQuadTransformation(
677  IsoparametricTransformation &Transf, int i)
678 {
679  const int *qv, *so;
680  const IntegrationRule *QuadVert;
681  DenseMatrix &locpm = Transf.GetPointMat();
682  Transf.Reset();
683 
684  Transf.SetFE(&SegmentFE);
685  qv = quad_t::Edges[i/64]; // (i/64) is the local face no. in the quad
686  so = seg_t::Orient[i%64]; // (i%64) is the orientation of the segment
687  QuadVert = Geometries.GetVertices(Geometry::SQUARE);
688  locpm.SetSize(2, 2);
689  for (int j = 0; j < 2; j++)
690  {
691  locpm(0, so[j]) = QuadVert->IntPoint(qv[j]).x;
692  locpm(1, so[j]) = QuadVert->IntPoint(qv[j]).y;
693  }
694 }
695 
696 void Mesh::GetLocalTriToTetTransformation(
697  IsoparametricTransformation &Transf, int i)
698 {
699  DenseMatrix &locpm = Transf.GetPointMat();
700  Transf.Reset();
701 
702  Transf.SetFE(&TriangleFE);
703  // (i/64) is the local face no. in the tet
704  const int *tv = tet_t::FaceVert[i/64];
705  // (i%64) is the orientation of the tetrahedron face
706  // w.r.t. the face element
707  const int *to = tri_t::Orient[i%64];
708  const IntegrationRule *TetVert =
709  Geometries.GetVertices(Geometry::TETRAHEDRON);
710  locpm.SetSize(3, 3);
711  for (int j = 0; j < 3; j++)
712  {
713  const IntegrationPoint &vert = TetVert->IntPoint(tv[to[j]]);
714  locpm(0, j) = vert.x;
715  locpm(1, j) = vert.y;
716  locpm(2, j) = vert.z;
717  }
718 }
719 
720 void Mesh::GetLocalTriToWdgTransformation(
721  IsoparametricTransformation &Transf, int i)
722 {
723  DenseMatrix &locpm = Transf.GetPointMat();
724  Transf.Reset();
725 
726  Transf.SetFE(&TriangleFE);
727  // (i/64) is the local face no. in the pri
728  MFEM_VERIFY(i < 128, "Local face index " << i/64
729  << " is not a triangular face of a wedge.");
730  const int *pv = pri_t::FaceVert[i/64];
731  // (i%64) is the orientation of the wedge face
732  // w.r.t. the face element
733  const int *to = tri_t::Orient[i%64];
734  const IntegrationRule *PriVert =
735  Geometries.GetVertices(Geometry::PRISM);
736  locpm.SetSize(3, 3);
737  for (int j = 0; j < 3; j++)
738  {
739  const IntegrationPoint &vert = PriVert->IntPoint(pv[to[j]]);
740  locpm(0, j) = vert.x;
741  locpm(1, j) = vert.y;
742  locpm(2, j) = vert.z;
743  }
744 }
745 
746 void Mesh::GetLocalTriToPyrTransformation(
747  IsoparametricTransformation &Transf, int i)
748 {
749  DenseMatrix &locpm = Transf.GetPointMat();
750 
751  Transf.SetFE(&TriangleFE);
752  // (i/64) is the local face no. in the pyr
753  MFEM_VERIFY(i >= 64, "Local face index " << i/64
754  << " is not a triangular face of a pyramid.");
755  const int *pv = pyr_t::FaceVert[i/64];
756  // (i%64) is the orientation of the pyramid face
757  // w.r.t. the face element
758  const int *to = tri_t::Orient[i%64];
759  const IntegrationRule *PyrVert =
760  Geometries.GetVertices(Geometry::PYRAMID);
761  locpm.SetSize(3, 3);
762  for (int j = 0; j < 3; j++)
763  {
764  const IntegrationPoint &vert = PyrVert->IntPoint(pv[to[j]]);
765  locpm(0, j) = vert.x;
766  locpm(1, j) = vert.y;
767  locpm(2, j) = vert.z;
768  }
769 }
770 
771 void Mesh::GetLocalQuadToHexTransformation(
772  IsoparametricTransformation &Transf, int i)
773 {
774  DenseMatrix &locpm = Transf.GetPointMat();
775  Transf.Reset();
776 
777  Transf.SetFE(&QuadrilateralFE);
778  // (i/64) is the local face no. in the hex
779  const int *hv = hex_t::FaceVert[i/64];
780  // (i%64) is the orientation of the quad
781  const int *qo = quad_t::Orient[i%64];
782  const IntegrationRule *HexVert = Geometries.GetVertices(Geometry::CUBE);
783  locpm.SetSize(3, 4);
784  for (int j = 0; j < 4; j++)
785  {
786  const IntegrationPoint &vert = HexVert->IntPoint(hv[qo[j]]);
787  locpm(0, j) = vert.x;
788  locpm(1, j) = vert.y;
789  locpm(2, j) = vert.z;
790  }
791 }
792 
793 void Mesh::GetLocalQuadToWdgTransformation(
794  IsoparametricTransformation &Transf, int i)
795 {
796  DenseMatrix &locpm = Transf.GetPointMat();
797  Transf.Reset();
798 
799  Transf.SetFE(&QuadrilateralFE);
800  // (i/64) is the local face no. in the pri
801  MFEM_VERIFY(i >= 128, "Local face index " << i/64
802  << " is not a quadrilateral face of a wedge.");
803  const int *pv = pri_t::FaceVert[i/64];
804  // (i%64) is the orientation of the quad
805  const int *qo = quad_t::Orient[i%64];
806  const IntegrationRule *PriVert = Geometries.GetVertices(Geometry::PRISM);
807  locpm.SetSize(3, 4);
808  for (int j = 0; j < 4; j++)
809  {
810  const IntegrationPoint &vert = PriVert->IntPoint(pv[qo[j]]);
811  locpm(0, j) = vert.x;
812  locpm(1, j) = vert.y;
813  locpm(2, j) = vert.z;
814  }
815 }
816 
817 void Mesh::GetLocalQuadToPyrTransformation(
818  IsoparametricTransformation &Transf, int i)
819 {
820  DenseMatrix &locpm = Transf.GetPointMat();
821 
822  Transf.SetFE(&QuadrilateralFE);
823  // (i/64) is the local face no. in the pyr
824  MFEM_VERIFY(i < 64, "Local face index " << i/64
825  << " is not a quadrilateral face of a pyramid.");
826  const int *pv = pyr_t::FaceVert[i/64];
827  // (i%64) is the orientation of the quad
828  const int *qo = quad_t::Orient[i%64];
829  const IntegrationRule *PyrVert = Geometries.GetVertices(Geometry::PYRAMID);
830  locpm.SetSize(3, 4);
831  for (int j = 0; j < 4; j++)
832  {
833  const IntegrationPoint &vert = PyrVert->IntPoint(pv[qo[j]]);
834  locpm(0, j) = vert.x;
835  locpm(1, j) = vert.y;
836  locpm(2, j) = vert.z;
837  }
838 }
839 
840 const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
841  const int flags,
842  MemoryType d_mt)
843 {
844  for (int i = 0; i < geom_factors.Size(); i++)
845  {
846  GeometricFactors *gf = geom_factors[i];
847  if (gf->IntRule == &ir && (gf->computed_factors & flags) == flags)
848  {
849  return gf;
850  }
851  }
852 
853  this->EnsureNodes();
854 
855  GeometricFactors *gf = new GeometricFactors(this, ir, flags, d_mt);
856  geom_factors.Append(gf);
857  return gf;
858 }
859 
860 const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
861  const IntegrationRule& ir,
862  const int flags, FaceType type, MemoryType d_mt)
863 {
864  for (int i = 0; i < face_geom_factors.Size(); i++)
865  {
866  FaceGeometricFactors *gf = face_geom_factors[i];
867  if (gf->IntRule == &ir && (gf->computed_factors & flags) == flags &&
868  gf->type==type)
869  {
870  return gf;
871  }
872  }
873 
874  this->EnsureNodes();
875 
876  FaceGeometricFactors *gf = new FaceGeometricFactors(this, ir, flags, type,
877  d_mt);
878  face_geom_factors.Append(gf);
879  return gf;
880 }
881 
882 void Mesh::DeleteGeometricFactors()
883 {
884  for (int i = 0; i < geom_factors.Size(); i++)
885  {
886  delete geom_factors[i];
887  }
888  geom_factors.SetSize(0);
889  for (int i = 0; i < face_geom_factors.Size(); i++)
890  {
891  delete face_geom_factors[i];
892  }
893  face_geom_factors.SetSize(0);
894 }
895 
896 void Mesh::GetLocalFaceTransformation(
897  int face_type, int elem_type, IsoparametricTransformation &Transf, int info)
898 {
899  switch (face_type)
900  {
901  case Element::POINT:
902  GetLocalPtToSegTransformation(Transf, info);
903  break;
904 
905  case Element::SEGMENT:
906  if (elem_type == Element::TRIANGLE)
907  {
908  GetLocalSegToTriTransformation(Transf, info);
909  }
910  else
911  {
912  MFEM_ASSERT(elem_type == Element::QUADRILATERAL, "");
913  GetLocalSegToQuadTransformation(Transf, info);
914  }
915  break;
916 
917  case Element::TRIANGLE:
918  if (elem_type == Element::TETRAHEDRON)
919  {
920  GetLocalTriToTetTransformation(Transf, info);
921  }
922  else if (elem_type == Element::WEDGE)
923  {
924  GetLocalTriToWdgTransformation(Transf, info);
925  }
926  else if (elem_type == Element::PYRAMID)
927  {
928  GetLocalTriToPyrTransformation(Transf, info);
929  }
930  else
931  {
932  MFEM_ABORT("Mesh::GetLocalFaceTransformation not defined for "
933  "face type " << face_type
934  << " and element type " << elem_type << "\n");
935  }
936  break;
937 
938  case Element::QUADRILATERAL:
939  if (elem_type == Element::HEXAHEDRON)
940  {
941  GetLocalQuadToHexTransformation(Transf, info);
942  }
943  else if (elem_type == Element::WEDGE)
944  {
945  GetLocalQuadToWdgTransformation(Transf, info);
946  }
947  else if (elem_type == Element::PYRAMID)
948  {
949  GetLocalQuadToPyrTransformation(Transf, info);
950  }
951  else
952  {
953  MFEM_ABORT("Mesh::GetLocalFaceTransformation not defined for "
954  "face type " << face_type
955  << " and element type " << elem_type << "\n");
956  }
957  break;
958  }
959 }
960 
961 FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
962  int mask)
963 {
964  FaceInfo &face_info = faces_info[FaceNo];
965 
966  int cmask = 0;
967  FaceElemTr.SetConfigurationMask(cmask);
968  FaceElemTr.Elem1 = NULL;
969  FaceElemTr.Elem2 = NULL;
970 
971  // setup the transformation for the first element
972  FaceElemTr.Elem1No = face_info.Elem1No;
973  if (mask & FaceElementTransformations::HAVE_ELEM1)
974  {
975  GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
976  FaceElemTr.Elem1 = &Transformation;
977  cmask |= 1;
978  }
979 
980  // setup the transformation for the second element
981  // return NULL in the Elem2 field if there's no second element, i.e.
982  // the face is on the "boundary"
983  FaceElemTr.Elem2No = face_info.Elem2No;
984  if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
985  FaceElemTr.Elem2No >= 0)
986  {
987 #ifdef MFEM_DEBUG
988  if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
989  { MFEM_ABORT("NURBS mesh not supported!"); }
990 #endif
991  GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
992  FaceElemTr.Elem2 = &Transformation2;
993  cmask |= 2;
994  }
995 
996  // setup the face transformation
997  if (mask & FaceElementTransformations::HAVE_FACE)
998  {
999  GetFaceTransformation(FaceNo, &FaceElemTr);
1000  cmask |= 16;
1001  }
1002  else
1003  {
1004  FaceElemTr.SetGeometryType(GetFaceGeometry(FaceNo));
1005  }
1006 
1007  // setup Loc1 & Loc2
1008  int face_type = GetFaceElementType(FaceNo);
1009  if (mask & FaceElementTransformations::HAVE_LOC1)
1010  {
1011  int elem_type = GetElementType(face_info.Elem1No);
1012  GetLocalFaceTransformation(face_type, elem_type,
1013  FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
1014  cmask |= 4;
1015  }
1016  if ((mask & FaceElementTransformations::HAVE_LOC2) &&
1017  FaceElemTr.Elem2No >= 0)
1018  {
1019  int elem_type = GetElementType(face_info.Elem2No);
1020  GetLocalFaceTransformation(face_type, elem_type,
1021  FaceElemTr.Loc2.Transf, face_info.Elem2Inf);
1022 
1023  // NC meshes: prepend slave edge/face transformation to Loc2
1024  if (Nonconforming() && IsSlaveFace(face_info))
1025  {
1026  ApplyLocalSlaveTransformation(FaceElemTr, face_info, false);
1027  }
1028  cmask |= 8;
1029  }
1030 
1031  FaceElemTr.SetConfigurationMask(cmask);
1032 
1033  // This check can be useful for internal debugging, however it will fail on
1034  // periodic boundary faces, so we keep it disabled in general.
1035 #if 0
1036 #ifdef MFEM_DEBUG
1037  double dist = FaceElemTr.CheckConsistency();
1038  if (dist >= 1e-12)
1039  {
1040  mfem::out << "\nInternal error: face id = " << FaceNo
1041  << ", dist = " << dist << '\n';
1042  FaceElemTr.CheckConsistency(1); // print coordinates
1043  MFEM_ABORT("internal error");
1044  }
1045 #endif
1046 #endif
1047 
1048  return &FaceElemTr;
1049 }
1050 
1051 bool Mesh::IsSlaveFace(const FaceInfo &fi) const
1052 {
1053  return fi.NCFace >= 0 && nc_faces_info[fi.NCFace].Slave;
1054 }
1055 
1056 void Mesh::ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
1057  const FaceInfo &fi, bool is_ghost)
1058 {
1059 #ifdef MFEM_THREAD_SAFE
1060  DenseMatrix composition;
1061 #else
1062  static DenseMatrix composition;
1063 #endif
1064  MFEM_ASSERT(fi.NCFace >= 0, "");
1065  MFEM_ASSERT(nc_faces_info[fi.NCFace].Slave, "internal error");
1066  if (!is_ghost)
1067  {
1068  // side 1 -> child side, side 2 -> parent side
1070  LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
1071  // In 2D, we need to flip the point matrix since it is aligned with the
1072  // parent side.
1073  if (Dim == 2)
1074  {
1075  // swap points (columns) 0 and 1
1076  std::swap(composition(0,0), composition(0,1));
1077  std::swap(composition(1,0), composition(1,1));
1078  }
1079  LT.SetPointMat(composition);
1080  }
1081  else // is_ghost == true
1082  {
1083  // side 1 -> parent side, side 2 -> child side
1085  LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
1086  // In 2D, there is no need to flip the point matrix since it is already
1087  // aligned with the parent side, see also ParNCMesh::GetFaceNeighbors.
1088  // In 3D the point matrix was flipped during construction in
1089  // ParNCMesh::GetFaceNeighbors and due to that it is already aligned with
1090  // the parent side.
1091  LT.SetPointMat(composition);
1092  }
1093 }
1094 
1095 FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
1096 {
1098  int fn = GetBdrFace(BdrElemNo);
1099 
1100  // Check if the face is interior, shared, or nonconforming.
1101  if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1102  {
1103  return NULL;
1104  }
1105  tr = GetFaceElementTransformations(fn, 21);
1106  tr->Attribute = boundary[BdrElemNo]->GetAttribute();
1107  tr->ElementNo = BdrElemNo;
1108  tr->ElementType = ElementTransformation::BDR_FACE;
1109  tr->mesh = this;
1110  return tr;
1111 }
1112 
1113 int Mesh::GetBdrFace(int BdrElemNo) const
1114 {
1115  int fn;
1116  if (Dim == 3)
1117  {
1118  fn = be_to_face[BdrElemNo];
1119  }
1120  else if (Dim == 2)
1121  {
1122  fn = be_to_edge[BdrElemNo];
1123  }
1124  else
1125  {
1126  fn = boundary[BdrElemNo]->GetVertices()[0];
1127  }
1128  return fn;
1129 }
1130 
1131 Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
1132 {
1133  FaceInformation face;
1134  int e1, e2;
1135  int inf1, inf2;
1136  int ncface;
1137  GetFaceElements(f, &e1, &e2);
1138  GetFaceInfos(f, &inf1, &inf2, &ncface);
1139  face.element[0].index = e1;
1140  face.element[0].location = ElementLocation::Local;
1141  face.element[0].orientation = inf1%64;
1142  face.element[0].local_face_id = inf1/64;
1143  face.element[1].local_face_id = inf2/64;
1144  face.ncface = ncface;
1145  face.point_matrix = nullptr;
1146  // The following figures out face.location, face.conformity,
1147  // face.element[1].index, and face.element[1].orientation.
1148  if (f < GetNumFaces()) // Non-ghost face
1149  {
1150  if (e2>=0)
1151  {
1152  if (ncface==-1)
1153  {
1154  face.tag = FaceInfoTag::LocalConforming;
1155  face.topology = FaceTopology::Conforming;
1156  face.element[1].location = ElementLocation::Local;
1157  face.element[0].conformity = ElementConformity::Coincident;
1158  face.element[1].conformity = ElementConformity::Coincident;
1159  face.element[1].index = e2;
1160  face.element[1].orientation = inf2%64;
1161  }
1162  else // ncface >= 0
1163  {
1164  face.tag = FaceInfoTag::LocalSlaveNonconforming;
1165  face.topology = FaceTopology::Nonconforming;
1166  face.element[1].location = ElementLocation::Local;
1167  face.element[0].conformity = ElementConformity::Coincident;
1168  face.element[1].conformity = ElementConformity::Superset;
1169  face.element[1].index = e2;
1170  MFEM_ASSERT(inf2%64==0, "unexpected slave face orientation.");
1171  face.element[1].orientation = inf2%64;
1172  face.point_matrix = nc_faces_info[ncface].PointMatrix;
1173  }
1174  }
1175  else // e2<0
1176  {
1177  if (ncface==-1)
1178  {
1179  if (inf2<0)
1180  {
1181  face.tag = FaceInfoTag::Boundary;
1182  face.topology = FaceTopology::Boundary;
1183  face.element[1].location = ElementLocation::NA;
1184  face.element[0].conformity = ElementConformity::Coincident;
1185  face.element[1].conformity = ElementConformity::NA;
1186  face.element[1].index = -1;
1187  face.element[1].orientation = -1;
1188  }
1189  else // inf2 >= 0
1190  {
1191  face.tag = FaceInfoTag::SharedConforming;
1192  face.topology = FaceTopology::Conforming;
1193  face.element[0].conformity = ElementConformity::Coincident;
1194  face.element[1].conformity = ElementConformity::Coincident;
1195  face.element[1].location = ElementLocation::FaceNbr;
1196  face.element[1].index = -1 - e2;
1197  face.element[1].orientation = inf2%64;
1198  }
1199  }
1200  else // ncface >= 0
1201  {
1202  if (inf2 < 0)
1203  {
1204  face.tag = FaceInfoTag::MasterNonconforming;
1205  face.topology = FaceTopology::Nonconforming;
1206  face.element[1].location = ElementLocation::NA;
1207  face.element[0].conformity = ElementConformity::Coincident;
1208  face.element[1].conformity = ElementConformity::Subset;
1209  face.element[1].index = -1;
1210  face.element[1].orientation = -1;
1211  }
1212  else
1213  {
1214  face.tag = FaceInfoTag::SharedSlaveNonconforming;
1215  face.topology = FaceTopology::Nonconforming;
1216  face.element[1].location = ElementLocation::FaceNbr;
1217  face.element[0].conformity = ElementConformity::Coincident;
1218  face.element[1].conformity = ElementConformity::Superset;
1219  face.element[1].index = -1 - e2;
1220  face.element[1].orientation = inf2%64;
1221  }
1222  face.point_matrix = nc_faces_info[ncface].PointMatrix;
1223  }
1224  }
1225  }
1226  else // Ghost face
1227  {
1228  if (e1==-1)
1229  {
1230  face.tag = FaceInfoTag::GhostMaster;
1231  face.topology = FaceTopology::NA;
1232  face.element[1].location = ElementLocation::NA;
1233  face.element[0].conformity = ElementConformity::NA;
1234  face.element[1].conformity = ElementConformity::NA;
1235  face.element[1].index = -1;
1236  face.element[1].orientation = -1;
1237  }
1238  else
1239  {
1240  face.tag = FaceInfoTag::GhostSlave;
1241  face.topology = FaceTopology::Nonconforming;
1242  face.element[1].location = ElementLocation::FaceNbr;
1243  face.element[0].conformity = ElementConformity::Superset;
1244  face.element[1].conformity = ElementConformity::Coincident;
1245  face.element[1].index = -1 - e2;
1246  face.element[1].orientation = inf2%64;
1247  face.point_matrix = nc_faces_info[ncface].PointMatrix;
1248  }
1249  }
1250  return face;
1251 }
1252 
1253 Mesh::FaceInformation::operator Mesh::FaceInfo() const
1254 {
1255  FaceInfo res {-1, -1, -1, -1, -1};
1256  switch (tag)
1257  {
1258  case FaceInfoTag::LocalConforming:
1259  res.Elem1No = element[0].index;
1260  res.Elem2No = element[1].index;
1261  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1262  res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1263  res.NCFace = ncface;
1264  break;
1265  case FaceInfoTag::LocalSlaveNonconforming:
1266  res.Elem1No = element[0].index;
1267  res.Elem2No = element[1].index;
1268  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1269  res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1270  res.NCFace = ncface;
1271  break;
1272  case FaceInfoTag::Boundary:
1273  res.Elem1No = element[0].index;
1274  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1275  break;
1276  case FaceInfoTag::SharedConforming:
1277  res.Elem1No = element[0].index;
1278  res.Elem2No = -1 - element[1].index;
1279  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1280  res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1281  break;
1282  case FaceInfoTag::MasterNonconforming:
1283  res.Elem1No = element[0].index;
1284  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1285  break;
1286  case FaceInfoTag::SharedSlaveNonconforming:
1287  res.Elem1No = element[0].index;
1288  res.Elem2No = -1 - element[1].index;
1289  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1290  res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1291  break;
1292  case FaceInfoTag::GhostMaster:
1293  break;
1294  case FaceInfoTag::GhostSlave:
1295  res.Elem1No = element[0].index;
1296  res.Elem2No = -1 - element[1].index;
1297  res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1298  res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1299  break;
1300  }
1301  return res;
1302 }
1303 
1304 std::ostream& operator<<(std::ostream& os, const Mesh::FaceInformation& info)
1305 {
1306  os << "face topology=";
1307  switch (info.topology)
1308  {
1309  case Mesh::FaceTopology::Boundary:
1310  os << "Boundary";
1311  break;
1312  case Mesh::FaceTopology::Conforming:
1313  os << "Conforming";
1314  break;
1315  case Mesh::FaceTopology::Nonconforming:
1316  os << "Non-conforming";
1317  break;
1318  case Mesh::FaceTopology::NA:
1319  os << "NA";
1320  break;
1321  }
1322  os << "element[0].location=";
1323  switch (info.element[0].location)
1324  {
1325  case Mesh::ElementLocation::Local:
1326  os << "Local";
1327  break;
1328  case Mesh::ElementLocation::FaceNbr:
1329  os << "FaceNbr";
1330  break;
1331  case Mesh::ElementLocation::NA:
1332  os << "NA";
1333  break;
1334  }
1335  os << std::endl;
1336  os << "element[1].location=";
1337  switch (info.element[1].location)
1338  {
1339  case Mesh::ElementLocation::Local:
1340  os << "Local";
1341  break;
1342  case Mesh::ElementLocation::FaceNbr:
1343  os << "FaceNbr";
1344  break;
1345  case Mesh::ElementLocation::NA:
1346  os << "NA";
1347  break;
1348  }
1349  os << std::endl;
1350  os << "element[0].conformity=";
1351  switch (info.element[0].conformity)
1352  {
1353  case Mesh::ElementConformity::Coincident:
1354  os << "Coincident";
1355  break;
1356  case Mesh::ElementConformity::Superset:
1357  os << "Superset";
1358  break;
1359  case Mesh::ElementConformity::Subset:
1360  os << "Subset";
1361  break;
1362  case Mesh::ElementConformity::NA:
1363  os << "NA";
1364  break;
1365  }
1366  os << std::endl;
1367  os << "element[1].conformity=";
1368  switch (info.element[1].conformity)
1369  {
1370  case Mesh::ElementConformity::Coincident:
1371  os << "Coincident";
1372  break;
1373  case Mesh::ElementConformity::Superset:
1374  os << "Superset";
1375  break;
1376  case Mesh::ElementConformity::Subset:
1377  os << "Subset";
1378  break;
1379  case Mesh::ElementConformity::NA:
1380  os << "NA";
1381  break;
1382  }
1383  os << std::endl;
1384  os << "element[0].index=" << info.element[0].index << std::endl
1385  << "element[1].index=" << info.element[1].index << std::endl
1386  << "element[0].local_face_id=" << info.element[0].local_face_id << std::endl
1387  << "element[1].local_face_id=" << info.element[1].local_face_id << std::endl
1388  << "element[0].orientation=" << info.element[0].orientation << std::endl
1389  << "element[1].orientation=" << info.element[1].orientation << std::endl
1390  << "ncface=" << info.ncface << std::endl;
1391  return os;
1392 }
1393 
1394 void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2) const
1395 {
1396  *Elem1 = faces_info[Face].Elem1No;
1397  *Elem2 = faces_info[Face].Elem2No;
1398 }
1399 
1400 void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2) const
1401 {
1402  *Inf1 = faces_info[Face].Elem1Inf;
1403  *Inf2 = faces_info[Face].Elem2Inf;
1404 }
1405 
1406 void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2, int *NCFace) const
1407 {
1408  *Inf1 = faces_info[Face].Elem1Inf;
1409  *Inf2 = faces_info[Face].Elem2Inf;
1410  *NCFace = faces_info[Face].NCFace;
1411 }
1412 
1413 Geometry::Type Mesh::GetFaceGeometry(int Face) const
1414 {
1415  switch (Dim)
1416  {
1417  case 1: return Geometry::POINT;
1418  case 2: return Geometry::SEGMENT;
1419  case 3:
1420  if (Face < NumOfFaces) // local (non-ghost) face
1421  {
1422  return faces[Face]->GetGeometryType();
1423  }
1424  // ghost face
1425  const int nc_face_id = faces_info[Face].NCFace;
1426  MFEM_ASSERT(nc_face_id >= 0, "parent ghost faces are not supported");
1427  return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1428  }
1429  return Geometry::INVALID;
1430 }
1431 
1432 Element::Type Mesh::GetFaceElementType(int Face) const
1433 {
1434  return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1435 }
1436 
1437 Array<int> Mesh::GetFaceToBdrElMap() const
1438 {
1439  Array<int> face_to_be(NumOfFaces);
1440  face_to_be = -1;
1441  for (int i = 0; i < NumOfBdrElements; i++)
1442  {
1443  face_to_be[GetBdrElementEdgeIndex(i)] = i;
1444  }
1445  return face_to_be;
1446 }
1447 
1448 void Mesh::Init()
1449 {
1450  // in order of declaration:
1451  Dim = spaceDim = 0;
1452  NumOfVertices = -1;
1453  NumOfElements = NumOfBdrElements = 0;
1454  NumOfEdges = NumOfFaces = 0;
1455  nbInteriorFaces = -1;
1456  nbBoundaryFaces = -1;
1457  meshgen = mesh_geoms = 0;
1458  sequence = 0;
1459  Nodes = NULL;
1460  own_nodes = 1;
1461  NURBSext = NULL;
1462  ncmesh = NULL;
1463  last_operation = Mesh::NONE;
1464 }
1465 
1466 void Mesh::InitTables()
1467 {
1468  el_to_edge =
1469  el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1470  face_to_elem = NULL;
1471 }
1472 
1473 void Mesh::SetEmpty()
1474 {
1475  Init();
1476  InitTables();
1477 }
1478 
1479 void Mesh::DestroyTables()
1480 {
1481  delete el_to_edge;
1482  delete el_to_face;
1483  delete el_to_el;
1484  DeleteGeometricFactors();
1485 
1486  if (Dim == 3)
1487  {
1488  delete bel_to_edge;
1489  }
1490 
1491  delete face_edge;
1492  delete edge_vertex;
1493 
1494  delete face_to_elem;
1495  face_to_elem = NULL;
1496 }
1497 
1498 void Mesh::DestroyPointers()
1499 {
1500  if (own_nodes) { delete Nodes; }
1501 
1502  delete ncmesh;
1503 
1504  delete NURBSext;
1505 
1506  for (int i = 0; i < NumOfElements; i++)
1507  {
1508  FreeElement(elements[i]);
1509  }
1510 
1511  for (int i = 0; i < NumOfBdrElements; i++)
1512  {
1513  FreeElement(boundary[i]);
1514  }
1515 
1516  for (int i = 0; i < faces.Size(); i++)
1517  {
1518  FreeElement(faces[i]);
1519  }
1520 
1521  DestroyTables();
1522 }
1523 
1524 void Mesh::Destroy()
1525 {
1526  DestroyPointers();
1527 
1528  elements.DeleteAll();
1529  vertices.DeleteAll();
1530  boundary.DeleteAll();
1531  faces.DeleteAll();
1532  faces_info.DeleteAll();
1533  nc_faces_info.DeleteAll();
1534  be_to_edge.DeleteAll();
1535  be_to_face.DeleteAll();
1536 
1537  // TODO:
1538  // IsoparametricTransformations
1539  // Transformation, Transformation2, BdrTransformation, FaceTransformation,
1540  // EdgeTransformation;
1541  // FaceElementTransformations FaceElemTr;
1542 
1543  CoarseFineTr.Clear();
1544 
1545 #ifdef MFEM_USE_MEMALLOC
1546  TetMemory.Clear();
1547 #endif
1548 
1549  attributes.DeleteAll();
1550  bdr_attributes.DeleteAll();
1551 }
1552 
1553 void Mesh::ResetLazyData()
1554 {
1555  delete el_to_el; el_to_el = NULL;
1556  delete face_edge; face_edge = NULL;
1557  delete face_to_elem; face_to_elem = NULL;
1558  delete edge_vertex; edge_vertex = NULL;
1559  DeleteGeometricFactors();
1560  nbInteriorFaces = -1;
1561  nbBoundaryFaces = -1;
1562 }
1563 
1564 void Mesh::SetAttributes()
1565 {
1566  Array<int> attribs;
1567 
1568  attribs.SetSize(GetNBE());
1569  for (int i = 0; i < attribs.Size(); i++)
1570  {
1571  attribs[i] = GetBdrAttribute(i);
1572  }
1573  attribs.Sort();
1574  attribs.Unique();
1575  attribs.Copy(bdr_attributes);
1576  if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1577  {
1578  MFEM_WARNING("Non-positive attributes on the boundary!");
1579  }
1580 
1581  attribs.SetSize(GetNE());
1582  for (int i = 0; i < attribs.Size(); i++)
1583  {
1584  attribs[i] = GetAttribute(i);
1585  }
1586  attribs.Sort();
1587  attribs.Unique();
1588  attribs.Copy(attributes);
1589  if (attributes.Size() > 0 && attributes[0] <= 0)
1590  {
1591  MFEM_WARNING("Non-positive attributes in the domain!");
1592  }
1593 }
1594 
1595 void Mesh::InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
1596 {
1597  SetEmpty();
1598 
1599  Dim = Dim_;
1600  spaceDim = spaceDim_;
1601 
1602  NumOfVertices = 0;
1603  vertices.SetSize(NVert); // just allocate space for vertices
1604 
1605  NumOfElements = 0;
1606  elements.SetSize(NElem); // just allocate space for Element *
1607 
1608  NumOfBdrElements = 0;
1609  boundary.SetSize(NBdrElem); // just allocate space for Element *
1610 }
1611 
1612 template<typename T>
1613 static void CheckEnlarge(Array<T> &array, int size)
1614 {
1615  if (size >= array.Size()) { array.SetSize(size + 1); }
1616 }
1617 
1618 int Mesh::AddVertex(double x, double y, double z)
1619 {
1620  CheckEnlarge(vertices, NumOfVertices);
1621  double *v = vertices[NumOfVertices]();
1622  v[0] = x;
1623  v[1] = y;
1624  v[2] = z;
1625  return NumOfVertices++;
1626 }
1627 
1628 int Mesh::AddVertex(const double *coords)
1629 {
1630  CheckEnlarge(vertices, NumOfVertices);
1631  vertices[NumOfVertices].SetCoords(spaceDim, coords);
1632  return NumOfVertices++;
1633 }
1634 
1635 int Mesh::AddVertex(const Vector &coords)
1636 {
1637  MFEM_ASSERT(coords.Size() >= spaceDim,
1638  "invalid 'coords' size: " << coords.Size());
1639  return AddVertex(coords.GetData());
1640 }
1641 
1642 void Mesh::AddVertexParents(int i, int p1, int p2)
1643 {
1644  tmp_vertex_parents.Append(Triple<int, int, int>(i, p1, p2));
1645 
1646  // if vertex coordinates are defined, make sure the hanging vertex has the
1647  // correct position
1648  if (i < vertices.Size())
1649  {
1650  double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1651  for (int j = 0; j < 3; j++)
1652  {
1653  vi[j] = (vp1[j] + vp2[j]) * 0.5;
1654  }
1655  }
1656 }
1657 
1658 int Mesh::AddSegment(int v1, int v2, int attr)
1659 {
1660  CheckEnlarge(elements, NumOfElements);
1661  elements[NumOfElements] = new Segment(v1, v2, attr);
1662  return NumOfElements++;
1663 }
1664 
1665 int Mesh::AddSegment(const int *vi, int attr)
1666 {
1667  CheckEnlarge(elements, NumOfElements);
1668  elements[NumOfElements] = new Segment(vi, attr);
1669  return NumOfElements++;
1670 }
1671 
1672 int Mesh::AddTriangle(int v1, int v2, int v3, int attr)
1673 {
1674  CheckEnlarge(elements, NumOfElements);
1675  elements[NumOfElements] = new Triangle(v1, v2, v3, attr);
1676  return NumOfElements++;
1677 }
1678 
1679 int Mesh::AddTriangle(const int *vi, int attr)
1680 {
1681  CheckEnlarge(elements, NumOfElements);
1682  elements[NumOfElements] = new Triangle(vi, attr);
1683  return NumOfElements++;
1684 }
1685 
1686 int Mesh::AddQuad(int v1, int v2, int v3, int v4, int attr)
1687 {
1688  CheckEnlarge(elements, NumOfElements);
1689  elements[NumOfElements] = new Quadrilateral(v1, v2, v3, v4, attr);
1690  return NumOfElements++;
1691 }
1692 
1693 int Mesh::AddQuad(const int *vi, int attr)
1694 {
1695  CheckEnlarge(elements, NumOfElements);
1696  elements[NumOfElements] = new Quadrilateral(vi, attr);
1697  return NumOfElements++;
1698 }
1699 
1700 int Mesh::AddTet(int v1, int v2, int v3, int v4, int attr)
1701 {
1702  int vi[4] = {v1, v2, v3, v4};
1703  return AddTet(vi, attr);
1704 }
1705 
1706 int Mesh::AddTet(const int *vi, int attr)
1707 {
1708  CheckEnlarge(elements, NumOfElements);
1709 #ifdef MFEM_USE_MEMALLOC
1710  Tetrahedron *tet;
1711  tet = TetMemory.Alloc();
1712  tet->SetVertices(vi);
1713  tet->SetAttribute(attr);
1714  elements[NumOfElements] = tet;
1715 #else
1716  elements[NumOfElements] = new Tetrahedron(vi, attr);
1717 #endif
1718  return NumOfElements++;
1719 }
1720 
1721 int Mesh::AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr)
1722 {
1723  CheckEnlarge(elements, NumOfElements);
1724  elements[NumOfElements] = new Wedge(v1, v2, v3, v4, v5, v6, attr);
1725  return NumOfElements++;
1726 }
1727 
1728 int Mesh::AddWedge(const int *vi, int attr)
1729 {
1730  CheckEnlarge(elements, NumOfElements);
1731  elements[NumOfElements] = new Wedge(vi, attr);
1732  return NumOfElements++;
1733 }
1734 
1735 int Mesh::AddPyramid(int v1, int v2, int v3, int v4, int v5, int attr)
1736 {
1737  CheckEnlarge(elements, NumOfElements);
1738  elements[NumOfElements] = new Pyramid(v1, v2, v3, v4, v5, attr);
1739  return NumOfElements++;
1740 }
1741 
1742 int Mesh::AddPyramid(const int *vi, int attr)
1743 {
1744  CheckEnlarge(elements, NumOfElements);
1745  elements[NumOfElements] = new Pyramid(vi, attr);
1746  return NumOfElements++;
1747 }
1748 
1749 int Mesh::AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8,
1750  int attr)
1751 {
1752  CheckEnlarge(elements, NumOfElements);
1753  elements[NumOfElements] =
1754  new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1755  return NumOfElements++;
1756 }
1757 
1758 int Mesh::AddHex(const int *vi, int attr)
1759 {
1760  CheckEnlarge(elements, NumOfElements);
1761  elements[NumOfElements] = new Hexahedron(vi, attr);
1762  return NumOfElements++;
1763 }
1764 
1765 void Mesh::AddHexAsTets(const int *vi, int attr)
1766 {
1767  static const int hex_to_tet[6][4] =
1768  {
1769  { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1770  { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1771  };
1772  int ti[4];
1773 
1774  for (int i = 0; i < 6; i++)
1775  {
1776  for (int j = 0; j < 4; j++)
1777  {
1778  ti[j] = vi[hex_to_tet[i][j]];
1779  }
1780  AddTet(ti, attr);
1781  }
1782 }
1783 
1784 void Mesh::AddHexAsWedges(const int *vi, int attr)
1785 {
1786  static const int hex_to_wdg[2][6] =
1787  {
1788  { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1789  };
1790  int ti[6];
1791 
1792  for (int i = 0; i < 2; i++)
1793  {
1794  for (int j = 0; j < 6; j++)
1795  {
1796  ti[j] = vi[hex_to_wdg[i][j]];
1797  }
1798  AddWedge(ti, attr);
1799  }
1800 }
1801 
1802 void Mesh::AddHexAsPyramids(const int *vi, int attr)
1803 {
1804  static const int hex_to_pyr[6][5] =
1805  {
1806  { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1807  { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1808  };
1809  int ti[5];
1810 
1811  for (int i = 0; i < 6; i++)
1812  {
1813  for (int j = 0; j < 5; j++)
1814  {
1815  ti[j] = vi[hex_to_pyr[i][j]];
1816  }
1817  AddPyramid(ti, attr);
1818  }
1819 }
1820 
1821 int Mesh::AddElement(Element *elem)
1822 {
1823  CheckEnlarge(elements, NumOfElements);
1824  elements[NumOfElements] = elem;
1825  return NumOfElements++;
1826 }
1827 
1828 int Mesh::AddBdrElement(Element *elem)
1829 {
1830  CheckEnlarge(boundary, NumOfBdrElements);
1831  boundary[NumOfBdrElements] = elem;
1832  return NumOfBdrElements++;
1833 }
1834 
1835 int Mesh::AddBdrSegment(int v1, int v2, int attr)
1836 {
1837  CheckEnlarge(boundary, NumOfBdrElements);
1838  boundary[NumOfBdrElements] = new Segment(v1, v2, attr);
1839  return NumOfBdrElements++;
1840 }
1841 
1842 int Mesh::AddBdrSegment(const int *vi, int attr)
1843 {
1844  CheckEnlarge(boundary, NumOfBdrElements);
1845  boundary[NumOfBdrElements] = new Segment(vi, attr);
1846  return NumOfBdrElements++;
1847 }
1848 
1849 int Mesh::AddBdrTriangle(int v1, int v2, int v3, int attr)
1850 {
1851  CheckEnlarge(boundary, NumOfBdrElements);
1852  boundary[NumOfBdrElements] = new Triangle(v1, v2, v3, attr);
1853  return NumOfBdrElements++;
1854 }
1855 
1856 int Mesh::AddBdrTriangle(const int *vi, int attr)
1857 {
1858  CheckEnlarge(boundary, NumOfBdrElements);
1859  boundary[NumOfBdrElements] = new Triangle(vi, attr);
1860  return NumOfBdrElements++;
1861 }
1862 
1863 int Mesh::AddBdrQuad(int v1, int v2, int v3, int v4, int attr)
1864 {
1865  CheckEnlarge(boundary, NumOfBdrElements);
1866  boundary[NumOfBdrElements] = new Quadrilateral(v1, v2, v3, v4, attr);
1867  return NumOfBdrElements++;
1868 }
1869 
1870 int Mesh::AddBdrQuad(const int *vi, int attr)
1871 {
1872  CheckEnlarge(boundary, NumOfBdrElements);
1873  boundary[NumOfBdrElements] = new Quadrilateral(vi, attr);
1874  return NumOfBdrElements++;
1875 }
1876 
1877 void Mesh::AddBdrQuadAsTriangles(const int *vi, int attr)
1878 {
1879  static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1880  int ti[3];
1881 
1882  for (int i = 0; i < 2; i++)
1883  {
1884  for (int j = 0; j < 3; j++)
1885  {
1886  ti[j] = vi[quad_to_tri[i][j]];
1887  }
1888  AddBdrTriangle(ti, attr);
1889  }
1890 }
1891 
1892 int Mesh::AddBdrPoint(int v, int attr)
1893 {
1894  CheckEnlarge(boundary, NumOfBdrElements);
1895  boundary[NumOfBdrElements] = new Point(&v, attr);
1896  return NumOfBdrElements++;
1897 }
1898 
1899 void Mesh::GenerateBoundaryElements()
1900 {
1901  int i, j;
1902  Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1903 
1904  // GenerateFaces();
1905 
1906  for (i = 0; i < boundary.Size(); i++)
1907  {
1908  FreeElement(boundary[i]);
1909  }
1910 
1911  if (Dim == 3)
1912  {
1913  delete bel_to_edge;
1914  bel_to_edge = NULL;
1915  }
1916 
1917  // count the 'NumOfBdrElements'
1918  NumOfBdrElements = 0;
1919  for (i = 0; i < faces_info.Size(); i++)
1920  {
1921  if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1922  }
1923 
1924  boundary.SetSize(NumOfBdrElements);
1925  be2face.SetSize(NumOfBdrElements);
1926  for (j = i = 0; i < faces_info.Size(); i++)
1927  {
1928  if (faces_info[i].Elem2No < 0)
1929  {
1930  boundary[j] = faces[i]->Duplicate(this);
1931  be2face[j++] = i;
1932  }
1933  }
1934  // In 3D, 'bel_to_edge' is destroyed but it's not updated.
1935 }
1936 
1937 void Mesh::FinalizeCheck()
1938 {
1939  MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1940  vertices.Size() == 0,
1941  "incorrect number of vertices: preallocated: " << vertices.Size()
1942  << ", actually added: " << NumOfVertices);
1943  MFEM_VERIFY(elements.Size() == NumOfElements,
1944  "incorrect number of elements: preallocated: " << elements.Size()
1945  << ", actually added: " << NumOfElements);
1946  MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1947  "incorrect number of boundary elements: preallocated: "
1948  << boundary.Size() << ", actually added: " << NumOfBdrElements);
1949 }
1950 
1951 void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
1952 {
1953  FinalizeCheck();
1954  CheckElementOrientation(fix_orientation);
1955 
1956  if (refine)
1957  {
1958  MarkTriMeshForRefinement();
1959  }
1960 
1961  if (generate_edges)
1962  {
1963  el_to_edge = new Table;
1964  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1965  GenerateFaces();
1966  CheckBdrElementOrientation();
1967  }
1968  else
1969  {
1970  NumOfEdges = 0;
1971  }
1972 
1973  NumOfFaces = 0;
1974 
1975  SetAttributes();
1976 
1977  SetMeshGen();
1978 }
1979 
1980 void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
1981  bool fix_orientation)
1982 {
1983  FinalizeCheck();
1984  if (fix_orientation)
1985  {
1986  CheckElementOrientation(fix_orientation);
1987  }
1988 
1989  if (generate_edges)
1990  {
1991  el_to_edge = new Table;
1992  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1993  GenerateFaces();
1994  CheckBdrElementOrientation();
1995  }
1996  else
1997  {
1998  NumOfEdges = 0;
1999  }
2000 
2001  NumOfFaces = 0;
2002 
2003  SetAttributes();
2004 
2005  SetMeshGen();
2006 }
2007 
2008 
2009 class GeckoProgress : public Gecko::Progress
2010 {
2011  double limit;
2012  mutable StopWatch sw;
2013 public:
2014  GeckoProgress(double limit) : limit(limit) { sw.Start(); }
2015  virtual bool quit() const { return limit > 0 && sw.UserTime() > limit; }
2016 };
2017 
2018 class GeckoVerboseProgress : public GeckoProgress
2019 {
2020  using Float = Gecko::Float;
2021  using Graph = Gecko::Graph;
2022  using uint = Gecko::uint;
2023 public:
2024  GeckoVerboseProgress(double limit) : GeckoProgress(limit) {}
2025 
2026  virtual void beginorder(const Graph* graph, Float cost) const
2027  { mfem::out << "Begin Gecko ordering, cost = " << cost << std::endl; }
2028  virtual void endorder(const Graph* graph, Float cost) const
2029  { mfem::out << "End ordering, cost = " << cost << std::endl; }
2030 
2031  virtual void beginiter(const Graph* graph,
2032  uint iter, uint maxiter, uint window) const
2033  {
2034  mfem::out << "Iteration " << iter << "/" << maxiter << ", window "
2035  << window << std::flush;
2036  }
2037  virtual void enditer(const Graph* graph, Float mincost, Float cost) const
2038  { mfem::out << ", cost = " << cost << endl; }
2039 };
2040 
2041 
2042 double Mesh::GetGeckoElementOrdering(Array<int> &ordering,
2043  int iterations, int window,
2044  int period, int seed, bool verbose,
2045  double time_limit)
2046 {
2047  Gecko::Graph graph;
2048  Gecko::FunctionalGeometric functional; // edge product cost
2049 
2050  GeckoProgress progress(time_limit);
2051  GeckoVerboseProgress vprogress(time_limit);
2052 
2053  // insert elements as nodes in the graph
2054  for (int elemid = 0; elemid < GetNE(); ++elemid)
2055  {
2056  graph.insert_node();
2057  }
2058 
2059  // insert graph edges for element neighbors
2060  // NOTE: indices in Gecko are 1 based hence the +1 on insertion
2061  const Table &my_el_to_el = ElementToElementTable();
2062  for (int elemid = 0; elemid < GetNE(); ++elemid)
2063  {
2064  const int *neighid = my_el_to_el.GetRow(elemid);
2065  for (int i = 0; i < my_el_to_el.RowSize(elemid); ++i)
2066  {
2067  graph.insert_arc(elemid + 1, neighid[i] + 1);
2068  }
2069  }
2070 
2071  // get the ordering from Gecko and copy it into the Array<int>
2072  graph.order(&functional, iterations, window, period, seed,
2073  verbose ? &vprogress : &progress);
2074 
2075  ordering.SetSize(GetNE());
2076  Gecko::Node::Index NE = GetNE();
2077  for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
2078  {
2079  ordering[gnodeid - 1] = graph.rank(gnodeid);
2080  }
2081 
2082  return graph.cost();
2083 }
2084 
2085 
2086 struct HilbertCmp
2087 {
2088  int coord;
2089  bool dir;
2090  const Array<double> &points;
2091  double mid;
2092 
2093  HilbertCmp(int coord, bool dir, const Array<double> &points, double mid)
2094  : coord(coord), dir(dir), points(points), mid(mid) {}
2095 
2096  bool operator()(int i) const
2097  {
2098  return (points[3*i + coord] < mid) != dir;
2099  }
2100 };
2101 
2102 static void HilbertSort2D(int coord1, // major coordinate to sort points by
2103  bool dir1, // sort coord1 ascending/descending?
2104  bool dir2, // sort coord2 ascending/descending?
2105  const Array<double> &points, int *beg, int *end,
2106  double xmin, double ymin, double xmax, double ymax)
2107 {
2108  if (end - beg <= 1) { return; }
2109 
2110  double xmid = (xmin + xmax)*0.5;
2111  double ymid = (ymin + ymax)*0.5;
2112 
2113  int coord2 = (coord1 + 1) % 2; // the 'other' coordinate
2114 
2115  // sort (partition) points into four quadrants
2116  int *p0 = beg, *p4 = end;
2117  int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2118  int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2119  int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2120 
2121  if (p1 != p4)
2122  {
2123  HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2124  ymin, xmin, ymid, xmid);
2125  }
2126  if (p1 != p0 || p2 != p4)
2127  {
2128  HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2129  xmin, ymid, xmid, ymax);
2130  }
2131  if (p2 != p0 || p3 != p4)
2132  {
2133  HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2134  xmid, ymid, xmax, ymax);
2135  }
2136  if (p3 != p0)
2137  {
2138  HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2139  ymid, xmax, ymin, xmid);
2140  }
2141 }
2142 
2143 static void HilbertSort3D(int coord1, bool dir1, bool dir2, bool dir3,
2144  const Array<double> &points, int *beg, int *end,
2145  double xmin, double ymin, double zmin,
2146  double xmax, double ymax, double zmax)
2147 {
2148  if (end - beg <= 1) { return; }
2149 
2150  double xmid = (xmin + xmax)*0.5;
2151  double ymid = (ymin + ymax)*0.5;
2152  double zmid = (zmin + zmax)*0.5;
2153 
2154  int coord2 = (coord1 + 1) % 3;
2155  int coord3 = (coord1 + 2) % 3;
2156 
2157  // sort (partition) points into eight octants
2158  int *p0 = beg, *p8 = end;
2159  int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2160  int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2161  int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2162  int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2163  int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2164  int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2165  int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2166 
2167  if (p1 != p8)
2168  {
2169  HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2170  zmin, xmin, ymin, zmid, xmid, ymid);
2171  }
2172  if (p1 != p0 || p2 != p8)
2173  {
2174  HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2175  ymin, zmid, xmin, ymid, zmax, xmid);
2176  }
2177  if (p2 != p0 || p3 != p8)
2178  {
2179  HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2180  ymid, zmid, xmin, ymax, zmax, xmid);
2181  }
2182  if (p3 != p0 || p4 != p8)
2183  {
2184  HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2185  xmin, ymax, zmid, xmid, ymid, zmin);
2186  }
2187  if (p4 != p0 || p5 != p8)
2188  {
2189  HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2190  xmid, ymax, zmid, xmax, ymid, zmin);
2191  }
2192  if (p5 != p0 || p6 != p8)
2193  {
2194  HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2195  ymax, zmid, xmax, ymid, zmax, xmid);
2196  }
2197  if (p6 != p0 || p7 != p8)
2198  {
2199  HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2200  ymid, zmid, xmax, ymin, zmax, xmid);
2201  }
2202  if (p7 != p0)
2203  {
2204  HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2205  zmid, xmax, ymin, zmin, xmid, ymid);
2206  }
2207 }
2208 
2209 void Mesh::GetHilbertElementOrdering(Array<int> &ordering)
2210 {
2211  MFEM_VERIFY(spaceDim <= 3, "");
2212 
2213  Vector min, max, center;
2214  GetBoundingBox(min, max);
2215 
2216  Array<int> indices(GetNE());
2217  Array<double> points(3*GetNE());
2218 
2219  if (spaceDim < 3) { points = 0.0; }
2220 
2221  // calculate element centers
2222  for (int i = 0; i < GetNE(); i++)
2223  {
2224  GetElementCenter(i, center);
2225  for (int j = 0; j < spaceDim; j++)
2226  {
2227  points[3*i + j] = center(j);
2228  }
2229  indices[i] = i;
2230  }
2231 
2232  if (spaceDim == 1)
2233  {
2234  indices.Sort([&](int a, int b)
2235  { return points[3*a] < points[3*b]; });
2236  }
2237  else if (spaceDim == 2)
2238  {
2239  // recursively partition the points in 2D
2240  HilbertSort2D(0, false, false,
2241  points, indices.begin(), indices.end(),
2242  min(0), min(1), max(0), max(1));
2243  }
2244  else
2245  {
2246  // recursively partition the points in 3D
2247  HilbertSort3D(0, false, false, false,
2248  points, indices.begin(), indices.end(),
2249  min(0), min(1), min(2), max(0), max(1), max(2));
2250  }
2251 
2252  // return ordering in the format required by ReorderElements
2253  ordering.SetSize(GetNE());
2254  for (int i = 0; i < GetNE(); i++)
2255  {
2256  ordering[indices[i]] = i;
2257  }
2258 }
2259 
2260 
2261 void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
2262 {
2263  if (NURBSext)
2264  {
2265  MFEM_WARNING("element reordering of NURBS meshes is not supported.");
2266  return;
2267  }
2268  if (ncmesh)
2269  {
2270  MFEM_WARNING("element reordering of non-conforming meshes is not"
2271  " supported.");
2272  return;
2273  }
2274  MFEM_VERIFY(ordering.Size() == GetNE(), "invalid reordering array.")
2275 
2276  // Data members that need to be updated:
2277 
2278  // - elements - reorder of the pointers and the vertex ids if reordering
2279  // the vertices
2280  // - vertices - if reordering the vertices
2281  // - boundary - update the vertex ids, if reordering the vertices
2282  // - faces - regenerate
2283  // - faces_info - regenerate
2284 
2285  // Deleted by DeleteTables():
2286  // - el_to_edge - rebuild in 2D and 3D only
2287  // - el_to_face - rebuild in 3D only
2288  // - bel_to_edge - rebuild in 3D only
2289  // - el_to_el - no need to rebuild
2290  // - face_edge - no need to rebuild
2291  // - edge_vertex - no need to rebuild
2292  // - geom_factors - no need to rebuild
2293 
2294  // - be_to_edge - 2D only
2295  // - be_to_face - 3D only
2296 
2297  // - Nodes
2298 
2299  // Save the locations of the Nodes so we can rebuild them later
2300  Array<Vector*> old_elem_node_vals;
2301  FiniteElementSpace *nodes_fes = NULL;
2302  if (Nodes)
2303  {
2304  old_elem_node_vals.SetSize(GetNE());
2305  nodes_fes = Nodes->FESpace();
2306  Array<int> old_dofs;
2307  Vector vals;
2308  for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
2309  {
2310  nodes_fes->GetElementVDofs(old_elid, old_dofs);
2311  Nodes->GetSubVector(old_dofs, vals);
2312  old_elem_node_vals[old_elid] = new Vector(vals);
2313  }
2314  }
2315 
2316  // Get the newly ordered elements
2317  Array<Element *> new_elements(GetNE());
2318  for (int old_elid = 0; old_elid < ordering.Size(); ++old_elid)
2319  {
2320  int new_elid = ordering[old_elid];
2321  new_elements[new_elid] = elements[old_elid];
2322  }
2323  mfem::Swap(elements, new_elements);
2324  new_elements.DeleteAll();
2325 
2326  if (reorder_vertices)
2327  {
2328  // Get the new vertex ordering permutation vectors and fill the new
2329  // vertices
2330  Array<int> vertex_ordering(GetNV());
2331  vertex_ordering = -1;
2332  Array<Vertex> new_vertices(GetNV());
2333  int new_vertex_ind = 0;
2334  for (int new_elid = 0; new_elid < GetNE(); ++new_elid)
2335  {
2336  int *elem_vert = elements[new_elid]->GetVertices();
2337  int nv = elements[new_elid]->GetNVertices();
2338  for (int vi = 0; vi < nv; ++vi)
2339  {
2340  int old_vertex_ind = elem_vert[vi];
2341  if (vertex_ordering[old_vertex_ind] == -1)
2342  {
2343  vertex_ordering[old_vertex_ind] = new_vertex_ind;
2344  new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
2345  new_vertex_ind++;
2346  }
2347  }
2348  }
2349  mfem::Swap(vertices, new_vertices);
2350  new_vertices.DeleteAll();
2351 
2352  // Replace the vertex ids in the elements with the reordered vertex
2353  // numbers
2354  for (int new_elid = 0; new_elid < GetNE(); ++new_elid)
2355  {
2356  int *elem_vert = elements[new_elid]->GetVertices();
2357  int nv = elements[new_elid]->GetNVertices();
2358  for (int vi = 0; vi < nv; ++vi)
2359  {
2360  elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2361  }
2362  }
2363 
2364  // Replace the vertex ids in the boundary with reordered vertex numbers
2365  for (int belid = 0; belid < GetNBE(); ++belid)
2366  {
2367  int *be_vert = boundary[belid]->GetVertices();
2368  int nv = boundary[belid]->GetNVertices();
2369  for (int vi = 0; vi < nv; ++vi)
2370  {
2371  be_vert[vi] = vertex_ordering[be_vert[vi]];
2372  }
2373  }
2374  }
2375 
2376  // Destroy tables that need to be rebuild
2377  DeleteTables();
2378 
2379  if (Dim > 1)
2380  {
2381  // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
2382  el_to_edge = new Table;
2383  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2384  }
2385  if (Dim > 2)
2386  {
2387  // generate el_to_face, be_to_face
2388  GetElementToFaceTable();
2389  }
2390  // Update faces and faces_info
2391  GenerateFaces();
2392 
2393  // Build the nodes from the saved locations if they were around before
2394  if (Nodes)
2395  {
2396  // To force FE space update, we need to increase 'sequence':
2397  sequence++;
2398  last_operation = Mesh::NONE;
2399  nodes_fes->Update(false); // want_transform = false
2400  Nodes->Update(); // just needed to update Nodes->sequence
2401  Array<int> new_dofs;
2402  for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
2403  {
2404  int new_elid = ordering[old_elid];
2405  nodes_fes->GetElementVDofs(new_elid, new_dofs);
2406  Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2407  delete old_elem_node_vals[old_elid];
2408  }
2409  }
2410 }
2411 
2412 
2413 void Mesh::MarkForRefinement()
2414 {
2415  if (meshgen & 1)
2416  {
2417  if (Dim == 2)
2418  {
2419  MarkTriMeshForRefinement();
2420  }
2421  else if (Dim == 3)
2422  {
2423  DSTable v_to_v(NumOfVertices);
2424  GetVertexToVertexTable(v_to_v);
2425  MarkTetMeshForRefinement(v_to_v);
2426  }
2427  }
2428 }
2429 
2430 void Mesh::MarkTriMeshForRefinement()
2431 {
2432  // Mark the longest triangle edge by rotating the indices so that
2433  // vertex 0 - vertex 1 is the longest edge in the triangle.
2434  DenseMatrix pmat;
2435  for (int i = 0; i < NumOfElements; i++)
2436  {
2437  if (elements[i]->GetType() == Element::TRIANGLE)
2438  {
2439  GetPointMatrix(i, pmat);
2440  static_cast<Triangle*>(elements[i])->MarkEdge(pmat);
2441  }
2442  }
2443 }
2444 
2445 void Mesh::GetEdgeOrdering(DSTable &v_to_v, Array<int> &order)
2446 {
2447  NumOfEdges = v_to_v.NumberOfEntries();
2448  order.SetSize(NumOfEdges);
2449  Array<Pair<double, int> > length_idx(NumOfEdges);
2450 
2451  for (int i = 0; i < NumOfVertices; i++)
2452  {
2453  for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
2454  {
2455  int j = it.Index();
2456  length_idx[j].one = GetLength(i, it.Column());
2457  length_idx[j].two = j;
2458  }
2459  }
2460 
2461  // Sort by increasing edge-length.
2462  length_idx.Sort();
2463 
2464  for (int i = 0; i < NumOfEdges; i++)
2465  {
2466  order[length_idx[i].two] = i;
2467  }
2468 }
2469 
2470 void Mesh::MarkTetMeshForRefinement(DSTable &v_to_v)
2471 {
2472  // Mark the longest tetrahedral edge by rotating the indices so that
2473  // vertex 0 - vertex 1 is the longest edge in the element.
2474  Array<int> order;
2475  GetEdgeOrdering(v_to_v, order);
2476 
2477  for (int i = 0; i < NumOfElements; i++)
2478  {
2479  if (elements[i]->GetType() == Element::TETRAHEDRON)
2480  {
2481  elements[i]->MarkEdge(v_to_v, order);
2482  }
2483  }
2484  for (int i = 0; i < NumOfBdrElements; i++)
2485  {
2486  if (boundary[i]->GetType() == Element::TRIANGLE)
2487  {
2488  boundary[i]->MarkEdge(v_to_v, order);
2489  }
2490  }
2491 }
2492 
2493 void Mesh::PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
2494 {
2495  if (*old_v_to_v && *old_elem_vert)
2496  {
2497  return;
2498  }
2499 
2500  FiniteElementSpace *fes = Nodes->FESpace();
2501 
2502  if (*old_v_to_v == NULL)
2503  {
2504  bool need_v_to_v = false;
2505  Array<int> dofs;
2506  for (int i = 0; i < GetNEdges(); i++)
2507  {
2508  // Since edge indices may change, we need to permute edge interior dofs
2509  // any time an edge index changes and there is at least one dof on that
2510  // edge.
2511  fes->GetEdgeInteriorDofs(i, dofs);
2512  if (dofs.Size() > 0)
2513  {
2514  need_v_to_v = true;
2515  break;
2516  }
2517  }
2518  if (need_v_to_v)
2519  {
2520  *old_v_to_v = new DSTable(NumOfVertices);
2521  GetVertexToVertexTable(*(*old_v_to_v));
2522  }
2523  }
2524  if (*old_elem_vert == NULL)
2525  {
2526  bool need_elem_vert = false;
2527  Array<int> dofs;
2528  for (int i = 0; i < GetNE(); i++)
2529  {
2530  // Since element indices do not change, we need to permute element
2531  // interior dofs only when there are at least 2 interior dofs in an
2532  // element (assuming the nodal dofs are non-directional).
2533  fes->GetElementInteriorDofs(i, dofs);
2534  if (dofs.Size() > 1)
2535  {
2536  need_elem_vert = true;
2537  break;
2538  }
2539  }
2540  if (need_elem_vert)
2541  {
2542  *old_elem_vert = new Table;
2543  (*old_elem_vert)->MakeI(GetNE());
2544  for (int i = 0; i < GetNE(); i++)
2545  {
2546  (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2547  }
2548  (*old_elem_vert)->MakeJ();
2549  for (int i = 0; i < GetNE(); i++)
2550  {
2551  (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2552  elements[i]->GetNVertices());
2553  }
2554  (*old_elem_vert)->ShiftUpI();
2555  }
2556  }
2557 }
2558 
2559 void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
2560 {
2561  FiniteElementSpace *fes = Nodes->FESpace();
2562  const FiniteElementCollection *fec = fes->FEColl();
2563  Array<int> old_dofs, new_dofs;
2564 
2565  // assuming that all edges have the same number of dofs
2566  if (NumOfEdges) { fes->GetEdgeInteriorDofs(0, old_dofs); }
2567  const int num_edge_dofs = old_dofs.Size();
2568 
2569  // Save the original nodes
2570  const Vector onodes = *Nodes;
2571 
2572  // vertex dofs do not need to be moved
2573  fes->GetVertexDofs(0, old_dofs);
2574  int offset = NumOfVertices * old_dofs.Size();
2575 
2576  // edge dofs:
2577  // edge enumeration may be different but edge orientation is the same
2578  if (num_edge_dofs > 0)
2579  {
2580  DSTable new_v_to_v(NumOfVertices);
2581  GetVertexToVertexTable(new_v_to_v);
2582 
2583  for (int i = 0; i < NumOfVertices; i++)
2584  {
2585  for (DSTable::RowIterator it(new_v_to_v, i); !it; ++it)
2586  {
2587  const int old_i = (*old_v_to_v)(i, it.Column());
2588  const int new_i = it.Index();
2589  if (new_i == old_i) { continue; }
2590 
2591  old_dofs.SetSize(num_edge_dofs);
2592  new_dofs.SetSize(num_edge_dofs);
2593  for (int j = 0; j < num_edge_dofs; j++)
2594  {
2595  old_dofs[j] = offset + old_i * num_edge_dofs + j;
2596  new_dofs[j] = offset + new_i * num_edge_dofs + j;
2597  }
2598  fes->DofsToVDofs(old_dofs);
2599  fes->DofsToVDofs(new_dofs);
2600  for (int j = 0; j < old_dofs.Size(); j++)
2601  {
2602  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2603  }
2604  }
2605  }
2606  offset += NumOfEdges * num_edge_dofs;
2607  }
2608 
2609  // face dofs:
2610  // both enumeration and orientation of the faces may be different
2611  if (fes->GetNFDofs() > 0)
2612  {
2613  // generate the old face-vertex table using the unmodified 'faces'
2614  Table old_face_vertex;
2615  old_face_vertex.MakeI(NumOfFaces);
2616  for (int i = 0; i < NumOfFaces; i++)
2617  {
2618  old_face_vertex.AddColumnsInRow(i, faces[i]->GetNVertices());
2619  }
2620  old_face_vertex.MakeJ();
2621  for (int i = 0; i < NumOfFaces; i++)
2622  old_face_vertex.AddConnections(i, faces[i]->GetVertices(),
2623  faces[i]->GetNVertices());
2624  old_face_vertex.ShiftUpI();
2625 
2626  // update 'el_to_face', 'be_to_face', 'faces', and 'faces_info'
2627  STable3D *faces_tbl = GetElementToFaceTable(1);
2628  GenerateFaces();
2629 
2630  // compute the new face dof offsets
2631  Array<int> new_fdofs(NumOfFaces+1);
2632  new_fdofs[0] = 0;
2633  for (int i = 0; i < NumOfFaces; i++) // i = old face index
2634  {
2635  const int *old_v = old_face_vertex.GetRow(i);
2636  int new_i; // new face index
2637  switch (old_face_vertex.RowSize(i))
2638  {
2639  case 3:
2640  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2641  break;
2642  case 4:
2643  default:
2644  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2645  break;
2646  }
2647  fes->GetFaceInteriorDofs(i, old_dofs);
2648  new_fdofs[new_i+1] = old_dofs.Size();
2649  }
2650  new_fdofs.PartialSum();
2651 
2652  // loop over the old face numbers
2653  for (int i = 0; i < NumOfFaces; i++)
2654  {
2655  const int *old_v = old_face_vertex.GetRow(i), *new_v;
2656  const int *dof_ord;
2657  int new_i, new_or;
2658  switch (old_face_vertex.RowSize(i))
2659  {
2660  case 3:
2661  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2662  new_v = faces[new_i]->GetVertices();
2663  new_or = GetTriOrientation(old_v, new_v);
2664  dof_ord = fec->DofOrderForOrientation(Geometry::TRIANGLE, new_or);
2665  break;
2666  case 4:
2667  default:
2668  new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2669  new_v = faces[new_i]->GetVertices();
2670  new_or = GetQuadOrientation(old_v, new_v);
2671  dof_ord = fec->DofOrderForOrientation(Geometry::SQUARE, new_or);
2672  break;
2673  }
2674 
2675  fes->GetFaceInteriorDofs(i, old_dofs);
2676  new_dofs.SetSize(old_dofs.Size());
2677  for (int j = 0; j < old_dofs.Size(); j++)
2678  {
2679  // we assume the dofs are non-directional, i.e. dof_ord[j] is >= 0
2680  const int old_j = dof_ord[j];
2681  new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2682  }
2683  fes->DofsToVDofs(old_dofs);
2684  fes->DofsToVDofs(new_dofs);
2685  for (int j = 0; j < old_dofs.Size(); j++)
2686  {
2687  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2688  }
2689  }
2690 
2691  offset += fes->GetNFDofs();
2692  delete faces_tbl;
2693  }
2694 
2695  // element dofs:
2696  // element orientation may be different
2697  if (old_elem_vert) // have elements with 2 or more dofs
2698  {
2699  // matters when the 'fec' is
2700  // (this code is executed only for triangles/tets)
2701  // - Pk on triangles, k >= 4
2702  // - Qk on quads, k >= 3
2703  // - Pk on tets, k >= 5
2704  // - Qk on hexes, k >= 3
2705  // - DG spaces
2706  // - ...
2707 
2708  // loop over all elements
2709  for (int i = 0; i < GetNE(); i++)
2710  {
2711  const int *old_v = old_elem_vert->GetRow(i);
2712  const int *new_v = elements[i]->GetVertices();
2713  const int *dof_ord;
2714  int new_or;
2715  const Geometry::Type geom = elements[i]->GetGeometryType();
2716  switch (geom)
2717  {
2718  case Geometry::SEGMENT:
2719  new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2720  break;
2721  case Geometry::TRIANGLE:
2722  new_or = GetTriOrientation(old_v, new_v);
2723  break;
2724  case Geometry::SQUARE:
2725  new_or = GetQuadOrientation(old_v, new_v);
2726  break;
2727  case Geometry::TETRAHEDRON:
2728  new_or = GetTetOrientation(old_v, new_v);
2729  break;
2730  default:
2731  new_or = 0;
2732  MFEM_ABORT(Geometry::Name[geom] << " elements (" << fec->Name()
2733  << " FE collection) are not supported yet!");
2734  break;
2735  }
2736  dof_ord = fec->DofOrderForOrientation(geom, new_or);
2737  MFEM_VERIFY(dof_ord != NULL,
2738  "FE collection '" << fec->Name()
2739  << "' does not define reordering for "
2740  << Geometry::Name[geom] << " elements!");
2741  fes->GetElementInteriorDofs(i, old_dofs);
2742  new_dofs.SetSize(old_dofs.Size());
2743  for (int j = 0; j < new_dofs.Size(); j++)
2744  {
2745  // we assume the dofs are non-directional, i.e. dof_ord[j] is >= 0
2746  const int old_j = dof_ord[j];
2747  new_dofs[old_j] = offset + j;
2748  }
2749  offset += new_dofs.Size();
2750  fes->DofsToVDofs(old_dofs);
2751  fes->DofsToVDofs(new_dofs);
2752  for (int j = 0; j < old_dofs.Size(); j++)
2753  {
2754  (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2755  }
2756  }
2757  }
2758 
2759  // Update Tables, faces, etc
2760  if (Dim > 2)
2761  {
2762  if (fes->GetNFDofs() == 0)
2763  {
2764  // needed for FE spaces that have face dofs, even if
2765  // the 'Nodes' do not have face dofs.
2766  GetElementToFaceTable();
2767  GenerateFaces();
2768  }
2769  CheckBdrElementOrientation();
2770  }
2771  if (el_to_edge)
2772  {
2773  // update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
2774  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2775  if (Dim == 2)
2776  {
2777  // update 'faces' and 'faces_info'
2778  GenerateFaces();
2779  CheckBdrElementOrientation();
2780  }
2781  }
2782  // To force FE space update, we need to increase 'sequence':
2783  sequence++;
2784  last_operation = Mesh::NONE;
2785  fes->Update(false); // want_transform = false
2786  Nodes->Update(); // just needed to update Nodes->sequence
2787 }
2788 
2789 void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
2790 {
2791  FinalizeCheck();
2792  CheckElementOrientation(fix_orientation);
2793 
2794  if (NumOfBdrElements == 0)
2795  {
2796  GetElementToFaceTable();
2797  GenerateFaces();
2798  GenerateBoundaryElements();
2799  }
2800 
2801  if (refine)
2802  {
2803  DSTable v_to_v(NumOfVertices);
2804  GetVertexToVertexTable(v_to_v);
2805  MarkTetMeshForRefinement(v_to_v);
2806  }
2807 
2808  GetElementToFaceTable();
2809  GenerateFaces();
2810 
2811  CheckBdrElementOrientation();
2812 
2813  if (generate_edges == 1)
2814  {
2815  el_to_edge = new Table;
2816  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2817  }
2818  else
2819  {
2820  el_to_edge = NULL; // Not really necessary -- InitTables was called
2821  bel_to_edge = NULL;
2822  NumOfEdges = 0;
2823  }
2824 
2825  SetAttributes();
2826 
2827  SetMeshGen();
2828 }
2829 
2830 void Mesh::FinalizeWedgeMesh(int generate_edges, int refine,
2831  bool fix_orientation)
2832 {
2833  FinalizeCheck();
2834  CheckElementOrientation(fix_orientation);
2835 
2836  if (NumOfBdrElements == 0)
2837  {
2838  GetElementToFaceTable();
2839  GenerateFaces();
2840  GenerateBoundaryElements();
2841  }
2842 
2843  GetElementToFaceTable();
2844  GenerateFaces();
2845 
2846  CheckBdrElementOrientation();
2847 
2848  if (generate_edges == 1)
2849  {
2850  el_to_edge = new Table;
2851  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2852  }
2853  else
2854  {
2855  el_to_edge = NULL; // Not really necessary -- InitTables was called
2856  bel_to_edge = NULL;
2857  NumOfEdges = 0;
2858  }
2859 
2860  SetAttributes();
2861 
2862  SetMeshGen();
2863 }
2864 
2865 void Mesh::FinalizeHexMesh(int generate_edges, int refine, bool fix_orientation)
2866 {
2867  FinalizeCheck();
2868  CheckElementOrientation(fix_orientation);
2869 
2870  GetElementToFaceTable();
2871  GenerateFaces();
2872 
2873  if (NumOfBdrElements == 0)
2874  {
2875  GenerateBoundaryElements();
2876  }
2877 
2878  CheckBdrElementOrientation();
2879 
2880  if (generate_edges)
2881  {
2882  el_to_edge = new Table;
2883  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2884  }
2885  else
2886  {
2887  NumOfEdges = 0;
2888  }
2889 
2890  SetAttributes();
2891 
2892  SetMeshGen();
2893 }
2894 
2895 void Mesh::FinalizeMesh(int refine, bool fix_orientation)
2896 {
2897  FinalizeTopology();
2898 
2899  Finalize(refine, fix_orientation);
2900 }
2901 
2902 void Mesh::FinalizeTopology(bool generate_bdr)
2903 {
2904  // Requirements: the following should be defined:
2905  // 1) Dim
2906  // 2) NumOfElements, elements
2907  // 3) NumOfBdrElements, boundary
2908  // 4) NumOfVertices
2909  // Optional:
2910  // 2) ncmesh may be defined
2911  // 3) el_to_edge may be allocated (it will be re-computed)
2912 
2913  FinalizeCheck();
2914  bool generate_edges = true;
2915 
2916  if (spaceDim == 0) { spaceDim = Dim; }
2917  if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2918 
2919  // if the user defined any hanging nodes (see AddVertexParent),
2920  // we're initializing a non-conforming mesh
2921  if (tmp_vertex_parents.Size())
2922  {
2923  MFEM_VERIFY(ncmesh == NULL, "");
2924  ncmesh = new NCMesh(this);
2925 
2926  // we need to recreate the Mesh because NCMesh reorders the vertices
2927  // (see NCMesh::UpdateVertices())
2928  InitFromNCMesh(*ncmesh);
2929  ncmesh->OnMeshUpdated(this);
2930  GenerateNCFaceInfo();
2931 
2932  SetAttributes();
2933 
2934  tmp_vertex_parents.DeleteAll();
2935  return;
2936  }
2937 
2938  // set the mesh type: 'meshgen', ...
2939  SetMeshGen();
2940 
2941  // generate the faces
2942  if (Dim > 2)
2943  {
2944  GetElementToFaceTable();
2945  GenerateFaces();
2946  if (NumOfBdrElements == 0 && generate_bdr)
2947  {
2948  GenerateBoundaryElements();
2949  GetElementToFaceTable(); // update be_to_face
2950  }
2951  }
2952  else
2953  {
2954  NumOfFaces = 0;
2955  }
2956 
2957  // generate edges if requested
2958  if (Dim > 1 && generate_edges)
2959  {
2960  // el_to_edge may already be allocated (P2 VTK meshes)
2961  if (!el_to_edge) { el_to_edge = new Table; }
2962  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2963  if (Dim == 2)
2964  {
2965  GenerateFaces(); // 'Faces' in 2D refers to the edges
2966  if (NumOfBdrElements == 0 && generate_bdr)
2967  {
2968  GenerateBoundaryElements();
2969  }
2970  }
2971  }
2972  else
2973  {
2974  NumOfEdges = 0;
2975  }
2976 
2977  if (Dim == 1)
2978  {
2979  GenerateFaces();
2980  }
2981 
2982  if (ncmesh)
2983  {
2984  // tell NCMesh the numbering of edges/faces
2985  ncmesh->OnMeshUpdated(this);
2986 
2987  // update faces_info with NC relations
2988  GenerateNCFaceInfo();
2989  }
2990 
2991  // generate the arrays 'attributes' and 'bdr_attributes'
2992  SetAttributes();
2993 }
2994 
2995 void Mesh::Finalize(bool refine, bool fix_orientation)
2996 {
2997  if (NURBSext || ncmesh)
2998  {
2999  MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
3000  MFEM_ASSERT(CheckBdrElementOrientation() == 0, "");
3001  return;
3002  }
3003 
3004  // Requirements:
3005  // 1) FinalizeTopology() or equivalent was called
3006  // 2) if (Nodes == NULL), vertices must be defined
3007  // 3) if (Nodes != NULL), Nodes must be defined
3008 
3009  const bool check_orientation = true; // for regular elements, not boundary
3010  const bool curved = (Nodes != NULL);
3011  const bool may_change_topology =
3012  ( refine && (Dim > 1 && (meshgen & 1)) ) ||
3013  ( check_orientation && fix_orientation &&
3014  (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
3015 
3016  DSTable *old_v_to_v = NULL;
3017  Table *old_elem_vert = NULL;
3018 
3019  if (curved && may_change_topology)
3020  {
3021  PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3022  }
3023 
3024  if (check_orientation)
3025  {
3026  // check and optionally fix element orientation
3027  CheckElementOrientation(fix_orientation);
3028  }
3029  if (refine)
3030  {
3031  MarkForRefinement(); // may change topology!
3032  }
3033 
3034  if (may_change_topology)
3035  {
3036  if (curved)
3037  {
3038  DoNodeReorder(old_v_to_v, old_elem_vert); // updates the mesh topology
3039  delete old_elem_vert;
3040  delete old_v_to_v;
3041  }
3042  else
3043  {
3044  FinalizeTopology(); // Re-computes some data unnecessarily.
3045  }
3046 
3047  // TODO: maybe introduce Mesh::NODE_REORDER operation and FESpace::
3048  // NodeReorderMatrix and do Nodes->Update() instead of DoNodeReorder?
3049  }
3050 
3051  // check and fix boundary element orientation
3052  CheckBdrElementOrientation();
3053 
3054 #ifdef MFEM_DEBUG
3055  // For non-orientable surfaces/manifolds, the check below will fail, so we
3056  // only perform it when Dim == spaceDim.
3057  if (Dim >= 2 && Dim == spaceDim)
3058  {
3059  const int num_faces = GetNumFaces();
3060  for (int i = 0; i < num_faces; i++)
3061  {
3062  MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
3063  faces_info[i].Elem2Inf%2 != 0, "Invalid mesh topology."
3064  " Interior face with incompatible orientations.");
3065  }
3066  }
3067 #endif
3068 }
3069 
3070 void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
3071  double sx, double sy, double sz, bool sfc_ordering)
3072 {
3073  int x, y, z;
3074 
3075  int NVert, NElem, NBdrElem;
3076 
3077  NVert = (nx+1) * (ny+1) * (nz+1);
3078  NElem = nx * ny * nz;
3079  NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3080  if (type == Element::TETRAHEDRON)
3081  {
3082  NElem *= 6;
3083  NBdrElem *= 2;
3084  }
3085  else if (type == Element::WEDGE)
3086  {
3087  NElem *= 2;
3088  NBdrElem += 2*nx*ny;
3089  }
3090  else if (type == Element::PYRAMID)
3091  {
3092  NElem *= 6;
3093  NVert += nx * ny * nz;
3094  }
3095 
3096  InitMesh(3, 3, NVert, NElem, NBdrElem);
3097 
3098  double coord[3];
3099  int ind[9];
3100 
3101  // Sets vertices and the corresponding coordinates
3102  for (z = 0; z <= nz; z++)
3103  {
3104  coord[2] = ((double) z / nz) * sz;
3105  for (y = 0; y <= ny; y++)
3106  {
3107  coord[1] = ((double) y / ny) * sy;
3108  for (x = 0; x <= nx; x++)
3109  {
3110  coord[0] = ((double) x / nx) * sx;
3111  AddVertex(coord);
3112  }
3113  }
3114  }
3115  if (type == Element::PYRAMID)
3116  {
3117  for (z = 0; z < nz; z++)
3118  {
3119  coord[2] = (((double) z + 0.5) / nz) * sz;
3120  for (y = 0; y < ny; y++)
3121  {
3122  coord[1] = (((double) y + 0.5 ) / ny) * sy;
3123  for (x = 0; x < nx; x++)
3124  {
3125  coord[0] = (((double) x + 0.5 ) / nx) * sx;
3126  AddVertex(coord);
3127  }
3128  }
3129  }
3130  }
3131 
3132 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3133 #define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3134 
3135  // Sets elements and the corresponding indices of vertices
3136  if (sfc_ordering && type == Element::HEXAHEDRON)
3137  {
3138  Array<int> sfc;
3139  NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
3140  MFEM_VERIFY(sfc.Size() == 3*nx*ny*nz, "");
3141 
3142  for (int k = 0; k < nx*ny*nz; k++)
3143  {
3144  x = sfc[3*k + 0];
3145  y = sfc[3*k + 1];
3146  z = sfc[3*k + 2];
3147 
3148  // *INDENT-OFF*
3149  ind[0] = VTX(x , y , z );
3150  ind[1] = VTX(x+1, y , z );
3151  ind[2] = VTX(x+1, y+1, z );
3152  ind[3] = VTX(x , y+1, z );
3153  ind[4] = VTX(x , y , z+1);
3154  ind[5] = VTX(x+1, y , z+1);
3155  ind[6] = VTX(x+1, y+1, z+1);
3156  ind[7] = VTX(x , y+1, z+1);
3157  // *INDENT-ON*
3158 
3159  AddHex(ind, 1);
3160  }
3161  }
3162  else
3163  {
3164  for (z = 0; z < nz; z++)
3165  {
3166  for (y = 0; y < ny; y++)
3167  {
3168  for (x = 0; x < nx; x++)
3169  {
3170  // *INDENT-OFF*
3171  ind[0] = VTX(x , y , z );
3172  ind[1] = VTX(x+1, y , z );
3173  ind[2] = VTX(x+1, y+1, z );
3174  ind[3] = VTX(x , y+1, z );
3175  ind[4] = VTX(x , y , z+1);
3176  ind[5] = VTX(x+1, y , z+1);
3177  ind[6] = VTX(x+1, y+1, z+1);
3178  ind[7] = VTX( x, y+1, z+1);
3179  // *INDENT-ON*
3180  if (type == Element::TETRAHEDRON)
3181  {
3182  AddHexAsTets(ind, 1);
3183  }
3184  else if (type == Element::WEDGE)
3185  {
3186  AddHexAsWedges(ind, 1);
3187  }
3188  else if (type == Element::PYRAMID)
3189  {
3190  ind[8] = VTXP(x, y, z);
3191  AddHexAsPyramids(ind, 1);
3192  }
3193  else
3194  {
3195  AddHex(ind, 1);
3196  }
3197  }
3198  }
3199  }
3200  }
3201 
3202  // Sets boundary elements and the corresponding indices of vertices
3203  // bottom, bdr. attribute 1
3204  for (y = 0; y < ny; y++)
3205  {
3206  for (x = 0; x < nx; x++)
3207  {
3208  // *INDENT-OFF*
3209  ind[0] = VTX(x , y , 0);
3210  ind[1] = VTX(x , y+1, 0);
3211  ind[2] = VTX(x+1, y+1, 0);
3212  ind[3] = VTX(x+1, y , 0);
3213  // *INDENT-ON*
3214  if (type == Element::TETRAHEDRON)
3215  {
3216  AddBdrQuadAsTriangles(ind, 1);
3217  }
3218  else if (type == Element::WEDGE)
3219  {
3220  AddBdrQuadAsTriangles(ind, 1);
3221  }
3222  else
3223  {
3224  AddBdrQuad(ind, 1);
3225  }
3226  }
3227  }
3228  // top, bdr. attribute 6
3229  for (y = 0; y < ny; y++)
3230  {
3231  for (x = 0; x < nx; x++)
3232  {
3233  // *INDENT-OFF*
3234  ind[0] = VTX(x , y , nz);
3235  ind[1] = VTX(x+1, y , nz);
3236  ind[2] = VTX(x+1, y+1, nz);
3237  ind[3] = VTX(x , y+1, nz);
3238  // *INDENT-ON*
3239  if (type == Element::TETRAHEDRON)
3240  {
3241  AddBdrQuadAsTriangles(ind, 6);
3242  }
3243  else if (type == Element::WEDGE)
3244  {
3245  AddBdrQuadAsTriangles(ind, 6);
3246  }
3247  else
3248  {
3249  AddBdrQuad(ind, 6);
3250  }
3251  }
3252  }
3253  // left, bdr. attribute 5
3254  for (z = 0; z < nz; z++)
3255  {
3256  for (y = 0; y < ny; y++)
3257  {
3258  // *INDENT-OFF*
3259  ind[0] = VTX(0 , y , z );
3260  ind[1] = VTX(0 , y , z+1);
3261  ind[2] = VTX(0 , y+1, z+1);
3262  ind[3] = VTX(0 , y+1, z );
3263  // *INDENT-ON*
3264  if (type == Element::TETRAHEDRON)
3265  {
3266  AddBdrQuadAsTriangles(ind, 5);
3267  }
3268  else
3269  {
3270  AddBdrQuad(ind, 5);
3271  }
3272  }
3273  }
3274  // right, bdr. attribute 3
3275  for (z = 0; z < nz; z++)
3276  {
3277  for (y = 0; y < ny; y++)
3278  {
3279  // *INDENT-OFF*
3280  ind[0] = VTX(nx, y , z );
3281  ind[1] = VTX(nx, y+1, z );
3282  ind[2] = VTX(nx, y+1, z+1);
3283  ind[3] = VTX(nx, y , z+1);
3284  // *INDENT-ON*
3285  if (type == Element::TETRAHEDRON)
3286  {
3287  AddBdrQuadAsTriangles(ind, 3);
3288  }
3289  else
3290  {
3291  AddBdrQuad(ind, 3);
3292  }
3293  }
3294  }
3295  // front, bdr. attribute 2
3296  for (x = 0; x < nx; x++)
3297  {
3298  for (z = 0; z < nz; z++)
3299  {
3300  // *INDENT-OFF*
3301  ind[0] = VTX(x , 0, z );
3302  ind[1] = VTX(x+1, 0, z );
3303  ind[2] = VTX(x+1, 0, z+1);
3304  ind[3] = VTX(x , 0, z+1);
3305  // *INDENT-ON*
3306  if (type == Element::TETRAHEDRON)
3307  {
3308  AddBdrQuadAsTriangles(ind, 2);
3309  }
3310  else
3311  {
3312  AddBdrQuad(ind, 2);
3313  }
3314  }
3315  }
3316  // back, bdr. attribute 4
3317  for (x = 0; x < nx; x++)
3318  {
3319  for (z = 0; z < nz; z++)
3320  {
3321  // *INDENT-OFF*
3322  ind[0] = VTX(x , ny, z );
3323  ind[1] = VTX(x , ny, z+1);
3324  ind[2] = VTX(x+1, ny, z+1);
3325  ind[3] = VTX(x+1, ny, z );
3326  // *INDENT-ON*
3327  if (type == Element::TETRAHEDRON)
3328  {
3329  AddBdrQuadAsTriangles(ind, 4);
3330  }
3331  else
3332  {
3333  AddBdrQuad(ind, 4);
3334  }
3335  }
3336  }
3337 
3338 #undef VTX
3339 
3340 #if 0
3341  ofstream test_stream("debug.mesh");
3342  Print(test_stream);
3343  test_stream.close();
3344 #endif
3345 
3346  FinalizeTopology();
3347 
3348  // Finalize(...) can be called after this method, if needed
3349 }
3350 
3351 void Mesh::Make2D(int nx, int ny, Element::Type type,
3352  double sx, double sy,
3353  bool generate_edges, bool sfc_ordering)
3354 {
3355  int i, j, k;
3356 
3357  SetEmpty();
3358 
3359  Dim = spaceDim = 2;
3360 
3361  // Creates quadrilateral mesh
3362  if (type == Element::QUADRILATERAL)
3363  {
3364  NumOfVertices = (nx+1) * (ny+1);
3365  NumOfElements = nx * ny;
3366  NumOfBdrElements = 2 * nx + 2 * ny;
3367 
3368  vertices.SetSize(NumOfVertices);
3369  elements.SetSize(NumOfElements);
3370  boundary.SetSize(NumOfBdrElements);
3371 
3372  double cx, cy;
3373  int ind[4];
3374 
3375  // Sets vertices and the corresponding coordinates
3376  k = 0;
3377  for (j = 0; j < ny+1; j++)
3378  {
3379  cy = ((double) j / ny) * sy;
3380  for (i = 0; i < nx+1; i++)
3381  {
3382  cx = ((double) i / nx) * sx;
3383  vertices[k](0) = cx;
3384  vertices[k](1) = cy;
3385  k++;
3386  }
3387  }
3388 
3389  // Sets elements and the corresponding indices of vertices
3390  if (sfc_ordering)
3391  {
3392  Array<int> sfc;
3393  NCMesh::GridSfcOrdering2D(nx, ny, sfc);
3394  MFEM_VERIFY(sfc.Size() == 2*nx*ny, "");
3395 
3396  for (k = 0; k < nx*ny; k++)
3397  {
3398  i = sfc[2*k + 0];
3399  j = sfc[2*k + 1];
3400  ind[0] = i + j*(nx+1);
3401  ind[1] = i + 1 +j*(nx+1);
3402  ind[2] = i + 1 + (j+1)*(nx+1);
3403  ind[3] = i + (j+1)*(nx+1);
3404  elements[k] = new Quadrilateral(ind);
3405  }
3406  }
3407  else
3408  {
3409  k = 0;
3410  for (j = 0; j < ny; j++)
3411  {
3412  for (i = 0; i < nx; i++)
3413  {
3414  ind[0] = i + j*(nx+1);
3415  ind[1] = i + 1 +j*(nx+1);
3416  ind[2] = i + 1 + (j+1)*(nx+1);
3417  ind[3] = i + (j+1)*(nx+1);
3418  elements[k] = new Quadrilateral(ind);
3419  k++;
3420  }
3421  }
3422  }
3423 
3424  // Sets boundary elements and the corresponding indices of vertices
3425  int m = (nx+1)*ny;
3426  for (i = 0; i < nx; i++)
3427  {
3428  boundary[i] = new Segment(i, i+1, 1);
3429  boundary[nx+i] = new Segment(m+i+1, m+i, 3);
3430  }
3431  m = nx+1;
3432  for (j = 0; j < ny; j++)
3433  {
3434  boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
3435  boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
3436  }
3437  }
3438  // Creates triangular mesh
3439  else if (type == Element::TRIANGLE)
3440  {
3441  NumOfVertices = (nx+1) * (ny+1);
3442  NumOfElements = 2 * nx * ny;
3443  NumOfBdrElements = 2 * nx + 2 * ny;
3444 
3445  vertices.SetSize(NumOfVertices);
3446  elements.SetSize(NumOfElements);
3447  boundary.SetSize(NumOfBdrElements);
3448 
3449  double cx, cy;
3450  int ind[3];
3451 
3452  // Sets vertices and the corresponding coordinates
3453  k = 0;
3454  for (j = 0; j < ny+1; j++)
3455  {
3456  cy = ((double) j / ny) * sy;
3457  for (i = 0; i < nx+1; i++)
3458  {
3459  cx = ((double) i / nx) * sx;
3460  vertices[k](0) = cx;
3461  vertices[k](1) = cy;
3462  k++;
3463  }
3464  }
3465 
3466  // Sets the elements and the corresponding indices of vertices
3467  k = 0;
3468  for (j = 0; j < ny; j++)
3469  {
3470  for (i = 0; i < nx; i++)
3471  {
3472  ind[0] = i + j*(nx+1);
3473  ind[1] = i + 1 + (j+1)*(nx+1);
3474  ind[2] = i + (j+1)*(nx+1);
3475  elements[k] = new Triangle(ind);
3476  k++;
3477  ind[1] = i + 1 + j*(nx+1);
3478  ind[2] = i + 1 + (j+1)*(nx+1);
3479  elements[k] = new Triangle(ind);
3480  k++;
3481  }
3482  }
3483 
3484  // Sets boundary elements and the corresponding indices of vertices
3485  int m = (nx+1)*ny;
3486  for (i = 0; i < nx; i++)
3487  {
3488  boundary[i] = new Segment(i, i+1, 1);
3489  boundary[nx+i] = new Segment(m+i+1, m+i, 3);
3490  }
3491  m = nx+1;
3492  for (j = 0; j < ny; j++)
3493  {
3494  boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
3495  boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
3496  }
3497 
3498  // MarkTriMeshForRefinement(); // done in Finalize(...)
3499  }
3500  else
3501  {
3502  MFEM_ABORT("Unsupported element type.");
3503  }
3504 
3505  SetMeshGen();
3506  CheckElementOrientation();
3507 
3508  if (generate_edges == 1)
3509  {
3510  el_to_edge = new Table;
3511  NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3512  GenerateFaces();
3513  CheckBdrElementOrientation();
3514  }
3515  else
3516  {
3517  NumOfEdges = 0;
3518  }
3519 
3520  NumOfFaces = 0;
3521 
3522  attributes.Append(1);
3523  bdr_attributes.Append(1); bdr_attributes.Append(2);
3524  bdr_attributes.Append(3); bdr_attributes.Append(4);
3525 
3526  // Finalize(...) can be called after this method, if needed
3527 }
3528 
3529 void Mesh::Make1D(int n, double sx)
3530 {
3531  int j, ind[1];
3532 
3533  SetEmpty();
3534 
3535  Dim = 1;
3536  spaceDim = 1;
3537 
3538  NumOfVertices = n + 1;
3539  NumOfElements = n;
3540  NumOfBdrElements = 2;
3541  vertices.SetSize(NumOfVertices);
3542  elements.SetSize(NumOfElements);
3543  boundary.SetSize(NumOfBdrElements);
3544 
3545  // Sets vertices and the corresponding coordinates
3546  for (j = 0; j < n+1; j++)
3547  {
3548  vertices[j](0) = ((double) j / n) * sx;
3549  }
3550 
3551  // Sets elements and the corresponding indices of vertices
3552  for (j = 0; j < n; j++)
3553  {
3554  elements[j] = new Segment(j, j+1, 1);
3555  }
3556 
3557  // Sets the boundary elements
3558  ind[0] = 0;
3559  boundary[0] = new Point(ind, 1);
3560  ind[0] = n;
3561  boundary[1] = new Point(ind, 2);
3562 
3563  NumOfEdges = 0;
3564  NumOfFaces = 0;
3565 
3566  SetMeshGen();
3567  GenerateFaces();
3568 
3569  attributes.Append(1);
3570  bdr_attributes.Append(1); bdr_attributes.Append(2);
3571 }
3572 
3573 Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
3574 {
3575  Dim = mesh.Dim;
3576  spaceDim = mesh.spaceDim;
3577 
3578  NumOfVertices = mesh.NumOfVertices;
3579  NumOfElements = mesh.NumOfElements;
3580  NumOfBdrElements = mesh.NumOfBdrElements;
3581  NumOfEdges = mesh.NumOfEdges;
3582  NumOfFaces = mesh.NumOfFaces;
3583  nbInteriorFaces = mesh.nbInteriorFaces;
3584  nbBoundaryFaces = mesh.nbBoundaryFaces;
3585 
3586  meshgen = mesh.meshgen;
3587  mesh_geoms = mesh.mesh_geoms;
3588 
3589  // Create the new Mesh instance without a record of its refinement history
3590  sequence = 0;
3591  last_operation = Mesh::NONE;
3592 
3593  // Duplicate the elements
3594  elements.SetSize(NumOfElements);
3595  for (int i = 0; i < NumOfElements; i++)
3596  {
3597  elements[i] = mesh.elements[i]->Duplicate(this);
3598  }
3599 
3600  // Copy the vertices
3601  mesh.vertices.Copy(vertices);
3602 
3603  // Duplicate the boundary
3604  boundary.SetSize(NumOfBdrElements);
3605  for (int i = 0; i < NumOfBdrElements; i++)
3606  {
3607  boundary[i] = mesh.boundary[i]->Duplicate(this);
3608  }
3609 
3610  // Copy the element-to-face Table, el_to_face
3611  el_to_face = (mesh.el_to_face) ? new Table(*mesh.el_to_face) : NULL;
3612 
3613  // Copy the boundary-to-face Array, be_to_face.
3614  mesh.be_to_face.Copy(be_to_face);
3615 
3616  // Copy the element-to-edge Table, el_to_edge
3617  el_to_edge = (mesh.el_to_edge) ? new Table(*mesh.el_to_edge) : NULL;
3618 
3619  // Copy the boundary-to-edge Table, bel_to_edge (3D)
3620  bel_to_edge = (mesh.bel_to_edge) ? new Table(*mesh.bel_to_edge) : NULL;
3621 
3622  // Copy the boundary-to-edge Array, be_to_edge (2D)
3623  mesh.be_to_edge.Copy(be_to_edge);
3624 
3625  // Duplicate the faces and faces_info.
3626  faces.SetSize(mesh.faces.Size());
3627  for (int i = 0; i < faces.Size(); i++)
3628  {
3629  Element *face = mesh.faces[i]; // in 1D the faces are NULL
3630  faces[i] = (face) ? face->Duplicate(this) : NULL;
3631  }
3632  mesh.faces_info.Copy(faces_info);
3633  mesh.nc_faces_info.Copy(nc_faces_info);
3634 
3635  // Do NOT copy the element-to-element Table, el_to_el
3636  el_to_el = NULL;
3637 
3638  // Do NOT copy the face-to-edge Table, face_edge
3639  face_edge = NULL;
3640  face_to_elem = NULL;
3641 
3642  // Copy the edge-to-vertex Table, edge_vertex
3643  edge_vertex = (mesh.edge_vertex) ? new Table(*mesh.edge_vertex) : NULL;
3644 
3645  // Copy the attributes and bdr_attributes
3646  mesh.attributes.Copy(attributes);
3647  mesh.bdr_attributes.Copy(bdr_attributes);
3648 
3649  // Deep copy the NURBSExtension.
3650 #ifdef MFEM_USE_MPI
3651  ParNURBSExtension *pNURBSext =
3652  dynamic_cast<ParNURBSExtension *>(mesh.NURBSext);
3653  if (pNURBSext)
3654  {
3655  NURBSext = new ParNURBSExtension(*pNURBSext);
3656  }
3657  else
3658 #endif
3659  {
3660  NURBSext = mesh.NURBSext ? new NURBSExtension(*mesh.NURBSext) : NULL;
3661  }
3662 
3663  // Deep copy the NCMesh.
3664 #ifdef MFEM_USE_MPI
3665  if (dynamic_cast<const ParMesh*>(&mesh))
3666  {
3667  ncmesh = NULL; // skip; will be done in ParMesh copy ctor
3668  }
3669  else
3670 #endif
3671  {
3672  ncmesh = mesh.ncmesh ? new NCMesh(*mesh.ncmesh) : NULL;
3673  }
3674 
3675  // Duplicate the Nodes, including the FiniteElementCollection and the
3676  // FiniteElementSpace
3677  if (mesh.Nodes && copy_nodes)
3678  {
3679  FiniteElementSpace *fes = mesh.Nodes->FESpace();
3680  const FiniteElementCollection *fec = fes->FEColl();
3681  FiniteElementCollection *fec_copy =
3682  FiniteElementCollection::New(fec->Name());
3683  FiniteElementSpace *fes_copy =
3684  new FiniteElementSpace(*fes, this, fec_copy);
3685  Nodes = new GridFunction(fes_copy);
3686  Nodes->MakeOwner(fec_copy);
3687  *Nodes = *mesh.Nodes;
3688  own_nodes = 1;
3689  }
3690  else
3691  {
3692  Nodes = mesh.Nodes;
3693  own_nodes = 0;
3694  }
3695 }
3696 
3697 Mesh::Mesh(Mesh &&mesh) : Mesh()
3698 {
3699  Swap(mesh, true);
3700 }
3701 
3703 {
3704  Swap(mesh, true);
3705  return *this;
3706 }
3707 
3708 Mesh Mesh::LoadFromFile(const char *filename, int generate_edges, int refine,
3709  bool fix_orientation)
3710 {
3711  Mesh mesh;
3712  named_ifgzstream imesh(filename);
3713  if (!imesh) { MFEM_ABORT("Mesh file not found: " << filename << '\n'); }
3714  else { mesh.Load(imesh, generate_edges, refine, fix_orientation); }
3715  return mesh;
3716 }
3717 
3718 Mesh Mesh::MakeCartesian1D(int n, double sx)
3719 {
3720  Mesh mesh;
3721  mesh.Make1D(n, sx);
3722  // mesh.Finalize(); not needed in this case
3723  return mesh;
3724 }
3725 
3727  int nx, int ny, Element::Type type, bool generate_edges,
3728  double sx, double sy, bool sfc_ordering)
3729 {
3730  Mesh mesh;
3731  mesh.Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3732  mesh.Finalize(true); // refine = true
3733  return mesh;
3734 }
3735 
3737  int nx, int ny, int nz, Element::Type type,
3738  double sx, double sy, double sz, bool sfc_ordering)
3739 {
3740  Mesh mesh;
3741  mesh.Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3742  mesh.Finalize(true); // refine = true
3743  return mesh;
3744 }
3745 
3746 Mesh Mesh::MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
3747 {
3748  Mesh mesh;
3749  Array<int> ref_factors(orig_mesh.GetNE());
3750  ref_factors = ref_factor;
3751  mesh.MakeRefined_(orig_mesh, ref_factors, ref_type);
3752  return mesh;
3753 }
3754 
3755 Mesh Mesh::MakeRefined(Mesh &orig_mesh, const Array<int> &ref_factors,
3756  int ref_type)
3757 {
3758  Mesh mesh;
3759  mesh.MakeRefined_(orig_mesh, ref_factors, ref_type);
3760  return mesh;
3761 }
3762 
3763 Mesh::Mesh(const char *filename, int generate_edges, int refine,
3764  bool fix_orientation)
3765 {
3766  // Initialization as in the default constructor
3767  SetEmpty();
3768 
3769  named_ifgzstream imesh(filename);
3770  if (!imesh)
3771  {
3772  // Abort with an error message.
3773  MFEM_ABORT("Mesh file not found: " << filename << '\n');
3774  }
3775  else
3776  {
3777  Load(imesh, generate_edges, refine, fix_orientation);
3778  }
3779 }
3780 
3781 Mesh::Mesh(std::istream &input, int generate_edges, int refine,
3782  bool fix_orientation)
3783 {
3784  SetEmpty();
3785  Load(input, generate_edges, refine, fix_orientation);
3786 }
3787 
3788 void Mesh::ChangeVertexDataOwnership(double *vertex_data, int len_vertex_data,
3789  bool zerocopy)
3790 {
3791  // A dimension of 3 is now required since we use mfem::Vertex objects as PODs
3792  // and these object have a hardcoded double[3] entry
3793  MFEM_VERIFY(len_vertex_data >= NumOfVertices * 3,
3794  "Not enough vertices in external array : "
3795  "len_vertex_data = "<< len_vertex_data << ", "
3796  "NumOfVertices * 3 = " << NumOfVertices * 3);
3797  // Allow multiple calls to this method with the same vertex_data
3798  if (vertex_data == (double *)(vertices.GetData()))
3799  {
3800  MFEM_ASSERT(!vertices.OwnsData(), "invalid ownership");
3801  return;
3802  }
3803  if (!zerocopy)
3804  {
3805  memcpy(vertex_data, vertices.GetData(),
3806  NumOfVertices * 3 * sizeof(double));
3807  }
3808  // Vertex is POD double[3]
3809  vertices.MakeRef(reinterpret_cast<Vertex*>(vertex_data), NumOfVertices);
3810 }
3811 
3812 Mesh::Mesh(double *vertices_, int num_vertices,
3813  int *element_indices, Geometry::Type element_type,
3814  int *element_attributes, int num_elements,
3815  int *boundary_indices, Geometry::Type boundary_type,
3816  int *boundary_attributes, int num_boundary_elements,
3817  int dimension, int space_dimension)
3818 {
3819  if (space_dimension == -1)
3820  {
3821  space_dimension = dimension;
3822  }
3823 
3824  InitMesh(dimension, space_dimension, /*num_vertices*/ 0, num_elements,
3825  num_boundary_elements);
3826 
3827  int element_index_stride = Geometry::NumVerts[element_type];
3828  int boundary_index_stride = num_boundary_elements > 0 ?
3829  Geometry::NumVerts[boundary_type] : 0;
3830 
3831  // assuming Vertex is POD
3832  vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3833  NumOfVertices = num_vertices;
3834 
3835  for (int i = 0; i < num_elements; i++)
3836  {
3837  elements[i] = NewElement(element_type);
3838  elements[i]->SetVertices(element_indices + i * element_index_stride);
3839  elements[i]->SetAttribute(element_attributes[i]);
3840  }
3841  NumOfElements = num_elements;
3842 
3843  for (int i = 0; i < num_boundary_elements; i++)
3844  {
3845  boundary[i] = NewElement(boundary_type);
3846  boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3847  boundary[i]->SetAttribute(boundary_attributes[i]);
3848  }
3849  NumOfBdrElements = num_boundary_elements;
3850 
3851  FinalizeTopology();
3852 }
3853 
3855 {
3856  switch (geom)
3857  {
3858  case Geometry::POINT: return (new Point);
3859  case Geometry::SEGMENT: return (new Segment);
3860  case Geometry::TRIANGLE: return (new Triangle);
3861  case Geometry::SQUARE: return (new Quadrilateral);
3862  case Geometry::TETRAHEDRON:
3863 #ifdef MFEM_USE_MEMALLOC
3864  return TetMemory.Alloc();
3865 #else
3866  return (new Tetrahedron);
3867 #endif
3868  case Geometry::CUBE: return (new Hexahedron);
3869  case Geometry::PRISM: return (new Wedge);
3870  case Geometry::PYRAMID: return (new Pyramid);
3871  default:
3872  MFEM_ABORT("invalid Geometry::Type, geom = " << geom);
3873  }
3874 
3875  return NULL;
3876 }
3877 
3879 {
3880  int geom, nv, *v;
3881  Element *el;
3882 
3883  input >> geom;
3884  el = NewElement(geom);
3885  MFEM_VERIFY(el, "Unsupported element type: " << geom);
3886  nv = el->GetNVertices();
3887  v = el->GetVertices();
3888  for (int i = 0; i < nv; i++)
3889  {
3890  input >> v[i];
3891  }
3892 
3893  return el;
3894 }
3895 
3896 void Mesh::PrintElementWithoutAttr(const Element *el, std::ostream &os)
3897 {
3898  os << el->GetGeometryType();
3899  const int nv = el->GetNVertices();
3900  const int *v = el->GetVertices();
3901  for (int j = 0; j < nv; j++)
3902  {
3903  os << ' ' << v[j];
3904  }
3905  os << '\n';
3906 }
3907 
3908 Element *Mesh::ReadElement(std::istream &input)
3909 {
3910  int attr;
3911  Element *el;
3912 
3913  input >> attr;
3914  el = ReadElementWithoutAttr(input);
3915  el->SetAttribute(attr);
3916 
3917  return el;
3918 }
3919 
3920 void Mesh::PrintElement(const Element *el, std::ostream &os)
3921 {
3922  os << el->GetAttribute() << ' ';
3923  PrintElementWithoutAttr(el, os);
3924 }
3925 
3927 {
3928  meshgen = mesh_geoms = 0;
3929  for (int i = 0; i < NumOfElements; i++)
3930  {
3931  const Element::Type type = GetElement(i)->GetType();
3932  switch (type)
3933  {
3934  case Element::TETRAHEDRON:
3936  case Element::TRIANGLE:
3937  mesh_geoms |= (1 << Geometry::TRIANGLE);
3938  case Element::SEGMENT:
3939  mesh_geoms |= (1 << Geometry::SEGMENT);
3940  case Element::POINT:
3941  mesh_geoms |= (1 << Geometry::POINT);
3942  meshgen |= 1;
3943  break;
3944 
3945  case Element::HEXAHEDRON:
3946  mesh_geoms |= (1 << Geometry::CUBE);
3948  mesh_geoms |= (1 << Geometry::SQUARE);
3949  mesh_geoms |= (1 << Geometry::SEGMENT);
3950  mesh_geoms |= (1 << Geometry::POINT);
3951  meshgen |= 2;
3952  break;
3953 
3954  case Element::WEDGE:
3955  mesh_geoms |= (1 << Geometry::PRISM);
3956  mesh_geoms |= (1 << Geometry::SQUARE);
3957  mesh_geoms |= (1 << Geometry::TRIANGLE);
3958  mesh_geoms |= (1 << Geometry::SEGMENT);
3959  mesh_geoms |= (1 << Geometry::POINT);
3960  meshgen |= 4;
3961  break;
3962 
3963  case Element::PYRAMID:
3964  mesh_geoms |= (1 << Geometry::PYRAMID);
3965  mesh_geoms |= (1 << Geometry::SQUARE);
3966  mesh_geoms |= (1 << Geometry::TRIANGLE);
3967  mesh_geoms |= (1 << Geometry::SEGMENT);
3968  mesh_geoms |= (1 << Geometry::POINT);
3969  meshgen |= 8;
3970  break;
3971 
3972  default:
3973  MFEM_ABORT("invalid element type: " << type);
3974  break;
3975  }
3976  }
3977 }
3978 
3979 void Mesh::Loader(std::istream &input, int generate_edges,
3980  std::string parse_tag)
3981 {
3982  int curved = 0, read_gf = 1;
3983  bool finalize_topo = true;
3984 
3985  if (!input)
3986  {
3987  MFEM_ABORT("Input stream is not open");
3988  }
3989 
3990  Clear();
3991 
3992  string mesh_type;
3993  input >> ws;
3994  getline(input, mesh_type);
3995  filter_dos(mesh_type);
3996 
3997  // MFEM's conforming mesh formats
3998  int mfem_version = 0;
3999  if (mesh_type == "MFEM mesh v1.0") { mfem_version = 10; } // serial
4000  else if (mesh_type == "MFEM mesh v1.2") { mfem_version = 12; } // parallel
4001 
4002  // MFEM nonconforming mesh format
4003  // (NOTE: previous v1.1 is now under this branch for backward compatibility)
4004  int mfem_nc_version = 0;
4005  if (mesh_type == "MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4006  else if (mesh_type == "MFEM mesh v1.1") { mfem_nc_version = 1 /*legacy*/; }
4007 
4008  if (mfem_version)
4009  {
4010  // Formats mfem_v12 and newer have a tag indicating the end of the mesh
4011  // section in the stream. A user provided parse tag can also be provided
4012  // via the arguments. For example, if this is called from parallel mesh
4013  // object, it can indicate to read until parallel mesh section begins.
4014  if (mfem_version == 12 && parse_tag.empty())
4015  {
4016  parse_tag = "mfem_mesh_end";
4017  }
4018  ReadMFEMMesh(input, mfem_version, curved);
4019  }
4020  else if (mfem_nc_version)
4021  {
4022  MFEM_ASSERT(ncmesh == NULL, "internal error");
4023  int is_nc = 1;
4024 
4025 #ifdef MFEM_USE_MPI
4026  ParMesh *pmesh = dynamic_cast<ParMesh*>(this);
4027  if (pmesh)
4028  {
4029  MFEM_VERIFY(mfem_nc_version >= 10,
4030  "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4031  "used to load a parallel nonconforming mesh, sorry.");
4032 
4033  ncmesh = new ParNCMesh(pmesh->GetComm(),
4034  input, mfem_nc_version, curved, is_nc);
4035  }
4036  else
4037 #endif
4038  {
4039  ncmesh = new NCMesh(input, mfem_nc_version, curved, is_nc);
4040  }
4042 
4043  if (!is_nc)
4044  {
4045  // special case for backward compatibility with MFEM <=4.2:
4046  // if the "vertex_parents" section is missing in the v1.1 format,
4047  // the mesh is treated as conforming
4048  delete ncmesh;
4049  ncmesh = NULL;
4050  }
4051  }
4052  else if (mesh_type == "linemesh") // 1D mesh
4053  {
4054  ReadLineMesh(input);
4055  }
4056  else if (mesh_type == "areamesh2" || mesh_type == "curved_areamesh2")
4057  {
4058  if (mesh_type == "curved_areamesh2")
4059  {
4060  curved = 1;
4061  }
4062  ReadNetgen2DMesh(input, curved);
4063  }
4064  else if (mesh_type == "NETGEN" || mesh_type == "NETGEN_Neutral_Format")
4065  {
4066  ReadNetgen3DMesh(input);
4067  }
4068  else if (mesh_type == "TrueGrid")
4069  {
4070  ReadTrueGridMesh(input);
4071  }
4072  else if (mesh_type.rfind("# vtk DataFile Version") == 0)
4073  {
4074  int major_vtk_version = mesh_type[mesh_type.length()-3] - '0';
4075  // int minor_vtk_version = mesh_type[mesh_type.length()-1] - '0';
4076  MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4077  "Unsupported VTK format");
4078  ReadVTKMesh(input, curved, read_gf, finalize_topo);
4079  }
4080  else if (mesh_type.rfind("<VTKFile ") == 0 || mesh_type.rfind("<?xml") == 0)
4081  {
4082  ReadXML_VTKMesh(input, curved, read_gf, finalize_topo, mesh_type);
4083  }
4084  else if (mesh_type == "MFEM NURBS mesh v1.0")
4085  {
4086  ReadNURBSMesh(input, curved, read_gf);
4087  }
4088  else if (mesh_type == "MFEM INLINE mesh v1.0")
4089  {
4090  ReadInlineMesh(input, generate_edges);
4091  return; // done with inline mesh construction
4092  }
4093  else if (mesh_type == "$MeshFormat") // Gmsh
4094  {
4095  ReadGmshMesh(input, curved, read_gf);
4096  }
4097  else if
4098  ((mesh_type.size() > 2 &&
4099  mesh_type[0] == 'C' && mesh_type[1] == 'D' && mesh_type[2] == 'F') ||
4100  (mesh_type.size() > 3 &&
4101  mesh_type[1] == 'H' && mesh_type[2] == 'D' && mesh_type[3] == 'F'))
4102  {
4103  named_ifgzstream *mesh_input = dynamic_cast<named_ifgzstream *>(&input);
4104  if (mesh_input)
4105  {
4106 #ifdef MFEM_USE_NETCDF
4107  ReadCubit(mesh_input->filename.c_str(), curved, read_gf);
4108 #else
4109  MFEM_ABORT("NetCDF support requires configuration with"
4110  " MFEM_USE_NETCDF=YES");
4111  return;
4112 #endif
4113  }
4114  else
4115  {
4116  MFEM_ABORT("Can not determine Cubit mesh filename!"
4117  " Use mfem::named_ifgzstream for input.");
4118  return;
4119  }
4120  }
4121  else
4122  {
4123  MFEM_ABORT("Unknown input mesh format: " << mesh_type);
4124  return;
4125  }
4126 
4127  // at this point the following should be defined:
4128  // 1) Dim
4129  // 2) NumOfElements, elements
4130  // 3) NumOfBdrElements, boundary
4131  // 4) NumOfVertices, with allocated space in vertices
4132  // 5) curved
4133  // 5a) if curved == 0, vertices must be defined
4134  // 5b) if curved != 0 and read_gf != 0,
4135  // 'input' must point to a GridFunction
4136  // 5c) if curved != 0 and read_gf == 0,
4137  // vertices and Nodes must be defined
4138  // optional:
4139  // 1) el_to_edge may be allocated (as in the case of P2 VTK meshes)
4140  // 2) ncmesh may be allocated
4141 
4142  // FinalizeTopology() will:
4143  // - assume that generate_edges is true
4144  // - assume that refine is false
4145  // - does not check the orientation of regular and boundary elements
4146  if (finalize_topo)
4147  {
4148  // don't generate any boundary elements, especially in parallel
4149  bool generate_bdr = false;
4150 
4151  FinalizeTopology(generate_bdr);
4152  }
4153 
4154  if (curved && read_gf)
4155  {
4156  Nodes = new GridFunction(this, input);
4157 
4158  own_nodes = 1;
4159  spaceDim = Nodes->VectorDim();
4160  if (ncmesh) { ncmesh->spaceDim = spaceDim; }
4161 
4162  // Set vertex coordinates from the 'Nodes'
4164  }
4165 
4166  // If a parse tag was supplied, keep reading the stream until the tag is
4167  // encountered.
4168  if (mfem_version == 12)
4169  {
4170  string line;
4171  do
4172  {
4173  skip_comment_lines(input, '#');
4174  MFEM_VERIFY(input.good(), "Required mesh-end tag not found");
4175  getline(input, line);
4176  filter_dos(line);
4177  // mfem v1.2 may not have parse_tag in it, e.g. if trying to read a
4178  // serial mfem v1.2 mesh as parallel with "mfem_serial_mesh_end" as
4179  // parse_tag. That's why, regardless of parse_tag, we stop reading if
4180  // we find "mfem_mesh_end" which is required by mfem v1.2 format.
4181  if (line == "mfem_mesh_end") { break; }
4182  }
4183  while (line != parse_tag);
4184  }
4185  else if (mfem_nc_version >= 10)
4186  {
4187  string ident;
4188  skip_comment_lines(input, '#');
4189  input >> ident;
4190  MFEM_VERIFY(ident == "mfem_mesh_end",
4191  "invalid mesh: end of file tag not found");
4192  }
4193 
4194  // Finalize(...) should be called after this, if needed.
4195 }
4196 
4197 Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
4198 {
4199  int i, j, ie, ib, iv, *v, nv;
4200  Element *el;
4201  Mesh *m;
4202 
4203  SetEmpty();
4204 
4205  Dim = mesh_array[0]->Dimension();
4206  spaceDim = mesh_array[0]->SpaceDimension();
4207 
4208  if (mesh_array[0]->NURBSext)
4209  {
4210  // assuming the pieces form a partition of a NURBS mesh
4211  NURBSext = new NURBSExtension(mesh_array, num_pieces);
4212 
4215 
4217 
4218  // NumOfBdrElements = NURBSext->GetNBE();
4219  // NURBSext->GetBdrElementTopo(boundary);
4220 
4221  Array<int> lvert_vert, lelem_elem;
4222 
4223  // Here, for visualization purposes, we copy the boundary elements from
4224  // the individual pieces which include the interior boundaries. This
4225  // creates 'boundary' array that is different from the one generated by
4226  // the NURBSExtension which, in particular, makes the boundary-dof table
4227  // invalid. This, in turn, causes GetBdrElementTransformation to not
4228  // function properly.
4229  NumOfBdrElements = 0;
4230  for (i = 0; i < num_pieces; i++)
4231  {
4232  NumOfBdrElements += mesh_array[i]->GetNBE();
4233  }
4234  boundary.SetSize(NumOfBdrElements);
4235  vertices.SetSize(NumOfVertices);
4236  ib = 0;
4237  for (i = 0; i < num_pieces; i++)
4238  {
4239  m = mesh_array[i];
4240  m->NURBSext->GetVertexLocalToGlobal(lvert_vert);
4241  m->NURBSext->GetElementLocalToGlobal(lelem_elem);
4242  // copy the element attributes
4243  for (j = 0; j < m->GetNE(); j++)
4244  {
4245  elements[lelem_elem[j]]->SetAttribute(m->GetAttribute(j));
4246  }
4247  // copy the boundary
4248  for (j = 0; j < m->GetNBE(); j++)
4249  {
4250  el = m->GetBdrElement(j)->Duplicate(this);
4251  v = el->GetVertices();
4252  nv = el->GetNVertices();
4253  for (int k = 0; k < nv; k++)
4254  {
4255  v[k] = lvert_vert[v[k]];
4256  }
4257  boundary[ib++] = el;
4258  }
4259  // copy the vertices
4260  for (j = 0; j < m->GetNV(); j++)
4261  {
4262  vertices[lvert_vert[j]].SetCoords(m->SpaceDimension(),
4263  m->GetVertex(j));
4264  }
4265  }
4266  }
4267  else // not a NURBS mesh
4268  {
4269  NumOfElements = 0;
4270  NumOfBdrElements = 0;
4271  NumOfVertices = 0;
4272  for (i = 0; i < num_pieces; i++)
4273  {
4274  m = mesh_array[i];
4275  NumOfElements += m->GetNE();
4276  NumOfBdrElements += m->GetNBE();
4277  NumOfVertices += m->GetNV();
4278  }
4279  elements.SetSize(NumOfElements);
4280  boundary.SetSize(NumOfBdrElements);
4281  vertices.SetSize(NumOfVertices);
4282  ie = ib = iv = 0;
4283  for (i = 0; i < num_pieces; i++)
4284  {
4285  m = mesh_array[i];
4286  // copy the elements
4287  for (j = 0; j < m->GetNE(); j++)
4288  {
4289  el = m->GetElement(j)->Duplicate(this);
4290  v = el->GetVertices();
4291  nv = el->GetNVertices();
4292  for (int k = 0; k < nv; k++)
4293  {
4294  v[k] += iv;
4295  }
4296  elements[ie++] = el;
4297  }
4298  // copy the boundary elements
4299  for (j = 0; j < m->GetNBE(); j++)
4300  {
4301  el = m->GetBdrElement(j)->Duplicate(this);
4302  v = el->GetVertices();
4303  nv = el->GetNVertices();
4304  for (int k = 0; k < nv; k++)
4305  {
4306  v[k] += iv;
4307  }
4308  boundary[ib++] = el;
4309  }
4310  // copy the vertices
4311  for (j = 0; j < m->GetNV(); j++)
4312  {
4313  vertices[iv++].SetCoords(m->SpaceDimension(), m->GetVertex(j));
4314  }
4315  }
4316  }
4317 
4318  FinalizeTopology();
4319 
4320  // copy the nodes (curvilinear meshes)
4321  GridFunction *g = mesh_array[0]->GetNodes();
4322  if (g)
4323  {
4324  Array<GridFunction *> gf_array(num_pieces);
4325  for (i = 0; i < num_pieces; i++)
4326  {
4327  gf_array[i] = mesh_array[i]->GetNodes();
4328  }
4329  Nodes = new GridFunction(this, gf_array, num_pieces);
4330  own_nodes = 1;
4331  }
4332 
4333 #ifdef MFEM_DEBUG
4334  CheckElementOrientation(false);
4336 #endif
4337 }
4338 
4339 Mesh::Mesh(Mesh *orig_mesh, int ref_factor, int ref_type)
4340 {
4341  Array<int> ref_factors(orig_mesh->GetNE());
4342  ref_factors = ref_factor;
4343  MakeRefined_(*orig_mesh, ref_factors, ref_type);
4344 }
4345 
4346 void Mesh::MakeRefined_(Mesh &orig_mesh, const Array<int> ref_factors,
4347  int ref_type)
4348 {
4349  SetEmpty();
4350  Dim = orig_mesh.Dimension();
4351  spaceDim = orig_mesh.SpaceDimension();
4352 
4353  int orig_ne = orig_mesh.GetNE();
4354  MFEM_VERIFY(ref_factors.Size() == orig_ne && orig_ne > 0,
4355  "Number of refinement factors must equal number of elements")
4356  MFEM_VERIFY(ref_factors.Min() >= 1, "Refinement factor must be >= 1");
4357  const int q_type = BasisType::GetQuadrature1D(ref_type);
4358  MFEM_VERIFY(Quadrature1D::CheckClosed(q_type) != Quadrature1D::Invalid,
4359  "Invalid refinement type. Must use closed basis type.");
4360 
4361  int min_ref = ref_factors.Min();
4362  int max_ref = ref_factors.Max();
4363 
4364  bool var_order = (min_ref != max_ref);
4365 
4366  // variable order space can only be constructed over an NC mesh
4367  if (var_order) { orig_mesh.EnsureNCMesh(true); }
4368 
4369  // Construct a scalar H1 FE space of order ref_factor and use its dofs as
4370  // the indices of the new, refined vertices.
4371  H1_FECollection rfec(min_ref, Dim, ref_type);
4372  FiniteElementSpace rfes(&orig_mesh, &rfec);
4373 
4374  if (var_order)
4375  {
4376  rfes.SetRelaxedHpConformity(false);
4377  for (int i = 0; i < orig_ne; i++)
4378  {
4379  rfes.SetElementOrder(i, ref_factors[i]);
4380  }
4381  rfes.Update(false);
4382  }
4383 
4384  // Set the number of vertices, set the actual coordinates later
4385  NumOfVertices = rfes.GetNDofs();
4386  vertices.SetSize(NumOfVertices);
4387 
4388  Array<int> rdofs;
4389  DenseMatrix phys_pts;
4390 
4391  GeometryRefiner refiner;
4392  refiner.SetType(q_type);
4393 
4394  // Add refined elements and set vertex coordinates
4395  for (int el = 0; el < orig_ne; el++)
4396  {
4397  Geometry::Type geom = orig_mesh.GetElementGeometry(el);
4398  int attrib = orig_mesh.GetAttribute(el);
4399  int nvert = Geometry::NumVerts[geom];
4400  RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[el]);
4401 
4402  rfes.GetElementDofs(el, rdofs);
4403  MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
4404  const FiniteElement *rfe = rfes.GetFE(el);
4405  orig_mesh.GetElementTransformation(el)->Transform(rfe->GetNodes(),
4406  phys_pts);
4407  const int *c2h_map = rfec.GetDofMap(geom, ref_factors[el]);
4408  for (int i = 0; i < phys_pts.Width(); i++)
4409  {
4410  vertices[rdofs[i]].SetCoords(spaceDim, phys_pts.GetColumn(i));
4411  }
4412  for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
4413  {
4414  Element *elem = NewElement(geom);
4415  elem->SetAttribute(attrib);
4416  int *v = elem->GetVertices();
4417  for (int k = 0; k < nvert; k++)
4418  {
4419  int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
4420  v[k] = rdofs[c2h_map[cid]];
4421  }
4422  AddElement(elem);
4423  }
4424  }
4425 
4426  if (Dim > 2)
4427  {
4428  GetElementToFaceTable(false);
4429  GenerateFaces();
4430  }
4431 
4432  // Add refined boundary elements
4433  for (int el = 0; el < orig_mesh.GetNBE(); el++)
4434  {
4435  int i, info;
4436  orig_mesh.GetBdrElementAdjacentElement(el, i, info);
4437  Geometry::Type geom = orig_mesh.GetBdrElementBaseGeometry(el);
4438  int attrib = orig_mesh.GetBdrAttribute(el);
4439  int nvert = Geometry::NumVerts[geom];
4440  RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[i]);
4441 
4442  rfes.GetBdrElementDofs(el, rdofs);
4443  MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
4444  const int *c2h_map = rfec.GetDofMap(geom, ref_factors[i]);
4445  for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
4446  {
4447  Element *elem = NewElement(geom);
4448  elem->SetAttribute(attrib);
4449  int *v = elem->GetVertices();
4450  for (int k = 0; k < nvert; k++)
4451  {
4452  int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
4453  v[k] = rdofs[c2h_map[cid]];
4454  }
4455  AddBdrElement(elem);
4456  }
4457  }
4458  FinalizeTopology(false);
4459  sequence = orig_mesh.GetSequence() + 1;
4461 
4462  // Set up the nodes of the new mesh (if the original mesh has nodes). The new
4463  // mesh is always straight-sided (i.e. degree 1 finite element space), but
4464  // the nodes are required for e.g. periodic meshes.
4465  if (orig_mesh.GetNodes())
4466  {
4467  bool discont = orig_mesh.GetNodalFESpace()->IsDGSpace();
4468  Ordering::Type dof_ordering = orig_mesh.GetNodalFESpace()->GetOrdering();
4469  Mesh::SetCurvature(1, discont, spaceDim, dof_ordering);
4470  FiniteElementSpace *nodal_fes = Nodes->FESpace();
4471  const FiniteElementCollection *nodal_fec = nodal_fes->FEColl();
4472  H1_FECollection vertex_fec(1, Dim);
4473  Array<int> dofs;
4474  int el_counter = 0;
4475  for (int iel = 0; iel < orig_ne; iel++)
4476  {
4477  Geometry::Type geom = orig_mesh.GetElementBaseGeometry(iel);
4478  int nvert = Geometry::NumVerts[geom];
4479  RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[iel]);
4480  rfes.GetElementDofs(iel, rdofs);
4481  const FiniteElement *rfe = rfes.GetFE(iel);
4482  orig_mesh.GetElementTransformation(iel)->Transform(rfe->GetNodes(),
4483  phys_pts);
4484  const int *node_map = NULL;
4485  const H1_FECollection *h1_fec =
4486  dynamic_cast<const H1_FECollection *>(nodal_fec);
4487  if (h1_fec != NULL) { node_map = h1_fec->GetDofMap(geom); }
4488  const int *vertex_map = vertex_fec.GetDofMap(geom);
4489  const int *c2h_map = rfec.GetDofMap(geom, ref_factors[iel]);
4490  for (int jel = 0; jel < RG.RefGeoms.Size()/nvert; jel++)
4491  {
4492  nodal_fes->GetElementVDofs(el_counter++, dofs);
4493  for (int iv_lex=0; iv_lex<nvert; ++iv_lex)
4494  {
4495  // convert from lexicographic to vertex index
4496  int iv = vertex_map[iv_lex];
4497  // index of vertex of current element in phys_pts matrix
4498  int pt_idx = c2h_map[RG.RefGeoms[iv+nvert*jel]];
4499  // index of current vertex into DOF array
4500  int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4501  for (int d=0; d<spaceDim; ++d)
4502  {
4503  (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4504  }
4505  }
4506  }
4507  }
4508  }
4509 
4510  // Setup the data for the coarse-fine refinement transformations
4511  CoarseFineTr.embeddings.SetSize(GetNE());
4512  // First, compute total number of point matrices that we need per geometry
4513  // and the offsets into that array
4514  using GeomRef = std::pair<Geometry::Type, int>;
4515  std::map<GeomRef, int> point_matrices_offsets;
4516  int n_point_matrices[Geometry::NumGeom] = {}; // initialize to zero
4517  for (int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4518  {
4519  Geometry::Type geom = orig_mesh.GetElementBaseGeometry(el_coarse);
4520  // Have we seen this pair of (goemetry, refinement level) before?
4521  GeomRef id(geom, ref_factors[el_coarse]);
4522  if (point_matrices_offsets.find(id) == point_matrices_offsets.end())
4523  {
4524  RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[el_coarse]);
4525  int nvert = Geometry::NumVerts[geom];
4526  int nref_el = RG.RefGeoms.Size()/nvert;
4527  // If not, then store the offset and add to the size required
4528  point_matrices_offsets[id] = n_point_matrices[geom];
4529  n_point_matrices[geom] += nref_el;
4530  }
4531  }
4532 
4533  // Set up the sizes
4534  for (int geom = 0; geom < Geometry::NumGeom; ++geom)
4535  {
4536  int nmatrices = n_point_matrices[geom];
4537  int nvert = Geometry::NumVerts[geom];
4538  CoarseFineTr.point_matrices[geom].SetSize(Dim, nvert, nmatrices);
4539  }
4540 
4541  // Compute the point matrices and embeddings
4542  int el_fine = 0;
4543  for (int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4544  {
4545  Geometry::Type geom = orig_mesh.GetElementBaseGeometry(el_coarse);
4546  int ref = ref_factors[el_coarse];
4547  int offset = point_matrices_offsets[GeomRef(geom, ref)];
4548  int nvert = Geometry::NumVerts[geom];
4549  RefinedGeometry &RG = *refiner.Refine(geom, ref);
4550  for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
4551  {
4552  DenseMatrix &Pj = CoarseFineTr.point_matrices[geom](offset + j);
4553  for (int k = 0; k < nvert; k++)
4554  {
4555  int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
4556  const IntegrationPoint &ip = RG.RefPts[cid];
4557  ip.Get(Pj.GetColumn(k), Dim);
4558  }
4559 
4560  Embedding &emb = CoarseFineTr.embeddings[el_fine];
4561  emb.parent = el_coarse;
4562  emb.matrix = offset + j;
4563  ++el_fine;
4564  }
4565  }
4566 
4567  MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
4568 
4569  // The check below is disabled because is fails for parallel meshes with
4570  // interior "boundary" element that, when such "boundary" element is between
4571  // two elements on different processors.
4572  // MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
4573 }
4574 
4575 Mesh Mesh::MakeSimplicial(const Mesh &orig_mesh)
4576 {
4577  Mesh mesh;
4578  mesh.MakeSimplicial_(orig_mesh, NULL);
4579  return mesh;
4580 }
4581 
4582 void Mesh::MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
4583 {
4584  MFEM_VERIFY(const_cast<Mesh&>(orig_mesh).CheckElementOrientation(false) == 0,
4585  "Mesh::MakeSimplicial requires a properly oriented input mesh");
4586  MFEM_VERIFY(orig_mesh.Conforming(),
4587  "Mesh::MakeSimplicial does not support non-conforming meshes.")
4588 
4589  int dim = orig_mesh.Dimension();
4590  int sdim = orig_mesh.SpaceDimension();
4591 
4592  if (dim == 1)
4593  {
4594  Mesh copy(orig_mesh);
4595  Swap(copy, true);
4596  return;
4597  }
4598 
4599  int nv = orig_mesh.GetNV();
4600  int ne = orig_mesh.GetNE();
4601  int nbe = orig_mesh.GetNBE();
4602 
4603  static int num_subdivisions[Geometry::NUM_GEOMETRIES];
4604  num_subdivisions[Geometry::POINT] = 1;
4605  num_subdivisions[Geometry::SEGMENT] = 1;
4606  num_subdivisions[Geometry::TRIANGLE] = 1;
4607  num_subdivisions[Geometry::TETRAHEDRON] = 1;
4608  num_subdivisions[Geometry::SQUARE] = 2;
4609  num_subdivisions[Geometry::PRISM] = 3;
4610  num_subdivisions[Geometry::CUBE] = 6;
4611  // NOTE: some hexes may be subdivided into only 5 tets, so this is an
4612  // estimate only. The actual number of created tets may be less, so the
4613  // elements array will need to be shrunk after mesh creation.
4614  int new_ne = 0, new_nbe = 0;
4615  for (int i=0; i<ne; ++i)
4616  {
4617  new_ne += num_subdivisions[orig_mesh.GetElementBaseGeometry(i)];
4618  }
4619  for (int i=0; i<nbe; ++i)
4620  {
4621  new_nbe += num_subdivisions[orig_mesh.GetBdrElementBaseGeometry(i)];
4622  }
4623 
4624  InitMesh(dim, sdim, nv, new_ne, new_nbe);
4625 
4626  // Vertices of the new mesh are same as the original mesh
4627  NumOfVertices = nv;
4628  for (int i=0; i<nv; ++i)
4629  {
4630  vertices[i].SetCoords(dim, orig_mesh.vertices[i]());
4631  }
4632 
4633  // We need a global vertex numbering to identify which diagonals to split
4634  // (quad faces are split using the diagonal originating from the smallest
4635  // global vertex number). Use the supplied global numbering, if it is
4636  // non-NULL, otherwise use the local numbering.
4637  Array<int> vglobal_id;
4638  if (vglobal == NULL)
4639  {
4640  vglobal_id.SetSize(nv);
4641  for (int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4642  vglobal = vglobal_id.GetData();
4643  }
4644 
4645  constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4646  constexpr int quad_ntris = 2, prism_ntets = 3;
4647  static const int quad_trimap[2][nv_tri*quad_ntris] =
4648  {
4649  {
4650  0, 0,
4651  1, 2,
4652  2, 3
4653  },{
4654  0, 1,
4655  1, 2,
4656  3, 3
4657  }
4658  };
4659  static const int prism_rot[nv_prism*nv_prism] =
4660  {
4661  0, 1, 2, 3, 4, 5,
4662  1, 2, 0, 4, 5, 3,
4663  2, 0, 1, 5, 3, 4,
4664  3, 5, 4, 0, 2, 1,
4665  4, 3, 5, 1, 0, 2,
4666  5, 4, 3, 2, 1, 0
4667  };
4668  static const int prism_f[nv_quad] = {1, 2, 5, 4};
4669  static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4670  {
4671  {
4672  0, 0, 0,
4673  1, 1, 4,
4674  2, 5, 5,
4675  5, 4, 3
4676  },{
4677  0, 0, 0,
4678  1, 4, 4,
4679  2, 2, 5,
4680  4, 5, 3
4681  }
4682  };
4683  static const int hex_rot[nv_hex*nv_hex] =
4684  {
4685  0, 1, 2, 3, 4, 5, 6, 7,
4686  1, 0, 4, 5, 2, 3, 7, 6,
4687  2, 1, 5, 6, 3, 0, 4, 7,
4688  3, 0, 1, 2, 7, 4, 5, 6,
4689  4, 0, 3, 7, 5, 1, 2, 6,
4690  5, 1, 0, 4, 6, 2, 3, 7,
4691  6, 2, 1, 5, 7, 3, 0, 4,
4692  7, 3, 2, 6, 4, 0, 1, 5
4693  };
4694  static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4695  static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4696  static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4697  static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4698  static const int hex_tetmap0[nv_tet*5] =
4699  {
4700  0, 0, 0, 0, 2,
4701  1, 2, 2, 5, 7,
4702  2, 7, 3, 7, 5,
4703  5, 5, 7, 4, 6
4704  };
4705  static const int hex_tetmap1[nv_tet*6] =
4706  {
4707  0, 0, 1, 0, 0, 1,
4708  5, 1, 6, 7, 7, 7,
4709  7, 7, 7, 2, 1, 6,
4710  4, 5, 5, 3, 2, 2
4711  };
4712  static const int hex_tetmap2[nv_tet*6] =
4713  {
4714  0, 0, 0, 0, 0, 0,
4715  4, 3, 7, 1, 3, 6,
4716  5, 7, 4, 2, 6, 5,
4717  6, 6, 6, 5, 2, 2
4718  };
4719  static const int hex_tetmap3[nv_tet*6] =
4720  {
4721  0, 0, 0, 0, 1, 1,
4722  2, 3, 7, 5, 5, 6,
4723  3, 7, 4, 6, 6, 2,
4724  6, 6, 6, 4, 0, 0
4725  };
4726  static const int *hex_tetmaps[4] =
4727  {
4728  hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4729  };
4730 
4731  auto find_min = [](const int*a, int n) { return std::min_element(a,a+n)-a; };
4732 
4733  for (int i=0; i<ne; ++i)
4734  {
4735  const int *v = orig_mesh.elements[i]->GetVertices();
4736  const int attrib = orig_mesh.GetAttribute(i);
4737  const Geometry::Type orig_geom = orig_mesh.GetElementBaseGeometry(i);
4738 
4739  if (num_subdivisions[orig_geom] == 1)
4740  {
4741  // (num_subdivisions[orig_geom] == 1) implies that the element does
4742  // not need to be further split (it is either a segment, triangle,
4743  // or tetrahedron), and so it is left unchanged.
4744  Element *e = NewElement(orig_geom);
4745  e->SetAttribute(attrib);
4746  e->SetVertices(v);
4747  AddElement(e);
4748  }
4749  else if (orig_geom == Geometry::SQUARE)
4750  {
4751  for (int itri=0; itri<quad_ntris; ++itri)
4752  {
4754  e->SetAttribute(attrib);
4755  int *v2 = e->GetVertices();
4756  for (int iv=0; iv<nv_tri; ++iv)
4757  {
4758  v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4759  }
4760  AddElement(e);
4761  }
4762  }
4763  else if (orig_geom == Geometry::PRISM)
4764  {
4765  int vg[nv_prism];
4766  for (int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4767  // Rotate the vertices of the prism so that the smallest vertex index
4768  // is in the first place
4769  int irot = find_min(vg, nv_prism);
4770  for (int iv=0; iv<nv_prism; ++iv)
4771  {
4772  int jv = prism_rot[iv + irot*nv_prism];
4773  vg[iv] = v[jv];
4774  }
4775  // Two cases according to which diagonal splits third quad face
4776  int q[nv_quad];
4777  for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4778  int j = find_min(q, nv_quad);
4779  const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4780  for (int itet=0; itet<prism_ntets; ++itet)
4781  {
4783  e->SetAttribute(attrib);
4784  int *v2 = e->GetVertices();
4785  for (int iv=0; iv<nv_tet; ++iv)
4786  {
4787  v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4788  }
4789  AddElement(e);
4790  }
4791  }
4792  else if (orig_geom == Geometry::CUBE)
4793  {
4794  int vg[nv_hex];
4795  for (int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4796 
4797  // Rotate the vertices of the hex so that the smallest vertex index is
4798  // in the first place
4799  int irot = find_min(vg, nv_hex);
4800  for (int iv=0; iv<nv_hex; ++iv)
4801  {
4802  int jv = hex_rot[iv + irot*nv_hex];
4803  vg[iv] = v[jv];
4804  }
4805 
4806  int q[nv_quad];
4807  // Bitmask is three binary digits, each digit is 1 if the diagonal of
4808  // the corresponding face goes through the 7th vertex, and 0 if not.
4809  int bitmask = 0;
4810  int j;
4811  // First quad face
4812  for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4813  j = find_min(q, nv_quad);
4814  if (j == 0 || j == 2) { bitmask += 4; }
4815  // Second quad face
4816  for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4817  j = find_min(q, nv_quad);
4818  if (j == 1 || j == 3) { bitmask += 2; }
4819  // Third quad face
4820  for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4821  j = find_min(q, nv_quad);
4822  if (j == 0 || j == 2) { bitmask += 1; }
4823 
4824  // Apply rotations
4825  int nrot = num_rot[bitmask];
4826  for (int k=0; k<nrot; ++k)
4827  {
4828  int vtemp;
4829  vtemp = vg[1];
4830  vg[1] = vg[4];
4831  vg[4] = vg[3];
4832  vg[3] = vtemp;
4833  vtemp = vg[5];
4834  vg[5] = vg[7];
4835  vg[7] = vg[2];
4836  vg[2] = vtemp;
4837  }
4838 
4839  // Sum up nonzero bits in bitmask
4840  int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4841  int ntets = (ndiags == 0) ? 5 : 6;
4842  const int *tetmap = hex_tetmaps[ndiags];
4843  for (int itet=0; itet<ntets; ++itet)
4844  {
4846  e->SetAttribute(attrib);
4847  int *v2 = e->GetVertices();
4848  for (int iv=0; iv<nv_tet; ++iv)
4849  {
4850  v2[iv] = vg[tetmap[itet + iv*ntets]];
4851  }
4852  AddElement(e);
4853  }
4854  }
4855  }
4856  // In 3D, shrink the element array because some hexes have only 5 tets
4857  if (dim == 3) { elements.SetSize(NumOfElements); }
4858 
4859  for (int i=0; i<nbe; ++i)
4860  {
4861  const int *v = orig_mesh.boundary[i]->GetVertices();
4862  const int attrib = orig_mesh.GetBdrAttribute(i);
4863  const Geometry::Type orig_geom = orig_mesh.GetBdrElementBaseGeometry(i);
4864  if (num_subdivisions[orig_geom] == 1)
4865  {
4866  Element *be = NewElement(orig_geom);
4867  be->SetAttribute(attrib);
4868  be->SetVertices(v);
4869  AddBdrElement(be);
4870  }
4871  else if (orig_geom == Geometry::SQUARE)
4872  {
4873  int vg[nv_quad];
4874  for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4875  // Split quad according the smallest (global) vertex
4876  int iv_min = find_min(vg, nv_quad);
4877  int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4878  for (int itri=0; itri<quad_ntris; ++itri)
4879  {
4881  be->SetAttribute(attrib);
4882  int *v2 = be->GetVertices();
4883  for (int iv=0; iv<nv_tri; ++iv)
4884  {
4885  v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4886  }
4887  AddBdrElement(be);
4888  }
4889  }
4890  else
4891  {
4892  MFEM_ABORT("Unreachable");
4893  }
4894  }
4895 
4896  FinalizeTopology(false);
4897  sequence = orig_mesh.GetSequence();
4898  last_operation = orig_mesh.last_operation;
4899 
4900  MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
4901  MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
4902 }
4903 
4904 Mesh Mesh::MakePeriodic(const Mesh &orig_mesh, const std::vector<int> &v2v)
4905 {
4906  Mesh periodic_mesh(orig_mesh, true); // Make a copy of the original mesh
4907  const FiniteElementSpace *nodal_fes = orig_mesh.GetNodalFESpace();
4908  int nodal_order = nodal_fes ? nodal_fes->GetMaxElementOrder() : 1;
4909  periodic_mesh.SetCurvature(nodal_order, true);
4910 
4911  // renumber element vertices
4912  for (int i = 0; i < periodic_mesh.GetNE(); i++)
4913  {
4914  Element *el = periodic_mesh.GetElement(i);
4915  int *v = el->GetVertices();
4916  int nv = el->GetNVertices();
4917  for (int j = 0; j < nv; j++)
4918  {
4919  v[j] = v2v[v[j]];
4920  }
4921  }
4922  // renumber boundary element vertices
4923  for (int i = 0; i < periodic_mesh.GetNBE(); i++)
4924  {
4925  Element *el = periodic_mesh.GetBdrElement(i);
4926  int *v = el->GetVertices();
4927  int nv = el->GetNVertices();
4928  for (int j = 0; j < nv; j++)
4929  {
4930  v[j] = v2v[v[j]];
4931  }
4932  }
4933 
4934  periodic_mesh.RemoveUnusedVertices();
4935  return periodic_mesh;
4936 }
4937 
4939  const std::vector<Vector> &translations, double tol) const
4940 {
4941  int sdim = SpaceDimension();
4942 
4943  Vector coord(sdim), at(sdim), dx(sdim);
4944  Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4945  xMax = xMin = xDiff = 0.0;
4946 
4947  // Get a list of all vertices on the boundary
4948  set<int> bdr_v;
4949  for (int be = 0; be < GetNBE(); be++)
4950  {
4951  Array<int> dofs;
4952  GetBdrElementVertices(be,dofs);
4953 
4954  for (int i = 0; i < dofs.Size(); i++)
4955  {
4956  bdr_v.insert(dofs[i]);
4957 
4958  coord = GetVertex(dofs[i]);
4959  for (int j = 0; j < sdim; j++)
4960  {
4961  xMax[j] = max(xMax[j], coord[j]);
4962  xMin[j] = min(xMin[j], coord[j]);
4963  }
4964  }
4965  }
4966  add(xMax, -1.0, xMin, xDiff);
4967  double dia = xDiff.Norml2(); // compute mesh diameter
4968 
4969  // We now identify coincident vertices. Several originally distinct vertices
4970  // may become coincident under the periodic mapping. One of these vertices
4971  // will be identified as the "primary" vertex, and all other coincident
4972  // vertices will be considered as "replicas".
4973 
4974  // replica2primary[v] is the index of the primary vertex of replica `v`
4975  map<int, int> replica2primary;
4976  // primary2replicas[v] is a set of indices of replicas of primary vertex `v`
4977  map<int, set<int>> primary2replicas;
4978 
4979  // We begin with the assumption that all vertices are primary, and that there
4980  // are no replicas.
4981  for (int v : bdr_v) { primary2replicas[v]; }
4982 
4983  // Make `r` and all of `r`'s replicas be replicas of `p`. Delete `r` from the
4984  // list of primary vertices.
4985  auto make_replica = [&replica2primary, &primary2replicas](int r, int p)
4986  {
4987  if (r == p) { return; }
4988  primary2replicas[p].insert(r);
4989  replica2primary[r] = p;
4990  for (int s : primary2replicas[r])
4991  {
4992  primary2replicas[p].insert(s);
4993  replica2primary[s] = p;
4994  }
4995  primary2replicas.erase(r);
4996  };
4997 
4998  for (unsigned int i = 0; i < translations.size(); i++)
4999  {
5000  for (int vi : bdr_v)
5001  {
5002  coord = GetVertex(vi);
5003  add(coord, translations[i], at);
5004 
5005  for (int vj : bdr_v)
5006  {
5007  coord = GetVertex(vj);
5008  add(at, -1.0, coord, dx);
5009  if (dx.Norml2() > dia*tol) { continue; }
5010 
5011  // The two vertices vi and vj are coincident.
5012 
5013  // Are vertices `vi` and `vj` already primary?
5014  bool pi = primary2replicas.find(vi) != primary2replicas.end();
5015  bool pj = primary2replicas.find(vj) != primary2replicas.end();
5016 
5017  if (pi && pj)
5018  {
5019  // Both vertices are currently primary
5020  // Demote `vj` to be a replica of `vi`
5021  make_replica(vj, vi);
5022  }
5023  else if (pi && !pj)
5024  {
5025  // `vi` is primary and `vj` is a replica
5026  int owner_of_vj = replica2primary[vj];
5027  // Make `vi` and its replicas be replicas of `vj`'s owner
5028  make_replica(vi, owner_of_vj);
5029  }
5030  else if (!pi && pj)
5031  {
5032  // `vi` is currently a replica and `vj` is currently primary
5033  // Make `vj` and its replicas be replicas of `vi`'s owner
5034  int owner_of_vi = replica2primary[vi];
5035  make_replica(vj, owner_of_vi);
5036  }
5037  else
5038  {
5039  // Both vertices are currently replicas
5040  // Make `vj`'s owner and all of its owner's replicas be replicas
5041  // of `vi`'s owner
5042  int owner_of_vi = replica2primary[vi];
5043  int owner_of_vj = replica2primary[vj];
5044  make_replica(owner_of_vj, owner_of_vi);
5045  }
5046  break;
5047  }
5048  }
5049  }
5050 
5051  std::vector<int> v2v(GetNV());
5052  for (size_t i = 0; i < v2v.size(); i++)
5053  {
5054  v2v[i] = i;
5055  }
5056  for (auto &&r2p : replica2primary)
5057  {
5058  v2v[r2p.first] = r2p.second;
5059  }
5060  return v2v;
5061 }
5062 
5064 {
5065  if (NURBSext == NULL)
5066  {
5067  mfem_error("Mesh::KnotInsert : Not a NURBS mesh!");
5068  }
5069 
5070  if (kv.Size() != NURBSext->GetNKV())
5071  {
5072  mfem_error("Mesh::KnotInsert : KnotVector array size mismatch!");
5073  }
5074 
5076 
5077  NURBSext->KnotInsert(kv);
5078 
5079  last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
5080  sequence++;
5081 
5082  UpdateNURBS();
5083 }
5084 
5086 {
5087  if (NURBSext == NULL)
5088  {
5089  mfem_error("Mesh::KnotInsert : Not a NURBS mesh!");
5090  }
5091 
5092  if (kv.Size() != NURBSext->GetNKV())
5093  {
5094  mfem_error("Mesh::KnotInsert : KnotVector array size mismatch!");
5095  }
5096 
5098 
5099  NURBSext->KnotInsert(kv);
5100 
5101  last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
5102  sequence++;
5103 
5104  UpdateNURBS();
5105 }
5106 
5108 {
5109  // do not check for NURBSext since this method is protected
5111 
5113 
5114  last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
5115  sequence++;
5116 
5117  UpdateNURBS();
5118 }
5119 
5120 void Mesh::DegreeElevate(int rel_degree, int degree)
5121 {
5122  if (NURBSext == NULL)
5123  {
5124  mfem_error("Mesh::DegreeElevate : Not a NURBS mesh!");
5125  }
5126 
5128 
5129  NURBSext->DegreeElevate(rel_degree, degree);
5130 
5131  last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
5132  sequence++;
5133 
5134  UpdateNURBS();
5135 }
5136 
5138 {
5139  ResetLazyData();
5140 
5142 
5143  Dim = NURBSext->Dimension();
5144  spaceDim = Dim;
5145 
5146  if (NumOfElements != NURBSext->GetNE())
5147  {
5148  for (int i = 0; i < elements.Size(); i++)
5149  {
5150  FreeElement(elements[i]);
5151  }
5154  }
5155 
5156  if (NumOfBdrElements != NURBSext->GetNBE())
5157  {
5158  for (int i = 0; i < boundary.Size(); i++)
5159  {
5160  FreeElement(boundary[i]);
5161  }
5164  }
5165 
5166  Nodes->FESpace()->Update();
5167  Nodes->Update();
5169 
5170  if (NumOfVertices != NURBSext->GetNV())
5171  {
5173  vertices.SetSize(NumOfVertices);
5174  int vd = Nodes->VectorDim();
5175  for (int i = 0; i < vd; i++)
5176  {
5177  Vector vert_val;
5178  Nodes->GetNodalValues(vert_val, i+1);
5179  for (int j = 0; j < NumOfVertices; j++)
5180  {
5181  vertices[j](i) = vert_val(j);
5182  }
5183  }
5184  }
5185 
5186  if (el_to_edge)
5187  {
5189  if (Dim == 2)
5190  {
5191  GenerateFaces();
5192  }
5193  }
5194 
5195  if (el_to_face)
5196  {
5198  GenerateFaces();
5199  }
5200 }
5201 
5202 void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_knot)
5203 {
5204  SetEmpty();
5205 
5206  // Read MFEM NURBS mesh v1.0 format
5207  string ident;
5208 
5209  skip_comment_lines(input, '#');
5210 
5211  input >> ident; // 'dimension'
5212  input >> Dim;
5213  spaceDim = Dim;
5214 
5215  skip_comment_lines(input, '#');
5216 
5217  input >> ident; // 'elements'
5218  input >> NumOfElements;
5219  elements.SetSize(NumOfElements);
5220  for (int j = 0; j < NumOfElements; j++)
5221  {
5222  elements[j] = ReadElement(input);
5223  }
5224 
5225  skip_comment_lines(input, '#');
5226 
5227  input >> ident; // 'boundary'
5228  input >> NumOfBdrElements;
5229  boundary.SetSize(NumOfBdrElements);
5230  for (int j = 0; j < NumOfBdrElements; j++)
5231  {
5232  boundary[j] = ReadElement(input);
5233  }
5234 
5235  skip_comment_lines(input, '#');
5236 
5237  input >> ident; // 'edges'
5238  input >> NumOfEdges;
5239  edge_vertex = new Table(NumOfEdges, 2);
5240  edge_to_knot.SetSize(NumOfEdges);
5241  for (int j = 0; j < NumOfEdges; j++)
5242  {
5243  int *v = edge_vertex->GetRow(j);
5244  input >> edge_to_knot[j] >> v[0] >> v[1];
5245  if (v[0] > v[1])
5246  {
5247  edge_to_knot[j] = -1 - edge_to_knot[j];
5248  }
5249  }
5250 
5251  skip_comment_lines(input, '#');
5252 
5253  input >> ident; // 'vertices'
5254  input >> NumOfVertices;
5255  vertices.SetSize(0);
5256 
5257  FinalizeTopology();
5258  CheckBdrElementOrientation(); // check and fix boundary element orientation
5259 }
5260 
5262 {
5263  if (p.Size() >= v.Size())
5264  {
5265  for (int d = 0; d < v.Size(); d++)
5266  {
5267  v(d) = p(d);
5268  }
5269  }
5270  else
5271  {
5272  int d;
5273  for (d = 0; d < p.Size(); d++)
5274  {
5275  v(d) = p(d);
5276  }
5277  for ( ; d < v.Size(); d++)
5278  {
5279  v(d) = 0.0;
5280  }
5281  }
5282 }
5283 
5284 void Mesh::GetNodes(GridFunction &nodes) const
5285 {
5286  if (Nodes == NULL || Nodes->FESpace() != nodes.FESpace())
5287  {
5288  const int newSpaceDim = nodes.FESpace()->GetVDim();
5290  nodes.ProjectCoefficient(xyz);
5291  }
5292  else
5293  {
5294  nodes = *Nodes;
5295  }
5296 }
5297 
5299 {
5300  GridFunction *nodes = new GridFunction(nfes);
5301  SetNodalGridFunction(nodes, true);
5302 }
5303 
5305 {
5306  if (Nodes)
5307  {
5309  if (dynamic_cast<const H1_FECollection*>(fec)
5310  || dynamic_cast<const L2_FECollection*>(fec))
5311  {
5312  return;
5313  }
5314  else // Mesh using a legacy FE_Collection
5315  {
5316  const int order = GetNodalFESpace()->GetElementOrder(0);
5317  if (NURBSext)
5318  {
5319 #ifndef MFEM_USE_MPI
5320  const bool warn = true;
5321 #else
5322  ParMesh *pmesh = dynamic_cast<ParMesh*>(this);
5323  const bool warn = !pmesh || pmesh->GetMyRank() == 0;
5324 #endif
5325  if (warn)
5326  {
5327  MFEM_WARNING("converting NURBS mesh to order " << order <<
5328  " H1-continuous mesh!\n "
5329  "If this is the desired behavior, you can silence"
5330  " this warning by converting\n "
5331  "the NURBS mesh to high-order mesh in advance by"
5332  " calling the method\n "
5333  "Mesh::SetCurvature().");
5334  }
5335  }
5336  SetCurvature(order, false, -1, Ordering::byVDIM);
5337  }
5338  }
5339  else // First order H1 mesh
5340  {
5341  SetCurvature(1, false, -1, Ordering::byVDIM);
5342  }
5343 }
5344 
5345 void Mesh::SetNodalGridFunction(GridFunction *nodes, bool make_owner)
5346 {
5347  GetNodes(*nodes);
5348  NewNodes(*nodes, make_owner);
5349 }
5350 
5352 {
5353  return ((Nodes) ? Nodes->FESpace() : NULL);
5354 }
5355 
5356 void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
5357 {
5358  space_dim = (space_dim == -1) ? spaceDim : space_dim;
5360  if (discont)
5361  {
5362  const int type = 1; // Gauss-Lobatto points
5363  nfec = new L2_FECollection(order, Dim, type);
5364  }
5365  else
5366  {
5367  nfec = new H1_FECollection(order, Dim);
5368  }
5369  FiniteElementSpace* nfes = new FiniteElementSpace(this, nfec, space_dim,
5370  ordering);
5371  SetNodalFESpace(nfes);
5372  Nodes->MakeOwner(nfec);
5373 }
5374 
5376 {
5377  MFEM_ASSERT(nodes != NULL, "");
5378  for (int i = 0; i < spaceDim; i++)
5379  {
5380  Vector vert_val;
5381  nodes->GetNodalValues(vert_val, i+1);
5382  for (int j = 0; j < NumOfVertices; j++)
5383  {
5384  vertices[j](i) = vert_val(j);
5385  }
5386  }
5387 }
5388 
5390 {
5391  switch (Dim)
5392  {
5393  case 1: return GetNV();
5394  case 2: return GetNEdges();
5395  case 3: return GetNFaces();
5396  }
5397  return 0;
5398 }
5399 
5401 {
5402  return faces_info.Size();
5403 }
5404 
5406 {
5407  const bool isInt = type==FaceType::Interior;
5408  int &nf = isInt ? nbInteriorFaces : nbBoundaryFaces;
5409  if (nf<0)
5410  {
5411  nf = 0;
5412  for (int f = 0; f < GetNumFacesWithGhost(); ++f)
5413  {
5415  if ( face.IsOfFaceType(type) )
5416  {
5417  if (face.IsNonconformingCoarse())
5418  {
5419  // We don't count nonconforming coarse faces.
5420  continue;
5421  }
5422  nf++;
5423  }
5424  }
5425  }
5426  return nf;
5427 }
5428 
5429 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5430 static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
5431 #endif
5432 
5434 {
5435  int i, j, k, wo = 0, fo = 0;
5436  double *v[4];
5437 
5438  if (Dim == 2 && spaceDim == 2)
5439  {
5440  DenseMatrix J(2, 2);
5441 
5442  for (i = 0; i < NumOfElements; i++)
5443  {
5444  int *vi = elements[i]->GetVertices();
5445  if (Nodes == NULL)
5446  {
5447  for (j = 0; j < 3; j++)
5448  {
5449  v[j] = vertices[vi[j]]();
5450  }
5451  for (j = 0; j < 2; j++)
5452  for (k = 0; k < 2; k++)
5453  {
5454  J(j, k) = v[j+1][k] - v[0][k];
5455  }
5456  }
5457  else
5458  {
5459  // only check the Jacobian at the center of the element
5460  GetElementJacobian(i, J);
5461  }
5462  if (J.Det() < 0.0)
5463  {
5464  if (fix_it)
5465  {
5466  switch (GetElementType(i))
5467  {
5468  case Element::TRIANGLE:
5469  mfem::Swap(vi[0], vi[1]);
5470  break;
5472  mfem::Swap(vi[1], vi[3]);
5473  break;
5474  default:
5475  MFEM_ABORT("Invalid 2D element type \""
5476  << GetElementType(i) << "\"");
5477  break;
5478  }
5479  fo++;
5480  }
5481  wo++;
5482  }
5483  }
5484  }
5485 
5486  if (Dim == 3)
5487  {
5488  DenseMatrix J(3, 3);
5489 
5490  for (i = 0; i < NumOfElements; i++)
5491  {
5492  int *vi = elements[i]->GetVertices();
5493  switch (GetElementType(i))
5494  {
5495  case Element::TETRAHEDRON:
5496  if (Nodes == NULL)
5497  {
5498  for (j = 0; j < 4; j++)
5499  {
5500  v[j] = vertices[vi[j]]();
5501  }
5502  for (j = 0; j < 3; j++)
5503  for (k = 0; k < 3; k++)
5504  {
5505  J(j, k) = v[j+1][k] - v[0][k];
5506  }
5507  }
5508  else
5509  {
5510  // only check the Jacobian at the center of the element
5511  GetElementJacobian(i, J);
5512  }
5513  if (J.Det() < 0.0)
5514  {
5515  wo++;
5516  if (fix_it)
5517  {
5518  mfem::Swap(vi[0], vi[1]);
5519  fo++;
5520  }
5521  }
5522  break;
5523 
5524  case Element::WEDGE:
5525  // only check the Jacobian at the center of the element
5526  GetElementJacobian(i, J);
5527  if (J.Det() < 0.0)
5528  {
5529  wo++;
5530  if (fix_it)
5531  {
5532  // how?
5533  }
5534  }
5535  break;
5536 
5537  case Element::PYRAMID:
5538  // only check the Jacobian at the center of the element
5539  GetElementJacobian(i, J);
5540  if (J.Det() < 0.0)
5541  {
5542  wo++;
5543  if (fix_it)
5544  {
5545  // how?
5546  }
5547  }
5548  break;
5549 
5550  case Element::HEXAHEDRON:
5551  // only check the Jacobian at the center of the element
5552  GetElementJacobian(i, J);
5553  if (J.Det() < 0.0)
5554  {
5555  wo++;
5556  if (fix_it)
5557  {
5558  // how?
5559  }
5560  }
5561  break;
5562 
5563  default:
5564  MFEM_ABORT("Invalid 3D element type \""
5565  << GetElementType(i) << "\"");
5566  break;
5567  }
5568  }
5569  }
5570 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5571  if (wo > 0)
5572  {
5573  mfem::out << "Elements with wrong orientation: " << wo << " / "
5574  << NumOfElements << " (" << fixed_or_not[(wo == fo) ? 0 : 1]
5575  << ")" << endl;
5576  }
5577 #endif
5578  return wo;
5579 }
5580 
5581 int Mesh::GetTriOrientation(const int *base, const int *test)
5582 {
5583  // Static method.
5584  // This function computes the index 'j' of the permutation that transforms
5585  // test into base: test[tri_orientation[j][i]]=base[i].
5586  // tri_orientation = Geometry::Constants<Geometry::TRIANGLE>::Orient
5587  int orient;
5588 
5589  if (test[0] == base[0])
5590  if (test[1] == base[1])
5591  {
5592  orient = 0; // (0, 1, 2)
5593  }
5594  else
5595  {
5596  orient = 5; // (0, 2, 1)
5597  }
5598  else if (test[0] == base[1])
5599  if (test[1] == base[0])
5600  {
5601  orient = 1; // (1, 0, 2)
5602  }
5603  else
5604  {
5605  orient = 2; // (1, 2, 0)
5606  }
5607  else // test[0] == base[2]
5608  if (test[1] == base[0])
5609  {
5610  orient = 4; // (2, 0, 1)
5611  }
5612  else
5613  {
5614  orient = 3; // (2, 1, 0)
5615  }
5616 
5617 #ifdef MFEM_DEBUG
5618  const int *aor = tri_t::Orient[orient];
5619  for (int j = 0; j < 3; j++)
5620  if (test[aor[j]] != base[j])
5621  {
5622  mfem_error("Mesh::GetTriOrientation(...)");
5623  }
5624 #endif
5625 
5626  return orient;
5627 }
5628 
5629 int Mesh::GetQuadOrientation(const int *base, const int *test)
5630 {
5631  int i;
5632 
5633  for (i = 0; i < 4; i++)
5634  if (test[i] == base[0])
5635  {
5636  break;
5637  }
5638 
5639 #ifdef MFEM_DEBUG
5640  int orient;
5641  if (test[(i+1)%4] == base[1])
5642  {
5643  orient = 2*i;
5644  }
5645  else
5646  {
5647  orient = 2*i+1;
5648  }
5649  const int *aor = quad_t::Orient[orient];
5650  for (int j = 0; j < 4; j++)
5651  if (test[aor[j]] != base[j])
5652  {
5653  mfem::err << "Mesh::GetQuadOrientation(...)" << endl;
5654  mfem::err << " base = [";
5655  for (int k = 0; k < 4; k++)
5656  {
5657  mfem::err << " " << base[k];
5658  }
5659  mfem::err << " ]\n test = [";
5660  for (int k = 0; k < 4; k++)
5661  {
5662  mfem::err << " " << test[k];
5663  }
5664  mfem::err << " ]" << endl;
5665  mfem_error();
5666  }
5667 #endif
5668 
5669  if (test[(i+1)%4] == base[1])
5670  {
5671  return 2*i;
5672  }
5673 
5674  return 2*i+1;
5675 }
5676 
5677 int Mesh::GetTetOrientation(const int *base, const int *test)
5678 {
5679  // Static method.
5680  // This function computes the index 'j' of the permutation that transforms
5681  // test into base: test[tet_orientation[j][i]]=base[i].
5682  // tet_orientation = Geometry::Constants<Geometry::TETRAHEDRON>::Orient
5683  int orient;
5684 
5685  if (test[0] == base[0])
5686  if (test[1] == base[1])
5687  if (test[2] == base[2])
5688  {
5689  orient = 0; // (0, 1, 2, 3)
5690  }
5691  else
5692  {
5693  orient = 1; // (0, 1, 3, 2)
5694  }
5695  else if (test[2] == base[1])
5696  if (test[3] == base[2])
5697  {
5698  orient = 2; // (0, 2, 3, 1)
5699  }
5700  else
5701  {
5702  orient = 3; // (0, 2, 1, 3)
5703  }
5704  else // test[3] == base[1]
5705  if (test[1] == base[2])
5706  {
5707  orient = 4; // (0, 3, 1, 2)
5708  }
5709  else
5710  {
5711  orient = 5; // (0, 3, 2, 1)
5712  }
5713  else if (test[1] == base[0])
5714  if (test[2] == base[1])
5715  if (test[0] == base[2])
5716  {
5717  orient = 6; // (1, 2, 0, 3)
5718  }
5719  else
5720  {
5721  orient = 7; // (1, 2, 3, 0)
5722  }
5723  else if (test[3] == base[1])
5724  if (test[2] == base[2])
5725  {
5726  orient = 8; // (1, 3, 2, 0)
5727  }
5728  else
5729  {
5730  orient = 9; // (1, 3, 0, 2)
5731  }
5732  else // test[0] == base[1]
5733  if (test[3] == base[2])
5734  {
5735  orient = 10; // (1, 0, 3, 2)
5736  }
5737  else
5738  {
5739  orient = 11; // (1, 0, 2, 3)
5740  }
5741  else if (test[2] == base[0])
5742  if (test[3] == base[1])
5743  if (test[0] == base[2])
5744  {
5745  orient = 12; // (2, 3, 0, 1)
5746  }
5747  else
5748  {
5749  orient = 13; // (2, 3, 1, 0)
5750  }
5751  else if (test[0] == base[1])
5752  if (test[1] == base[2])
5753  {
5754  orient = 14; // (2, 0, 1, 3)
5755  }
5756  else
5757  {
5758  orient = 15; // (2, 0, 3, 1)
5759  }
5760  else // test[1] == base[1]
5761  if (test[3] == base[2])
5762  {
5763  orient = 16; // (2, 1, 3, 0)
5764  }
5765  else
5766  {
5767  orient = 17; // (2, 1, 0, 3)
5768  }
5769  else // (test[3] == base[0])
5770  if (test[0] == base[1])
5771  if (test[2] == base[2])
5772  {
5773  orient = 18; // (3, 0, 2, 1)
5774  }
5775  else
5776  {
5777  orient = 19; // (3, 0, 1, 2)
5778  }
5779  else if (test[1] == base[1])
5780  if (test[0] == base[2])
5781  {
5782  orient = 20; // (3, 1, 0, 2)
5783  }
5784  else
5785  {
5786  orient = 21; // (3, 1, 2, 0)
5787  }
5788  else // test[2] == base[1]
5789  if (test[1] == base[2])
5790  {
5791  orient = 22; // (3, 2, 1, 0)
5792  }
5793  else
5794  {
5795  orient = 23; // (3, 2, 0, 1)
5796  }
5797 
5798 #ifdef MFEM_DEBUG
5799  const int *aor = tet_t::Orient[orient];
5800  for (int j = 0; j < 4; j++)
5801  if (test[aor[j]] != base[j])
5802  {
5803  mfem_error("Mesh::GetTetOrientation(...)");
5804  }
5805 #endif
5806 
5807  return orient;
5808 }
5809 
5811 {
5812  int wo = 0; // count wrong orientations
5813 
5814  if (Dim == 2)
5815  {
5816  if (el_to_edge == NULL) // edges were not generated
5817  {
5818  el_to_edge = new Table;
5820  GenerateFaces(); // 'Faces' in 2D refers to the edges
5821  }
5822  for (int i = 0; i < NumOfBdrElements; i++)
5823  {
5824  if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
5825  {
5826  int *bv = boundary[i]->GetVertices();
5827  int *fv = faces[be_to_edge[i]]->GetVertices();
5828  if (bv[0] != fv[0])
5829  {
5830  if (fix_it)
5831  {
5832  mfem::Swap<int>(bv[0], bv[1]);
5833  }
5834  wo++;
5835  }
5836  }
5837  }
5838  }
5839 
5840  if (Dim == 3)
5841  {
5842  for (int i = 0; i < NumOfBdrElements; i++)
5843  {
5844  const int fi = be_to_face[i];
5845 
5846  if (faces_info[fi].Elem2No >= 0) { continue; }
5847 
5848  // boundary face
5849  int *bv = boundary[i]->GetVertices();
5850  // Make sure the 'faces' are generated:
5851  MFEM_ASSERT(fi < faces.Size(), "internal error");
5852  const int *fv = faces[fi]->GetVertices();
5853  int orientation; // orientation of the bdr. elem. w.r.t. the
5854  // corresponding face element (that's the base)
5855  const Element::Type bdr_type = GetBdrElementType(i);
5856  switch (bdr_type)
5857  {
5858  case Element::TRIANGLE:
5859  {
5860  orientation = GetTriOrientation(fv, bv);
5861  break;
5862  }
5864  {
5865  orientation = GetQuadOrientation(fv, bv);
5866  break;
5867  }
5868  default:
5869  MFEM_ABORT("Invalid 2D boundary element type \""
5870  << bdr_type << "\"");
5871  orientation = 0; // suppress a warning
5872  break;
5873  }
5874 
5875  if (orientation % 2 == 0) { continue; }
5876  wo++;
5877  if (!fix_it) { continue; }
5878 
5879  switch (bdr_type)
5880  {
5881  case Element::TRIANGLE:
5882  {
5883  // swap vertices 0 and 1 so that we don't change the marked edge:
5884  // (0,1,2) -> (1,0,2)
5885  mfem::Swap<int>(bv[0], bv[1]);
5886  if (bel_to_edge)
5887  {
5888  int *be = bel_to_edge->GetRow(i);
5889  mfem::Swap<int>(be[1], be[2]);
5890  }
5891  break;
5892  }
5894  {
5895  mfem::Swap<int>(bv[0], bv[2]);
5896  if (bel_to_edge)
5897  {
5898  int *be = bel_to_edge->GetRow(i);
5899  mfem::Swap<int>(be[0], be[1]);
5900  mfem::Swap<int>(be[2], be[3]);
5901  }
5902  break;
5903  }
5904  default: // unreachable
5905  break;
5906  }
5907  }
5908  }
5909  // #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5910 #ifdef MFEM_DEBUG
5911  if (wo > 0)
5912  {
5913  mfem::out << "Boundary elements with wrong orientation: " << wo << " / "
5914  << NumOfBdrElements << " (" << fixed_or_not[fix_it ? 0 : 1]
5915  << ")" << endl;
5916  }
5917 #endif
5918  return wo;
5919 }
5920 
5922 {
5923  MFEM_ASSERT(0 <= dim && dim <= Dim, "invalid dim: " << dim);
5924  int num_geoms = 0;
5925  for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
5926  {
5927  if (HasGeometry(Geometry::Type(g))) { num_geoms++; }
5928  }
5929  return num_geoms;
5930 }
5931 
5933 {
5934  MFEM_ASSERT(0 <= dim && dim <= Dim, "invalid dim: " << dim);
5935  el_geoms.SetSize(0);
5936  for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
5937  {
5938  if (HasGeometry(Geometry::Type(g)))
5939  {
5940  el_geoms.Append(Geometry::Type(g));
5941  }
5942  }
5943 }
5944 
5945 void Mesh::GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const
5946 {
5947  if (el_to_edge)
5948  {
5949  el_to_edge->GetRow(i, edges);
5950  }
5951  else
5952  {
5953  mfem_error("Mesh::GetElementEdges(...) element to edge table "
5954  "is not generated.");
5955  }
5956 
5957  const int *v = elements[i]->GetVertices();
5958  const int ne = elements[i]->GetNEdges();
5959  cor.SetSize(ne);
5960  for (int j = 0; j < ne; j++)
5961  {
5962  const int *e = elements[i]->GetEdgeVertices(j);
5963  cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5964  }
5965 }
5966 
5967 void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const
5968 {
5969  if (Dim == 2)
5970  {
5971  edges.SetSize(1);
5972  cor.SetSize(1);
5973  edges[0] = be_to_edge[i];
5974  const int *v = boundary[i]->GetVertices();
5975  cor[0] = (v[0] < v[1]) ? (1) : (-1);
5976  }
5977  else if (Dim == 3)
5978  {
5979  if (bel_to_edge)
5980  {
5981  bel_to_edge->GetRow(i, edges);
5982  }
5983  else
5984  {
5985  mfem_error("Mesh::GetBdrElementEdges(...)");
5986  }
5987 
5988  const int *v = boundary[i]->GetVertices();
5989  const int ne = boundary[i]->GetNEdges();
5990  cor.SetSize(ne);
5991  for (int j = 0; j < ne; j++)
5992  {
5993  const int *e = boundary[i]->GetEdgeVertices(j);
5994  cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5995  }
5996  }
5997 }
5998 
5999 void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
6000 {
6001  if (Dim == 2)
6002  {
6003  edges.SetSize(1);
6004  edges[0] = i;
6005  o.SetSize(1);
6006  const int *v = faces[i]->GetVertices();
6007  o[0] = (v[0] < v[1]) ? (1) : (-1);
6008  }
6009 
6010  if (Dim != 3)
6011  {
6012  return;
6013  }
6014 
6015  GetFaceEdgeTable(); // generate face_edge Table (if not generated)
6016 
6017  face_edge->GetRow(i, edges);
6018 
6019  const int *v = faces[i]->GetVertices();
6020  const int ne = faces[i]->GetNEdges();
6021  o.SetSize(ne);
6022  for (int j = 0; j < ne; j++)
6023  {
6024  const int *e = faces[i]->GetEdgeVertices(j);
6025  o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6026  }
6027 }
6028 
6029 void Mesh::GetEdgeVertices(int i, Array<int> &vert) const
6030 {
6031  // the two vertices are sorted: vert[0] < vert[1]
6032  // this is consistent with the global edge orientation
6033  // generate edge_vertex Table (if not generated)
6034  if (!edge_vertex) { GetEdgeVertexTable(); }
6035  edge_vertex->GetRow(i, vert);
6036 }
6037 
6039 {
6040  if (face_edge)
6041  {
6042  return face_edge;
6043  }
6044 
6045  if (Dim != 3)
6046  {
6047  return NULL;
6048  }
6049 
6050 #ifdef MFEM_DEBUG
6051  if (faces.Size() != NumOfFaces)
6052  {
6053  mfem_error("Mesh::GetFaceEdgeTable : faces were not generated!");
6054  }
6055 #endif
6056 
6057  DSTable v_to_v(NumOfVertices);
6058  GetVertexToVertexTable(v_to_v);
6059 
6060  face_edge = new Table;
6062 
6063  return (face_edge);
6064 }
6065 
6067 {
6068  if (edge_vertex)
6069  {
6070  return edge_vertex;
6071  }
6072 
6073  DSTable v_to_v(NumOfVertices);
6074  GetVertexToVertexTable(v_to_v);
6075 
6076  int nedges = v_to_v.NumberOfEntries();
6077  edge_vertex = new Table(nedges, 2);
6078  for (int i = 0; i < NumOfVertices; i++)
6079  {
6080  for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
6081  {
6082  int j = it.Index();
6083  edge_vertex->Push(j, i);
6084  edge_vertex->Push(j, it.Column());
6085  }
6086  }
6087  edge_vertex->Finalize();
6088 
6089  return edge_vertex;
6090 }
6091 
6093 {
6094  int i, j, nv, *v;
6095 
6096  Table *vert_elem = new Table;
6097 
6098  vert_elem->MakeI(NumOfVertices);
6099 
6100  for (i = 0; i < NumOfElements; i++)
6101  {
6102  nv = elements[i]->GetNVertices();
6103  v = elements[i]->GetVertices();
6104  for (j = 0; j < nv; j++)
6105  {
6106  vert_elem->AddAColumnInRow(v[j]);
6107  }
6108  }
6109 
6110  vert_elem->MakeJ();
6111 
6112  for (i = 0; i < NumOfElements; i++)
6113  {
6114  nv = elements[i]->GetNVertices();
6115  v = elements[i]->GetVertices();
6116  for (j = 0; j < nv; j++)
6117  {
6118  vert_elem->AddConnection(v[j], i);
6119  }
6120  }
6121 
6122  vert_elem->ShiftUpI();
6123 
6124  return vert_elem;
6125 }
6126 
6128 {
6129  Table *face_elem = new Table;
6130 
6131  face_elem->MakeI(faces_info.Size());
6132 
6133  for (int i = 0; i < faces_info.Size(); i++)
6134  {
6135  if (faces_info[i].Elem2No >= 0)
6136  {
6137  face_elem->AddColumnsInRow(i, 2);
6138  }
6139  else
6140  {
6141  face_elem->AddAColumnInRow(i);
6142  }
6143  }
6144 
6145  face_elem->MakeJ();
6146 
6147  for (int i = 0; i < faces_info.Size(); i++)
6148  {
6149  face_elem->AddConnection(i, faces_info[i].Elem1No);
6150  if (faces_info[i].Elem2No >= 0)
6151  {
6152  face_elem->AddConnection(i, faces_info[i].Elem2No);
6153  }
6154  }
6155 
6156  face_elem->ShiftUpI();
6157 
6158  return face_elem;
6159 }
6160 
6161 void Mesh::GetElementFaces(int i, Array<int> &el_faces, Array<int> &ori) const
6162 {
6163  MFEM_VERIFY(el_to_face != NULL, "el_to_face not generated");
6164 
6165  el_to_face->GetRow(i, el_faces);
6166 
6167  int n = el_faces.Size();
6168  ori.SetSize(n);
6169 
6170  for (int j = 0; j < n; j++)
6171  {
6172  if (faces_info[el_faces[j]].Elem1No == i)
6173  {
6174  ori[j] = faces_info[el_faces[j]].Elem1Inf % 64;
6175  }
6176  else
6177  {
6178  MFEM_ASSERT(faces_info[el_faces[j]].Elem2No == i, "internal error");
6179  ori[j] = faces_info[el_faces[j]].Elem2Inf % 64;
6180  }
6181  }
6182 }
6183 
6185 {
6186  if (face_to_elem == NULL)
6187  {
6189  }
6190 
6191  Array<int> elem_faces;
6192  Array<int> ori;
6193  GetElementFaces(elem, elem_faces, ori);
6194 
6195  Array<int> nghb;
6196  for (auto f : elem_faces)
6197  {
6198  Array<int> row;
6199  face_to_elem->GetRow(f, row);
6200  for (auto r : row)
6201  {
6202  nghb.Append(r);
6203  }
6204  }
6205 
6206  nghb.Sort();
6207  nghb.Unique();
6208 
6209  return nghb;
6210 }
6211 
6212 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
6213 {
6214  const int *bv, *fv;
6215 
6216  *f = be_to_face[i];
6217  bv = boundary[i]->GetVertices();
6218  fv = faces[be_to_face[i]]->GetVertices();
6219 
6220  // find the orientation of the bdr. elem. w.r.t.
6221  // the corresponding face element (that's the base)
6222  switch (GetBdrElementType(i))
6223  {
6224  case Element::TRIANGLE:
6225  *o = GetTriOrientation(fv, bv);
6226  break;
6228  *o = GetQuadOrientation(fv, bv);
6229  break;
6230  default:
6231  MFEM_ABORT("invalid geometry");
6232  }
6233 }
6234 
6236 {
6237  switch (Dim)
6238  {
6239  case 1: return boundary[i]->GetVertices()[0];
6240  case 2: return be_to_edge[i];
6241  case 3: return be_to_face[i];
6242  default: MFEM_ABORT("invalid dimension!");
6243  }
6244  return -1;
6245 }
6246 
6247 void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
6248 {
6249  int fid = GetBdrElementEdgeIndex(bdr_el);
6250 
6251  const FaceInfo &fi = faces_info[fid];
6252  MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
6253 
6254  const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
6255  const int *bv = boundary[bdr_el]->GetVertices();
6256  int ori;
6257  switch (GetBdrElementGeometry(bdr_el))
6258  {
6259  case Geometry::POINT: ori = 0; break;
6260  case Geometry::SEGMENT: ori = (fv[0] == bv[0]) ? 0 : 1; break;
6261  case Geometry::TRIANGLE: ori = GetTriOrientation(fv, bv); break;
6262  case Geometry::SQUARE: ori = GetQuadOrientation(fv, bv); break;
6263  default: MFEM_ABORT("boundary element type not implemented"); ori = 0;
6264  }
6265  el = fi.Elem1No;
6266  info = fi.Elem1Inf + ori;
6267 }
6268 
6269 void Mesh::GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
6270 {
6271  int fid = GetBdrElementEdgeIndex(bdr_el);
6272 
6273  const FaceInfo &fi = faces_info[fid];
6274  MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
6275 
6276  const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
6277  const int *bv = boundary[bdr_el]->GetVertices();
6278  int ori;
6279  switch (GetBdrElementGeometry(bdr_el))
6280  {
6281  case Geometry::POINT: ori = 0; break;
6282  case Geometry::SEGMENT: ori = (fv[0] == bv[0]) ? 0 : 1; break;
6283  case Geometry::TRIANGLE: ori = GetTriOrientation(bv, fv); break;
6284  case Geometry::SQUARE: ori = GetQuadOrientation(bv, fv); break;
6285  default: MFEM_ABORT("boundary element type not implemented"); ori = 0;
6286  }
6287  el = fi.Elem1No;
6288  info = fi.Elem1Inf + ori;
6289 }
6290 
6292 {
6293  return elements[i]->GetType();
6294 }
6295 
6297 {
6298  return boundary[i]->GetType();
6299 }
6300 
6301 void Mesh::GetPointMatrix(int i, DenseMatrix &pointmat) const
6302 {
6303  int k, j, nv;
6304  const int *v;
6305 
6306  v = elements[i]->GetVertices();
6307  nv = elements[i]->GetNVertices();
6308 
6309  pointmat.SetSize(spaceDim, nv);
6310  for (k = 0; k < spaceDim; k++)
6311  {
6312  for (j = 0; j < nv; j++)
6313  {
6314  pointmat(k, j) = vertices[v[j]](k);
6315  }
6316  }
6317 }
6318 
6319 void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
6320 {
6321  int k, j, nv;
6322  const int *v;
6323 
6324  v = boundary[i]->GetVertices();
6325  nv = boundary[i]->GetNVertices();
6326 
6327  pointmat.SetSize(spaceDim, nv);
6328  for (k = 0; k < spaceDim; k++)
6329  for (j = 0; j < nv; j++)
6330  {
6331  pointmat(k, j) = vertices[v[j]](k);
6332  }
6333 }
6334 
6335 double Mesh::GetLength(int i, int j) const
6336 {
6337  const double *vi = vertices[i]();
6338  const double *vj = vertices[j]();
6339  double length = 0.;
6340 
6341  for (int k = 0; k < spaceDim; k++)
6342  {
6343  length += (vi[k]-vj[k])*(vi[k]-vj[k]);
6344  }
6345 
6346  return sqrt(length);
6347 }
6348 
6349 // static method
6351  const DSTable &v_to_v, Table &el_to_edge)
6352 {
6353  el_to_edge.MakeI(elem_array.Size());
6354  for (int i = 0; i < elem_array.Size(); i++)
6355  {
6356  el_to_edge.AddColumnsInRow(i, elem_array[i]->GetNEdges());
6357  }
6358  el_to_edge.MakeJ();
6359  for (int i = 0; i < elem_array.Size(); i++)
6360  {
6361  const int *v = elem_array[i]->GetVertices();
6362  const int ne = elem_array[i]->GetNEdges();
6363  for (int j = 0; j < ne; j++)
6364  {
6365  const int *e = elem_array[i]->GetEdgeVertices(j);
6366  el_to_edge.AddConnection(i, v_to_v(v[e[0]], v[e[1]]));
6367  }
6368  }
6369  el_to_edge.ShiftUpI();
6370 }
6371 
6373 {
6374  if (edge_vertex)
6375  {
6376  for (int i = 0; i < edge_vertex->Size(); i++)
6377  {
6378  const int *v = edge_vertex->GetRow(i);
6379  v_to_v.Push(v[0], v[1]);
6380  }
6381  }
6382  else
6383  {
6384  for (int i = 0; i < NumOfElements; i++)
6385  {
6386  const int *v = elements[i]->GetVertices();
6387  const int ne = elements[i]->GetNEdges();
6388  for (int j = 0; j < ne; j++)
6389  {
6390  const int *e = elements[i]->GetEdgeVertices(j);
6391  v_to_v.Push(v[e[0]], v[e[1]]);
6392  }
6393  }
6394  }
6395 }
6396 
6398 {
6399  int i, NumberOfEdges;
6400 
6401  DSTable v_to_v(NumOfVertices);
6402  GetVertexToVertexTable(v_to_v);
6403 
6404  NumberOfEdges = v_to_v.NumberOfEntries();
6405 
6406  // Fill the element to edge table
6407  GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
6408 
6409  if (Dim == 2)
6410  {
6411  // Initialize the indices for the boundary elements.
6412  be_to_f.SetSize(NumOfBdrElements);
6413  for (i = 0; i < NumOfBdrElements; i++)
6414  {
6415  const int *v = boundary[i]->GetVertices();
6416  be_to_f[i] = v_to_v(v[0], v[1]);
6417  }
6418  }
6419  else if (Dim == 3)
6420  {
6421  if (bel_to_edge == NULL)
6422  {
6423  bel_to_edge = new Table;
6424  }
6426  }
6427  else
6428  {
6429  mfem_error("1D GetElementToEdgeTable is not yet implemented.");
6430  }
6431 
6432  // Return the number of edges
6433  return NumberOfEdges;
6434 }
6435 
6437 {
6438  if (el_to_el)
6439  {
6440  return *el_to_el;
6441  }
6442 
6443  // Note that, for ParNCMeshes, faces_info will contain also the ghost faces
6444  MFEM_ASSERT(faces_info.Size() >= GetNumFaces(), "faces were not generated!");
6445 
6446  Array<Connection> conn;
6447  conn.Reserve(2*faces_info.Size());
6448 
6449  for (int i = 0; i < faces_info.Size(); i++)
6450  {
6451  const FaceInfo &fi = faces_info[i];
6452  if (fi.Elem2No >= 0)
6453  {
6454  conn.Append(Connection(fi.Elem1No, fi.Elem2No));
6455  conn.Append(Connection(fi.Elem2No, fi.Elem1No));
6456  }
6457  else if (fi.Elem2Inf >= 0)
6458  {
6459  int nbr_elem_idx = NumOfElements - 1 - fi.Elem2No;
6460  conn.Append(Connection(fi.Elem1No, nbr_elem_idx));
6461  conn.Append(Connection(nbr_elem_idx, fi.Elem1No));
6462  }
6463  }
6464 
6465  conn.Sort();
6466  conn.Unique();
6467  el_to_el = new Table(NumOfElements, conn);
6468 
6469  return *el_to_el;
6470 }
6471 
6473 {
6474  if (el_to_face == NULL)
6475  {
6476  mfem_error("Mesh::ElementToFaceTable()");
6477  }
6478  return *el_to_face;
6479 }
6480 
6482 {
6483  if (el_to_edge == NULL)
6484  {
6485  mfem_error("Mesh::ElementToEdgeTable()");
6486  }
6487  return *el_to_edge;
6488 }
6489 
6490 void Mesh::AddPointFaceElement(int lf, int gf, int el)
6491 {
6492  if (faces_info[gf].Elem1No == -1) // this will be elem1
6493  {
6494  // faces[gf] = new Point(&gf);
6495  faces_info[gf].Elem1No = el;
6496  faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
6497  faces_info[gf].Elem2No = -1; // in case there's no other side
6498  faces_info[gf].Elem2Inf = -1; // face is not shared
6499  }
6500  else // this will be elem2
6501  {
6502  /* WARNING: Without the following check the mesh faces_info data structure
6503  may contain unreliable data. Normally, the order in which elements are
6504  processed could swap which elements appear as Elem1No and Elem2No. In
6505  branched meshes, where more than two elements can meet at a given node,
6506  the indices stored in Elem1No and Elem2No will be the first and last,
6507  respectively, elements found which touch a given node. This can lead to
6508  inconsistencies in any algorithms which rely on this data structure. To
6509  properly support branched meshes this data structure should be extended
6510  to support multiple elements per face. */
6511  /*
6512  MFEM_VERIFY(faces_info[gf].Elem2No < 0, "Invalid mesh topology. "
6513  "Interior point found connecting 1D elements "
6514  << faces_info[gf].Elem1No << ", " << faces_info[gf].Elem2No
6515  << " and " << el << ".");
6516  */
6517  faces_info[gf].Elem2No = el;
6518  faces_info[gf].Elem2Inf = 64 * lf + 1;
6519  }
6520 }
6521