mesh.cpp

Go to the documentation of this file.

// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.

// Implementation of data type mesh

#include "mesh_headers.hpp"
#include "../fem/fem.hpp"
#include "../general/sort_pairs.hpp"
#include "../general/binaryio.hpp"
#include "../general/text.hpp"
#include "../general/device.hpp"
#include "../general/tic_toc.hpp"
#include "../general/gecko.hpp"
#include "../fem/quadinterpolator.hpp"

#include <iostream>
#include <sstream>
#include <fstream>
#include <limits>
#include <cmath>
#include <cstring>
#include <ctime>
#include <functional>
#include <map>
#include <set>

// Include the METIS header, if using version 5. If using METIS 4, the needed
// declarations are inlined below, i.e. no header is needed.
#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
#include "metis.h"
#endif

// METIS 4 prototypes
#if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
typedef int idx_t;
typedef int idxtype;
extern "C" {
   void METIS_PartGraphRecursive(int*, idxtype*, idxtype*, idxtype*, idxtype*,
                                 int*, int*, int*, int*, int*, idxtype*);
   void METIS_PartGraphKway(int*, idxtype*, idxtype*, idxtype*, idxtype*,
                            int*, int*, int*, int*, int*, idxtype*);
   void METIS_PartGraphVKway(int*, idxtype*, idxtype*, idxtype*, idxtype*,
                             int*, int*, int*, int*, int*, idxtype*);
}
#endif

using namespace std;

namespace mfem
{

void Mesh::GetElementJacobian(int i, DenseMatrix &J)
{
   Geometry::Type geom = GetElementBaseGeometry(i);
   ElementTransformation *eltransf = GetElementTransformation(i);
   eltransf->SetIntPoint(&Geometries.GetCenter(geom));
   Geometries.JacToPerfJac(geom, eltransf->Jacobian(), J);
}

void Mesh::GetElementCenter(int i, Vector &center)
{
   center.SetSize(spaceDim);
   int geom = GetElementBaseGeometry(i);
   ElementTransformation *eltransf = GetElementTransformation(i);
   eltransf->Transform(Geometries.GetCenter(geom), center);
}

double Mesh::GetElementSize(ElementTransformation *T, int type)
{
   DenseMatrix J(Dim);

   Geometry::Type geom = T->GetGeometryType();
   T->SetIntPoint(&Geometries.GetCenter(geom));
   Geometries.JacToPerfJac(geom, T->Jacobian(), J);

   if (type == 0)
   {
      return pow(fabs(J.Det()), 1./Dim);
   }
   else if (type == 1)
   {
      return J.CalcSingularvalue(Dim-1);   // h_min
   }
   else
   {
      return J.CalcSingularvalue(0);   // h_max
   }
}

double Mesh::GetElementSize(int i, int type)
{
   return GetElementSize(GetElementTransformation(i), type);
}

double Mesh::GetElementSize(int i, const Vector &dir)
{
   DenseMatrix J(Dim);
   Vector d_hat(Dim);
   GetElementJacobian(i, J);
   J.MultTranspose(dir, d_hat);
   return sqrt((d_hat * d_hat) / (dir * dir));
}

double Mesh::GetElementVolume(int i)
{
   ElementTransformation *et = GetElementTransformation(i);
   const IntegrationRule &ir = IntRules.Get(GetElementBaseGeometry(i),
                                            et->OrderJ());
   double volume = 0.0;
   for (int j = 0; j < ir.GetNPoints(); j++)
   {
      const IntegrationPoint &ip = ir.IntPoint(j);
      et->SetIntPoint(&ip);
      volume += ip.weight * et->Weight();
   }

   return volume;
}

// Similar to VisualizationSceneSolution3d::FindNewBox in GLVis
void Mesh::GetBoundingBox(Vector &min, Vector &max, int ref)
{
   min.SetSize(spaceDim);
   max.SetSize(spaceDim);

   for (int d = 0; d < spaceDim; d++)
   {
      min(d) = infinity();
      max(d) = -infinity();
   }

   if (Nodes == NULL)
   {
      double *coord;
      for (int i = 0; i < NumOfVertices; i++)
      {
         coord = GetVertex(i);
         for (int d = 0; d < spaceDim; d++)
         {
            if (coord[d] < min(d)) { min(d) = coord[d]; }
            if (coord[d] > max(d)) { max(d) = coord[d]; }
         }
      }
   }
   else
   {
      const bool use_boundary = false; // make this a parameter?
      int ne = use_boundary ? GetNBE() : GetNE();
      int fn, fo;
      DenseMatrix pointmat;
      RefinedGeometry *RefG;
      IntegrationRule eir;
      FaceElementTransformations *Tr;
      ElementTransformation *T;

      for (int i = 0; i < ne; i++)
      {
         if (use_boundary)
         {
            GetBdrElementFace(i, &fn, &fo);
            RefG = GlobGeometryRefiner.Refine(GetFaceBaseGeometry(fn), ref);
            Tr = GetFaceElementTransformations(fn, 5);
            eir.SetSize(RefG->RefPts.GetNPoints());
            Tr->Loc1.Transform(RefG->RefPts, eir);
            Tr->Elem1->Transform(eir, pointmat);
         }
         else
         {
            T = GetElementTransformation(i);
            RefG = GlobGeometryRefiner.Refine(GetElementBaseGeometry(i), ref);
            T->Transform(RefG->RefPts, pointmat);
         }
         for (int j = 0; j < pointmat.Width(); j++)
         {
            for (int d = 0; d < pointmat.Height(); d++)
            {
               if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
               if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
            }
         }
      }
   }
}

void Mesh::GetCharacteristics(double &h_min, double &h_max,
                              double &kappa_min, double &kappa_max,
                              Vector *Vh, Vector *Vk)
{
   int i, dim, sdim;
   DenseMatrix J;
   double h, kappa;

   dim = Dimension();
   sdim = SpaceDimension();

   if (Vh) { Vh->SetSize(NumOfElements); }
   if (Vk) { Vk->SetSize(NumOfElements); }

   h_min = kappa_min = infinity();
   h_max = kappa_max = -h_min;
   if (dim == 0) { if (Vh) { *Vh = 1.0; } if (Vk) {*Vk = 1.0; } return; }
   J.SetSize(sdim, dim);
   for (i = 0; i < NumOfElements; i++)
   {
      GetElementJacobian(i, J);
      h = pow(fabs(J.Weight()), 1.0/double(dim));
      kappa = (dim == sdim) ?
              J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1) : -1.0;
      if (Vh) { (*Vh)(i) = h; }
      if (Vk) { (*Vk)(i) = kappa; }

      if (h < h_min) { h_min = h; }
      if (h > h_max) { h_max = h; }
      if (kappa < kappa_min) { kappa_min = kappa; }
      if (kappa > kappa_max) { kappa_max = kappa; }
   }
}

// static method
void Mesh::PrintElementsByGeometry(int dim,
                                   const Array<int> &num_elems_by_geom,
                                   std::ostream &out)
{
   for (int g = Geometry::DimStart[dim], first = 1;
        g < Geometry::DimStart[dim+1]; g++)
   {
      if (!num_elems_by_geom[g]) { continue; }
      if (!first) { out << " + "; }
      else { first = 0; }
      out << num_elems_by_geom[g] << ' ' << Geometry::Name[g] << "(s)";
   }
}

void Mesh::PrintCharacteristics(Vector *Vh, Vector *Vk, std::ostream &out)
{
   double h_min, h_max, kappa_min, kappa_max;

   out << "Mesh Characteristics:";

   this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);

   Array<int> num_elems_by_geom(Geometry::NumGeom);
   num_elems_by_geom = 0;
   for (int i = 0; i < GetNE(); i++)
   {
      num_elems_by_geom[GetElementBaseGeometry(i)]++;
   }

   out << '\n'
       << "Dimension          : " << Dimension() << '\n'
       << "Space dimension    : " << SpaceDimension();
   if (Dim == 0)
   {
      out << '\n'
          << "Number of vertices : " << GetNV() << '\n'
          << "Number of elements : " << GetNE() << '\n'
          << "Number of bdr elem : " << GetNBE() << '\n';
   }
   else if (Dim == 1)
   {
      out << '\n'
          << "Number of vertices : " << GetNV() << '\n'
          << "Number of elements : " << GetNE() << '\n'
          << "Number of bdr elem : " << GetNBE() << '\n'
          << "h_min              : " << h_min << '\n'
          << "h_max              : " << h_max << '\n';
   }
   else if (Dim == 2)
   {
      out << '\n'
          << "Number of vertices : " << GetNV() << '\n'
          << "Number of edges    : " << GetNEdges() << '\n'
          << "Number of elements : " << GetNE() << "  --  ";
      PrintElementsByGeometry(2, num_elems_by_geom, out);
      out << '\n'
          << "Number of bdr elem : " << GetNBE() << '\n'
          << "Euler Number       : " << EulerNumber2D() << '\n'
          << "h_min              : " << h_min << '\n'
          << "h_max              : " << h_max << '\n'
          << "kappa_min          : " << kappa_min << '\n'
          << "kappa_max          : " << kappa_max << '\n';
   }
   else
   {
      Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
      num_bdr_elems_by_geom = 0;
      for (int i = 0; i < GetNBE(); i++)
      {
         num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
      }
      Array<int> num_faces_by_geom(Geometry::NumGeom);
      num_faces_by_geom = 0;
      for (int i = 0; i < GetNFaces(); i++)
      {
         num_faces_by_geom[GetFaceBaseGeometry(i)]++;
      }

      out << '\n'
          << "Number of vertices : " << GetNV() << '\n'
          << "Number of edges    : " << GetNEdges() << '\n'
          << "Number of faces    : " << GetNFaces() << "  --  ";
      PrintElementsByGeometry(Dim-1, num_faces_by_geom, out);
      out << '\n'
          << "Number of elements : " << GetNE() << "  --  ";
      PrintElementsByGeometry(Dim, num_elems_by_geom, out);
      out << '\n'
          << "Number of bdr elem : " << GetNBE() << "  --  ";
      PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, out);
      out << '\n'
          << "Euler Number       : " << EulerNumber() << '\n'
          << "h_min              : " << h_min << '\n'
          << "h_max              : " << h_max << '\n'
          << "kappa_min          : " << kappa_min << '\n'
          << "kappa_max          : " << kappa_max << '\n';
   }
   out << '\n' << std::flush;
}

FiniteElement *Mesh::GetTransformationFEforElementType(Element::Type ElemType)
{
   switch (ElemType)
   {
      case Element::POINT :          return &PointFE;
      case Element::SEGMENT :        return &SegmentFE;
      case Element::TRIANGLE :       return &TriangleFE;
      case Element::QUADRILATERAL :  return &QuadrilateralFE;
      case Element::TETRAHEDRON :    return &TetrahedronFE;
      case Element::HEXAHEDRON :     return &HexahedronFE;
      case Element::WEDGE :          return &WedgeFE;
      default:
         MFEM_ABORT("Unknown element type \"" << ElemType << "\"");
         break;
   }
   MFEM_ABORT("Unknown element type");
   return NULL;
}


void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
{
   ElTr->Attribute = GetAttribute(i);
   ElTr->ElementNo = i;
   ElTr->ElementType = ElementTransformation::ELEMENT;
   ElTr->Reset();
   if (Nodes == NULL)
   {
      GetPointMatrix(i, ElTr->GetPointMat());
      ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
   }
   else
   {
      DenseMatrix &pm = ElTr->GetPointMat();
      Array<int> vdofs;
      Nodes->FESpace()->GetElementVDofs(i, vdofs);
      Nodes->HostRead();
      const GridFunction &nodes = *Nodes;
      int n = vdofs.Size()/spaceDim;
      pm.SetSize(spaceDim, n);
      for (int k = 0; k < spaceDim; k++)
      {
         for (int j = 0; j < n; j++)
         {
            pm(k,j) = nodes(vdofs[n*k+j]);
         }
      }
      ElTr->SetFE(Nodes->FESpace()->GetFE(i));
   }
}

void Mesh::GetElementTransformation(int i, const Vector &nodes,
                                    IsoparametricTransformation *ElTr)
{
   ElTr->Attribute = GetAttribute(i);
   ElTr->ElementNo = i;
   ElTr->ElementType = ElementTransformation::ELEMENT;
   DenseMatrix &pm = ElTr->GetPointMat();
   ElTr->Reset();
   nodes.HostRead();
   if (Nodes == NULL)
   {
      MFEM_ASSERT(nodes.Size() == spaceDim*GetNV(), "");
      int       nv = elements[i]->GetNVertices();
      const int *v = elements[i]->GetVertices();
      int n = vertices.Size();
      pm.SetSize(spaceDim, nv);
      for (int k = 0; k < spaceDim; k++)
      {
         for (int j = 0; j < nv; j++)
         {
            pm(k, j) = nodes(k*n+v[j]);
         }
      }
      ElTr->SetFE(GetTransformationFEforElementType(GetElementType(i)));
   }
   else
   {
      MFEM_ASSERT(nodes.Size() == Nodes->Size(), "");
      Array<int> vdofs;
      Nodes->FESpace()->GetElementVDofs(i, vdofs);
      int n = vdofs.Size()/spaceDim;
      pm.SetSize(spaceDim, n);
      for (int k = 0; k < spaceDim; k++)
      {
         for (int j = 0; j < n; j++)
         {
            pm(k,j) = nodes(vdofs[n*k+j]);
         }
      }
      ElTr->SetFE(Nodes->FESpace()->GetFE(i));
   }
}

ElementTransformation *Mesh::GetElementTransformation(int i)
{
   GetElementTransformation(i, &Transformation);

   return &Transformation;
}

ElementTransformation *Mesh::GetBdrElementTransformation(int i)
{
   GetBdrElementTransformation(i, &BdrTransformation);
   return &BdrTransformation;
}

void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
{
   ElTr->Attribute = GetBdrAttribute(i);
   ElTr->ElementNo = i; // boundary element number
   ElTr->ElementType = ElementTransformation::BDR_ELEMENT;
   DenseMatrix &pm = ElTr->GetPointMat();
   ElTr->Reset();
   if (Nodes == NULL)
   {
      GetBdrPointMatrix(i, pm);
      ElTr->SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
   }
   else
   {
      const FiniteElement *bdr_el = Nodes->FESpace()->GetBE(i);
      Nodes->HostRead();
      const GridFunction &nodes = *Nodes;
      if (bdr_el)
      {
         Array<int> vdofs;
         Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
         int n = vdofs.Size()/spaceDim;
         pm.SetSize(spaceDim, n);
         for (int k = 0; k < spaceDim; k++)
         {
            for (int j = 0; j < n; j++)
            {
               pm(k,j) = nodes(vdofs[n*k+j]);
            }
         }
         ElTr->SetFE(bdr_el);
      }
      else // L2 Nodes (e.g., periodic mesh)
      {
         int elem_id, face_info;
         GetBdrElementAdjacentElement(i, elem_id, face_info);

         GetLocalFaceTransformation(GetBdrElementType(i),
                                    GetElementType(elem_id),
                                    FaceElemTr.Loc1.Transf, face_info);
         // NOTE: FaceElemTr.Loc1 is overwritten here -- used as a temporary

         Geometry::Type face_geom = GetBdrElementBaseGeometry(i);
         const FiniteElement *face_el =
            Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
         MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
                     "Mesh requires nodal Finite Element.");
         IntegrationRule eir(face_el->GetDof());
         FaceElemTr.Loc1.Transf.ElementNo = elem_id;
         FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
         FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
         Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);

         ElTr->SetFE(face_el);
      }
   }
}

void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
{
   FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
   FTr->ElementNo = FaceNo;
   FTr->ElementType = ElementTransformation::FACE;
   DenseMatrix &pm = FTr->GetPointMat();
   FTr->Reset();
   if (Nodes == NULL)
   {
      const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
      const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
      pm.SetSize(spaceDim, nv);
      for (int i = 0; i < spaceDim; i++)
      {
         for (int j = 0; j < nv; j++)
         {
            pm(i, j) = vertices[v[j]](i);
         }
      }
      FTr->SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
   }
   else // curved mesh
   {
      const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
      Nodes->HostRead();
      const GridFunction &nodes = *Nodes;
      if (face_el)
      {
         Array<int> vdofs;
         Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
         int n = vdofs.Size()/spaceDim;
         pm.SetSize(spaceDim, n);
         for (int i = 0; i < spaceDim; i++)
         {
            for (int j = 0; j < n; j++)
            {
               pm(i, j) = nodes(vdofs[n*i+j]);
            }
         }
         FTr->SetFE(face_el);
      }
      else // L2 Nodes (e.g., periodic mesh), go through the volume of Elem1
      {
         FaceInfo &face_info = faces_info[FaceNo];

         Geometry::Type face_geom = GetFaceGeometryType(FaceNo);
         Element::Type  face_type = GetFaceElementType(FaceNo);

         GetLocalFaceTransformation(face_type,
                                    GetElementType(face_info.Elem1No),
                                    FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
         // NOTE: FaceElemTr.Loc1 is overwritten here -- used as a temporary

         face_el = Nodes->FESpace()->GetTraceElement(face_info.Elem1No,
                                                     face_geom);
         MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
                     "Mesh requires nodal Finite Element.");

         IntegrationRule eir(face_el->GetDof());
         FaceElemTr.Loc1.Transf.ElementNo = face_info.Elem1No;
         FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
         FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
         Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);

         FTr->SetFE(face_el);
      }
   }
}

ElementTransformation *Mesh::GetFaceTransformation(int FaceNo)
{
   GetFaceTransformation(FaceNo, &FaceTransformation);
   return &FaceTransformation;
}

void Mesh::GetEdgeTransformation(int EdgeNo, IsoparametricTransformation *EdTr)
{
   if (Dim == 2)
   {
      GetFaceTransformation(EdgeNo, EdTr);
      return;
   }
   if (Dim == 1)
   {
      mfem_error("Mesh::GetEdgeTransformation not defined in 1D \n");
   }

   EdTr->Attribute = 1;
   EdTr->ElementNo = EdgeNo;
   EdTr->ElementType = ElementTransformation::EDGE;
   DenseMatrix &pm = EdTr->GetPointMat();
   EdTr->Reset();
   if (Nodes == NULL)
   {
      Array<int> v;
      GetEdgeVertices(EdgeNo, v);
      const int nv = 2;
      pm.SetSize(spaceDim, nv);
      for (int i = 0; i < spaceDim; i++)
      {
         for (int j = 0; j < nv; j++)
         {
            pm(i, j) = vertices[v[j]](i);
         }
      }
      EdTr->SetFE(GetTransformationFEforElementType(Element::SEGMENT));
   }
   else
   {
      const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
      Nodes->HostRead();
      const GridFunction &nodes = *Nodes;
      if (edge_el)
      {
         Array<int> vdofs;
         Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
         int n = vdofs.Size()/spaceDim;
         pm.SetSize(spaceDim, n);
         for (int i = 0; i < spaceDim; i++)
         {
            for (int j = 0; j < n; j++)
            {
               pm(i, j) = nodes(vdofs[n*i+j]);
            }
         }
         EdTr->SetFE(edge_el);
      }
      else
      {
         MFEM_ABORT("Not implemented.");
      }
   }
}

ElementTransformation *Mesh::GetEdgeTransformation(int EdgeNo)
{
   GetEdgeTransformation(EdgeNo, &EdgeTransformation);
   return &EdgeTransformation;
}


void Mesh::GetLocalPtToSegTransformation(
   IsoparametricTransformation &Transf, int i)
{
   const IntegrationRule *SegVert;
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&PointFE);
   SegVert = Geometries.GetVertices(Geometry::SEGMENT);
   locpm.SetSize(1, 1);
   locpm(0, 0) = SegVert->IntPoint(i/64).x;
   //  (i/64) is the local face no. in the segment
   //  (i%64) is the orientation of the point (not used)
}

void Mesh::GetLocalSegToTriTransformation(
   IsoparametricTransformation &Transf, int i)
{
   const int *tv, *so;
   const IntegrationRule *TriVert;
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&SegmentFE);
   tv = tri_t::Edges[i/64];  //  (i/64) is the local face no. in the triangle
   so = seg_t::Orient[i%64]; //  (i%64) is the orientation of the segment
   TriVert = Geometries.GetVertices(Geometry::TRIANGLE);
   locpm.SetSize(2, 2);
   for (int j = 0; j < 2; j++)
   {
      locpm(0, so[j]) = TriVert->IntPoint(tv[j]).x;
      locpm(1, so[j]) = TriVert->IntPoint(tv[j]).y;
   }
}

void Mesh::GetLocalSegToQuadTransformation(
   IsoparametricTransformation &Transf, int i)
{
   const int *qv, *so;
   const IntegrationRule *QuadVert;
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&SegmentFE);
   qv = quad_t::Edges[i/64]; //  (i/64) is the local face no. in the quad
   so = seg_t::Orient[i%64]; //  (i%64) is the orientation of the segment
   QuadVert = Geometries.GetVertices(Geometry::SQUARE);
   locpm.SetSize(2, 2);
   for (int j = 0; j < 2; j++)
   {
      locpm(0, so[j]) = QuadVert->IntPoint(qv[j]).x;
      locpm(1, so[j]) = QuadVert->IntPoint(qv[j]).y;
   }
}

void Mesh::GetLocalTriToTetTransformation(
   IsoparametricTransformation &Transf, int i)
{
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&TriangleFE);
   //  (i/64) is the local face no. in the tet
   const int *tv = tet_t::FaceVert[i/64];
   //  (i%64) is the orientation of the tetrahedron face
   //         w.r.t. the face element
   const int *to = tri_t::Orient[i%64];
   const IntegrationRule *TetVert =
      Geometries.GetVertices(Geometry::TETRAHEDRON);
   locpm.SetSize(3, 3);
   for (int j = 0; j < 3; j++)
   {
      const IntegrationPoint &vert = TetVert->IntPoint(tv[to[j]]);
      locpm(0, j) = vert.x;
      locpm(1, j) = vert.y;
      locpm(2, j) = vert.z;
   }
}

void Mesh::GetLocalTriToWdgTransformation(
   IsoparametricTransformation &Transf, int i)
{
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&TriangleFE);
   //  (i/64) is the local face no. in the pri
   MFEM_VERIFY(i < 128, "Local face index " << i/64
               << " is not a triangular face of a wedge.");
   const int *pv = pri_t::FaceVert[i/64];
   //  (i%64) is the orientation of the wedge face
   //         w.r.t. the face element
   const int *to = tri_t::Orient[i%64];
   const IntegrationRule *PriVert =
      Geometries.GetVertices(Geometry::PRISM);
   locpm.SetSize(3, 3);
   for (int j = 0; j < 3; j++)
   {
      const IntegrationPoint &vert = PriVert->IntPoint(pv[to[j]]);
      locpm(0, j) = vert.x;
      locpm(1, j) = vert.y;
      locpm(2, j) = vert.z;
   }
}

void Mesh::GetLocalQuadToHexTransformation(
   IsoparametricTransformation &Transf, int i)
{
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&QuadrilateralFE);
   //  (i/64) is the local face no. in the hex
   const int *hv = hex_t::FaceVert[i/64];
   //  (i%64) is the orientation of the quad
   const int *qo = quad_t::Orient[i%64];
   const IntegrationRule *HexVert = Geometries.GetVertices(Geometry::CUBE);
   locpm.SetSize(3, 4);
   for (int j = 0; j < 4; j++)
   {
      const IntegrationPoint &vert = HexVert->IntPoint(hv[qo[j]]);
      locpm(0, j) = vert.x;
      locpm(1, j) = vert.y;
      locpm(2, j) = vert.z;
   }
}

void Mesh::GetLocalQuadToWdgTransformation(
   IsoparametricTransformation &Transf, int i)
{
   DenseMatrix &locpm = Transf.GetPointMat();
   Transf.Reset();

   Transf.SetFE(&QuadrilateralFE);
   //  (i/64) is the local face no. in the pri
   MFEM_VERIFY(i >= 128, "Local face index " << i/64
               << " is not a quadrilateral face of a wedge.");
   const int *pv = pri_t::FaceVert[i/64];
   //  (i%64) is the orientation of the quad
   const int *qo = quad_t::Orient[i%64];
   const IntegrationRule *PriVert = Geometries.GetVertices(Geometry::PRISM);
   locpm.SetSize(3, 4);
   for (int j = 0; j < 4; j++)
   {
      const IntegrationPoint &vert = PriVert->IntPoint(pv[qo[j]]);
      locpm(0, j) = vert.x;
      locpm(1, j) = vert.y;
      locpm(2, j) = vert.z;
   }
}

const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
                                                  const int flags,
                                                  MemoryType d_mt)
{
   for (int i = 0; i < geom_factors.Size(); i++)
   {
      GeometricFactors *gf = geom_factors[i];
      if (gf->IntRule == &ir && (gf->computed_factors & flags) == flags)
      {
         return gf;
      }
   }

   this->EnsureNodes();

   GeometricFactors *gf = new GeometricFactors(this, ir, flags, d_mt);
   geom_factors.Append(gf);
   return gf;
}

const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
   const IntegrationRule& ir,
   const int flags, FaceType type)
{
   for (int i = 0; i < face_geom_factors.Size(); i++)
   {
      FaceGeometricFactors *gf = face_geom_factors[i];
      if (gf->IntRule == &ir && (gf->computed_factors & flags) == flags &&
          gf->type==type)
      {
         return gf;
      }
   }

   this->EnsureNodes();

   FaceGeometricFactors *gf = new FaceGeometricFactors(this, ir, flags, type);
   face_geom_factors.Append(gf);
   return gf;
}

void Mesh::DeleteGeometricFactors()
{
   for (int i = 0; i < geom_factors.Size(); i++)
   {
      delete geom_factors[i];
   }
   geom_factors.SetSize(0);
   for (int i = 0; i < face_geom_factors.Size(); i++)
   {
      delete face_geom_factors[i];
   }
   face_geom_factors.SetSize(0);
}

void Mesh::GetLocalFaceTransformation(
   int face_type, int elem_type, IsoparametricTransformation &Transf, int info)
{
   switch (face_type)
   {
      case Element::POINT:
         GetLocalPtToSegTransformation(Transf, info);
         break;

      case Element::SEGMENT:
         if (elem_type == Element::TRIANGLE)
         {
            GetLocalSegToTriTransformation(Transf, info);
         }
         else
         {
            MFEM_ASSERT(elem_type == Element::QUADRILATERAL, "");
            GetLocalSegToQuadTransformation(Transf, info);
         }
         break;

      case Element::TRIANGLE:
         if (elem_type == Element::TETRAHEDRON)
         {
            GetLocalTriToTetTransformation(Transf, info);
         }
         else
         {
            MFEM_ASSERT(elem_type == Element::WEDGE, "");
            GetLocalTriToWdgTransformation(Transf, info);
         }
         break;

      case Element::QUADRILATERAL:
         if (elem_type == Element::HEXAHEDRON)
         {
            GetLocalQuadToHexTransformation(Transf, info);
         }
         else
         {
            MFEM_ASSERT(elem_type == Element::WEDGE, "");
            GetLocalQuadToWdgTransformation(Transf, info);
         }
         break;
   }
}

FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
                                                                int mask)
{
   FaceInfo &face_info = faces_info[FaceNo];

   int cmask = 0;
   FaceElemTr.SetConfigurationMask(cmask);
   FaceElemTr.Elem1 = NULL;
   FaceElemTr.Elem2 = NULL;

   // setup the transformation for the first element
   FaceElemTr.Elem1No = face_info.Elem1No;
   if (mask & FaceElementTransformations::HAVE_ELEM1)
   {
      GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
      FaceElemTr.Elem1 = &Transformation;
      cmask |= 1;
   }

   //  setup the transformation for the second element
   //     return NULL in the Elem2 field if there's no second element, i.e.
   //     the face is on the "boundary"
   FaceElemTr.Elem2No = face_info.Elem2No;
   if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
       FaceElemTr.Elem2No >= 0)
   {
#ifdef MFEM_DEBUG
      if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
      { MFEM_ABORT("NURBS mesh not supported!"); }
#endif
      GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
      FaceElemTr.Elem2 = &Transformation2;
      cmask |= 2;
   }

   // setup the face transformation
   if (mask & FaceElementTransformations::HAVE_FACE)
   {
      GetFaceTransformation(FaceNo, &FaceElemTr);
      cmask |= 16;
   }
   else
   {
      FaceElemTr.SetGeometryType(GetFaceGeometryType(FaceNo));
   }

   // setup Loc1 & Loc2
   int face_type = GetFaceElementType(FaceNo);
   if (mask & FaceElementTransformations::HAVE_LOC1)
   {
      int elem_type = GetElementType(face_info.Elem1No);
      GetLocalFaceTransformation(face_type, elem_type,
                                 FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
      cmask |= 4;
   }
   if ((mask & FaceElementTransformations::HAVE_LOC2) &&
       FaceElemTr.Elem2No >= 0)
   {
      int elem_type = GetElementType(face_info.Elem2No);
      GetLocalFaceTransformation(face_type, elem_type,
                                 FaceElemTr.Loc2.Transf, face_info.Elem2Inf);

      // NC meshes: prepend slave edge/face transformation to Loc2
      if (Nonconforming() && IsSlaveFace(face_info))
      {
         ApplyLocalSlaveTransformation(FaceElemTr, face_info, false);
      }
      cmask |= 8;
   }

   FaceElemTr.SetConfigurationMask(cmask);

   // This check can be useful for internal debugging, however it will fail on
   // periodic boundary faces, so we keep it disabled in general.
#if 0
#ifdef MFEM_DEBUG
   double dist = FaceElemTr.CheckConsistency();
   if (dist >= 1e-12)
   {
      mfem::out << "\nInternal error: face id = " << FaceNo
                << ", dist = " << dist << '\n';
      FaceElemTr.CheckConsistency(1); // print coordinates
      MFEM_ABORT("internal error");
   }
#endif
#endif

   return &FaceElemTr;
}

bool Mesh::IsSlaveFace(const FaceInfo &fi) const
{
   return fi.NCFace >= 0 && nc_faces_info[fi.NCFace].Slave;
}

void Mesh::ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
                                         const FaceInfo &fi, bool is_ghost)
{
#ifdef MFEM_THREAD_SAFE
   DenseMatrix composition;
#else
   static DenseMatrix composition;
#endif
   MFEM_ASSERT(fi.NCFace >= 0, "");
   MFEM_ASSERT(nc_faces_info[fi.NCFace].Slave, "internal error");
   if (!is_ghost)
   {
      // side 1 -> child side, side 2 -> parent side
      IsoparametricTransformation &LT = FT.Loc2.Transf;
      LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
      // In 2D, we need to flip the point matrix since it is aligned with the
      // parent side.
      if (Dim == 2)
      {
         // swap points (columns) 0 and 1
         std::swap(composition(0,0), composition(0,1));
         std::swap(composition(1,0), composition(1,1));
      }
      LT.SetPointMat(composition);
   }
   else // is_ghost == true
   {
      // side 1 -> parent side, side 2 -> child side
      IsoparametricTransformation &LT = FT.Loc1.Transf;
      LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
      // In 2D, there is no need to flip the point matrix since it is already
      // aligned with the parent side, see also ParNCMesh::GetFaceNeighbors.
      // In 3D the point matrix was flipped during construction in
      // ParNCMesh::GetFaceNeighbors and due to that it is already aligned with
      // the parent side.
      LT.SetPointMat(composition);
   }
}

FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
{
   FaceElementTransformations *tr;
   int fn = GetBdrFace(BdrElemNo);

   // Check if the face is interior, shared, or non-conforming.
   if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
   {
      return NULL;
   }
   tr = GetFaceElementTransformations(fn, 21);
   tr->Attribute = boundary[BdrElemNo]->GetAttribute();
   tr->ElementNo = BdrElemNo;
   tr->ElementType = ElementTransformation::BDR_FACE;
   return tr;
}

int Mesh::GetBdrFace(int BdrElemNo) const
{
   int fn;
   if (Dim == 3)
   {
      fn = be_to_face[BdrElemNo];
   }
   else if (Dim == 2)
   {
      fn = be_to_edge[BdrElemNo];
   }
   else
   {
      fn = boundary[BdrElemNo]->GetVertices()[0];
   }
   return fn;
}

void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2) const
{
   *Elem1 = faces_info[Face].Elem1No;
   *Elem2 = faces_info[Face].Elem2No;
}

void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2) const
{
   *Inf1 = faces_info[Face].Elem1Inf;
   *Inf2 = faces_info[Face].Elem2Inf;
}

void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2, int *NCFace) const
{
   *Inf1   = faces_info[Face].Elem1Inf;
   *Inf2   = faces_info[Face].Elem2Inf;
   *NCFace = faces_info[Face].NCFace;
}

Geometry::Type Mesh::GetFaceGeometryType(int Face) const
{
   switch (Dim)
   {
      case 1: return Geometry::POINT;
      case 2: return Geometry::SEGMENT;
      case 3:
         if (Face < NumOfFaces) // local (non-ghost) face
         {
            return faces[Face]->GetGeometryType();
         }
         // ghost face
         const int nc_face_id = faces_info[Face].NCFace;
         MFEM_ASSERT(nc_face_id >= 0, "parent ghost faces are not supported");
         return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
   }
   return Geometry::INVALID;
}

Element::Type Mesh::GetFaceElementType(int Face) const
{
   return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
}

void Mesh::Init()
{
   // in order of declaration:
   Dim = spaceDim = 0;
   NumOfVertices = -1;
   NumOfElements = NumOfBdrElements = 0;
   NumOfEdges = NumOfFaces = 0;
   nbInteriorFaces = -1;
   nbBoundaryFaces = -1;
   meshgen = mesh_geoms = 0;
   sequence = 0;
   Nodes = NULL;
   own_nodes = 1;
   NURBSext = NULL;
   ncmesh = NULL;
   last_operation = Mesh::NONE;
}

void Mesh::InitTables()
{
   el_to_edge =
      el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
}

void Mesh::SetEmpty()
{
   Init();
   InitTables();
}

void Mesh::DestroyTables()
{
   delete el_to_edge;
   delete el_to_face;
   delete el_to_el;
   DeleteGeometricFactors();

   if (Dim == 3)
   {
      delete bel_to_edge;
   }

   delete face_edge;
   delete edge_vertex;
}

void Mesh::DestroyPointers()
{
   if (own_nodes) { delete Nodes; }

   delete ncmesh;

   delete NURBSext;

   for (int i = 0; i < NumOfElements; i++)
   {
      FreeElement(elements[i]);
   }

   for (int i = 0; i < NumOfBdrElements; i++)
   {
      FreeElement(boundary[i]);
   }

   for (int i = 0; i < faces.Size(); i++)
   {
      FreeElement(faces[i]);
   }

   DestroyTables();
}

void Mesh::Destroy()
{
   DestroyPointers();

   elements.DeleteAll();
   vertices.DeleteAll();
   boundary.DeleteAll();
   faces.DeleteAll();
   faces_info.DeleteAll();
   nc_faces_info.DeleteAll();
   be_to_edge.DeleteAll();
   be_to_face.DeleteAll();

   // TODO:
   // IsoparametricTransformations
   // Transformation, Transformation2, BdrTransformation, FaceTransformation,
   // EdgeTransformation;
   // FaceElementTransformations FaceElemTr;

   CoarseFineTr.Clear();

#ifdef MFEM_USE_MEMALLOC
   TetMemory.Clear();
#endif

   attributes.DeleteAll();
   bdr_attributes.DeleteAll();
}

void Mesh::ResetLazyData()
{
   delete el_to_el;     el_to_el = NULL;
   delete face_edge;    face_edge = NULL;
   delete edge_vertex;  edge_vertex = NULL;
   DeleteGeometricFactors();
   nbInteriorFaces = -1;
   nbBoundaryFaces = -1;
}

void Mesh::SetAttributes()
{
   Array<int> attribs;

   attribs.SetSize(GetNBE());
   for (int i = 0; i < attribs.Size(); i++)
   {
      attribs[i] = GetBdrAttribute(i);
   }
   attribs.Sort();
   attribs.Unique();
   attribs.Copy(bdr_attributes);
   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
   {
      MFEM_WARNING("Non-positive attributes on the boundary!");
   }

   attribs.SetSize(GetNE());
   for (int i = 0; i < attribs.Size(); i++)
   {
      attribs[i] = GetAttribute(i);
   }
   attribs.Sort();
   attribs.Unique();
   attribs.Copy(attributes);
   if (attributes.Size() > 0 && attributes[0] <= 0)
   {
      MFEM_WARNING("Non-positive attributes in the domain!");
   }
}

void Mesh::InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
{
   SetEmpty();

   Dim = Dim_;
   spaceDim = spaceDim_;

   NumOfVertices = 0;
   vertices.SetSize(NVert);  // just allocate space for vertices

   NumOfElements = 0;
   elements.SetSize(NElem);  // just allocate space for Element *

   NumOfBdrElements = 0;
   boundary.SetSize(NBdrElem);  // just allocate space for Element *
}

template<typename T>
static void CheckEnlarge(Array<T> &array, int size)
{
   if (size >= array.Size()) { array.SetSize(size + 1); }
}

int Mesh::AddVertex(double x, double y, double z)
{
   CheckEnlarge(vertices, NumOfVertices);
   double *v = vertices[NumOfVertices]();
   v[0] = x;
   v[1] = y;
   v[2] = z;
   return NumOfVertices++;
}

int Mesh::AddVertex(const double *coords)
{
   CheckEnlarge(vertices, NumOfVertices);
   vertices[NumOfVertices].SetCoords(spaceDim, coords);
   return NumOfVertices++;
}

void Mesh::AddVertexParents(int i, int p1, int p2)
{
   tmp_vertex_parents.Append(Triple<int, int, int>(i, p1, p2));

   // if vertex coordinates are defined, make sure the hanging vertex has the
   // correct position
   if (i < vertices.Size())
   {
      double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
      for (int j = 0; j < 3; j++)
      {
         vi[j] = (vp1[j] + vp2[j]) * 0.5;
      }
   }
}

int Mesh::AddSegment(int v1, int v2, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Segment(v1, v2, attr);
   return NumOfElements++;
}

int Mesh::AddSegment(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Segment(vi, attr);
   return NumOfElements++;
}

int Mesh::AddTriangle(int v1, int v2, int v3, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Triangle(v1, v2, v3, attr);
   return NumOfElements++;
}

int Mesh::AddTriangle(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Triangle(vi, attr);
   return NumOfElements++;
}

int Mesh::AddQuad(int v1, int v2, int v3, int v4, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Quadrilateral(v1, v2, v3, v4, attr);
   return NumOfElements++;
}

int Mesh::AddQuad(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Quadrilateral(vi, attr);
   return NumOfElements++;
}

int Mesh::AddTet(int v1, int v2, int v3, int v4, int attr)
{
   int vi[4] = {v1, v2, v3, v4};
   return AddTet(vi, attr);
}

int Mesh::AddTet(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
#ifdef MFEM_USE_MEMALLOC
   Tetrahedron *tet;
   tet = TetMemory.Alloc();
   tet->SetVertices(vi);
   tet->SetAttribute(attr);
   elements[NumOfElements] = tet;
#else
   elements[NumOfElements] = new Tetrahedron(vi, attr);
#endif
   return NumOfElements++;
}

int Mesh::AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Wedge(v1, v2, v3, v4, v5, v6, attr);
   return NumOfElements++;
}

int Mesh::AddWedge(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Wedge(vi, attr);
   return NumOfElements++;
}

int Mesh::AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8,
                 int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] =
      new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
   return NumOfElements++;
}

int Mesh::AddHex(const int *vi, int attr)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = new Hexahedron(vi, attr);
   return NumOfElements++;
}

void Mesh::AddHexAsTets(const int *vi, int attr)
{
   static const int hex_to_tet[6][4] =
   {
      { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
      { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
   };
   int ti[4];

   for (int i = 0; i < 6; i++)
   {
      for (int j = 0; j < 4; j++)
      {
         ti[j] = vi[hex_to_tet[i][j]];
      }
      AddTet(ti, attr);
   }
}

void Mesh::AddHexAsWedges(const int *vi, int attr)
{
   static const int hex_to_wdg[2][6] =
   {
      { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
   };
   int ti[6];

   for (int i = 0; i < 2; i++)
   {
      for (int j = 0; j < 6; j++)
      {
         ti[j] = vi[hex_to_wdg[i][j]];
      }
      AddWedge(ti, attr);
   }
}

int Mesh::AddElement(Element *elem)
{
   CheckEnlarge(elements, NumOfElements);
   elements[NumOfElements] = elem;
   return NumOfElements++;
}

int Mesh::AddBdrElement(Element *elem)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = elem;
   return NumOfBdrElements++;
}

int Mesh::AddBdrSegment(int v1, int v2, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Segment(v1, v2, attr);
   return NumOfBdrElements++;
}

int Mesh::AddBdrSegment(const int *vi, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Segment(vi, attr);
   return NumOfBdrElements++;
}

int Mesh::AddBdrTriangle(int v1, int v2, int v3, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Triangle(v1, v2, v3, attr);
   return NumOfBdrElements++;
}

int Mesh::AddBdrTriangle(const int *vi, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Triangle(vi, attr);
   return NumOfBdrElements++;
}

int Mesh::AddBdrQuad(int v1, int v2, int v3, int v4, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Quadrilateral(v1, v2, v3, v4, attr);
   return NumOfBdrElements++;
}

int Mesh::AddBdrQuad(const int *vi, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Quadrilateral(vi, attr);
   return NumOfBdrElements++;
}

void Mesh::AddBdrQuadAsTriangles(const int *vi, int attr)
{
   static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
   int ti[3];

   for (int i = 0; i < 2; i++)
   {
      for (int j = 0; j < 3; j++)
      {
         ti[j] = vi[quad_to_tri[i][j]];
      }
      AddBdrTriangle(ti, attr);
   }
}

int Mesh::AddBdrPoint(int v, int attr)
{
   CheckEnlarge(boundary, NumOfBdrElements);
   boundary[NumOfBdrElements] = new Point(&v, attr);
   return NumOfBdrElements++;
}

void Mesh::GenerateBoundaryElements()
{
   int i, j;
   Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;

   // GenerateFaces();

   for (i = 0; i < boundary.Size(); i++)
   {
      FreeElement(boundary[i]);
   }

   if (Dim == 3)
   {
      delete bel_to_edge;
      bel_to_edge = NULL;
   }

   // count the 'NumOfBdrElements'
   NumOfBdrElements = 0;
   for (i = 0; i < faces_info.Size(); i++)
   {
      if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
   }

   boundary.SetSize(NumOfBdrElements);
   be2face.SetSize(NumOfBdrElements);
   for (j = i = 0; i < faces_info.Size(); i++)
   {
      if (faces_info[i].Elem2No < 0)
      {
         boundary[j] = faces[i]->Duplicate(this);
         be2face[j++] = i;
      }
   }
   // In 3D, 'bel_to_edge' is destroyed but it's not updated.
}

void Mesh::FinalizeCheck()
{
   MFEM_VERIFY(vertices.Size() == NumOfVertices ||
               vertices.Size() == 0,
               "incorrect number of vertices: preallocated: " << vertices.Size()
               << ", actually added: " << NumOfVertices);
   MFEM_VERIFY(elements.Size() == NumOfElements,
               "incorrect number of elements: preallocated: " << elements.Size()
               << ", actually added: " << NumOfElements);
   MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
               "incorrect number of boundary elements: preallocated: "
               << boundary.Size() << ", actually added: " << NumOfBdrElements);
}

void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
{
   FinalizeCheck();
   CheckElementOrientation(fix_orientation);

   if (refine)
   {
      MarkTriMeshForRefinement();
   }

   if (generate_edges)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      GenerateFaces();
      CheckBdrElementOrientation();
   }
   else
   {
      NumOfEdges = 0;
   }

   NumOfFaces = 0;

   SetAttributes();

   SetMeshGen();
}

void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
                            bool fix_orientation)
{
   FinalizeCheck();
   if (fix_orientation)
   {
      CheckElementOrientation(fix_orientation);
   }

   if (generate_edges)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      GenerateFaces();
      CheckBdrElementOrientation();
   }
   else
   {
      NumOfEdges = 0;
   }

   NumOfFaces = 0;

   SetAttributes();

   SetMeshGen();
}


class GeckoProgress : public Gecko::Progress
{
   double limit;
   mutable StopWatch sw;
public:
   GeckoProgress(double limit) : limit(limit) { sw.Start(); }
   virtual bool quit() const { return limit > 0 && sw.UserTime() > limit; }
};

class GeckoVerboseProgress : public GeckoProgress
{
   using Float = Gecko::Float;
   using Graph = Gecko::Graph;
   using uint = Gecko::uint;
public:
   GeckoVerboseProgress(double limit) : GeckoProgress(limit) {}

   virtual void beginorder(const Graph* graph, Float cost) const
   { mfem::out << "Begin Gecko ordering, cost = " << cost << std::endl; }
   virtual void endorder(const Graph* graph, Float cost) const
   { mfem::out << "End ordering, cost = " << cost << std::endl; }

   virtual void beginiter(const Graph* graph,
                          uint iter, uint maxiter, uint window) const
   {
      mfem::out << "Iteration " << iter << "/" << maxiter << ", window "
                << window << std::flush;
   }
   virtual void enditer(const Graph* graph, Float mincost, Float cost) const
   { mfem::out << ", cost = " << cost << endl; }
};


double Mesh::GetGeckoElementOrdering(Array<int> &ordering,
                                     int iterations, int window,
                                     int period, int seed, bool verbose,
                                     double time_limit)
{
   Gecko::Graph graph;
   Gecko::FunctionalGeometric functional; // edge product cost

   GeckoProgress progress(time_limit);
   GeckoVerboseProgress vprogress(time_limit);

   // insert elements as nodes in the graph
   for (int elemid = 0; elemid < GetNE(); ++elemid)
   {
      graph.insert_node();
   }

   // insert graph edges for element neighbors
   // NOTE: indices in Gecko are 1 based hence the +1 on insertion
   const Table &my_el_to_el = ElementToElementTable();
   for (int elemid = 0; elemid < GetNE(); ++elemid)
   {
      const int *neighid = my_el_to_el.GetRow(elemid);
      for (int i = 0; i < my_el_to_el.RowSize(elemid); ++i)
      {
         graph.insert_arc(elemid + 1,  neighid[i] + 1);
      }
   }

   // get the ordering from Gecko and copy it into the Array<int>
   graph.order(&functional, iterations, window, period, seed,
               verbose ? &vprogress : &progress);

   ordering.SetSize(GetNE());
   Gecko::Node::Index NE = GetNE();
   for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
   {
      ordering[gnodeid - 1] = graph.rank(gnodeid);
   }

   return graph.cost();
}


struct HilbertCmp
{
   int coord;
   bool dir;
   const Array<double> &points;
   double mid;

   HilbertCmp(int coord, bool dir, const Array<double> &points, double mid)
      : coord(coord), dir(dir), points(points), mid(mid) {}

   bool operator()(int i) const
   {
      return (points[3*i + coord] < mid) != dir;
   }
};

static void HilbertSort2D(int coord1, // major coordinate to sort points by
                          bool dir1,  // sort coord1 ascending/descending?
                          bool dir2,  // sort coord2 ascending/descending?
                          const Array<double> &points, int *beg, int *end,
                          double xmin, double ymin, double xmax, double ymax)
{
   if (end - beg <= 1) { return; }

   double xmid = (xmin + xmax)*0.5;
   double ymid = (ymin + ymax)*0.5;

   int coord2 = (coord1 + 1) % 2; // the 'other' coordinate

   // sort (partition) points into four quadrants
   int *p0 = beg, *p4 = end;
   int *p2 = std::partition(p0, p4, HilbertCmp(coord1,  dir1, points, xmid));
   int *p1 = std::partition(p0, p2, HilbertCmp(coord2,  dir2, points, ymid));
   int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));

   if (p1 != p4)
   {
      HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
                    ymin, xmin, ymid, xmid);
   }
   if (p1 != p0 || p2 != p4)
   {
      HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
                    xmin, ymid, xmid, ymax);
   }
   if (p2 != p0 || p3 != p4)
   {
      HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
                    xmid, ymid, xmax, ymax);
   }
   if (p3 != p0)
   {
      HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
                    ymid, xmax, ymin, xmid);
   }
}

static void HilbertSort3D(int coord1, bool dir1, bool dir2, bool dir3,
                          const Array<double> &points, int *beg, int *end,
                          double xmin, double ymin, double zmin,
                          double xmax, double ymax, double zmax)
{
   if (end - beg <= 1) { return; }

   double xmid = (xmin + xmax)*0.5;
   double ymid = (ymin + ymax)*0.5;
   double zmid = (zmin + zmax)*0.5;

   int coord2 = (coord1 + 1) % 3;
   int coord3 = (coord1 + 2) % 3;

   // sort (partition) points into eight octants
   int *p0 = beg, *p8 = end;
   int *p4 = std::partition(p0, p8, HilbertCmp(coord1,  dir1, points, xmid));
   int *p2 = std::partition(p0, p4, HilbertCmp(coord2,  dir2, points, ymid));
   int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
   int *p1 = std::partition(p0, p2, HilbertCmp(coord3,  dir3, points, zmid));
   int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
   int *p5 = std::partition(p4, p6, HilbertCmp(coord3,  dir3, points, zmid));
   int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));

   if (p1 != p8)
   {
      HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
                    zmin, xmin, ymin, zmid, xmid, ymid);
   }
   if (p1 != p0 || p2 != p8)
   {
      HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
                    ymin, zmid, xmin, ymid, zmax, xmid);
   }
   if (p2 != p0 || p3 != p8)
   {
      HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
                    ymid, zmid, xmin, ymax, zmax, xmid);
   }
   if (p3 != p0 || p4 != p8)
   {
      HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
                    xmin, ymax, zmid, xmid, ymid, zmin);
   }
   if (p4 != p0 || p5 != p8)
   {
      HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
                    xmid, ymax, zmid, xmax, ymid, zmin);
   }
   if (p5 != p0 || p6 != p8)
   {
      HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
                    ymax, zmid, xmax, ymid, zmax, xmid);
   }
   if (p6 != p0 || p7 != p8)
   {
      HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
                    ymid, zmid, xmax, ymin, zmax, xmid);
   }
   if (p7 != p0)
   {
      HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
                    zmid, xmax, ymin, zmin, xmid, ymid);
   }
}

void Mesh::GetHilbertElementOrdering(Array<int> &ordering)
{
   MFEM_VERIFY(spaceDim <= 3, "");

   Vector min, max, center;
   GetBoundingBox(min, max);

   Array<int> indices(GetNE());
   Array<double> points(3*GetNE());

   if (spaceDim < 3) { points = 0.0; }

   // calculate element centers
   for (int i = 0; i < GetNE(); i++)
   {
      GetElementCenter(i, center);
      for (int j = 0; j < spaceDim; j++)
      {
         points[3*i + j] = center(j);
      }
      indices[i] = i;
   }

   if (spaceDim == 1)
   {
      indices.Sort([&](int a, int b)
      { return points[3*a] < points[3*b]; });
   }
   else if (spaceDim == 2)
   {
      // recursively partition the points in 2D
      HilbertSort2D(0, false, false,
                    points, indices.begin(), indices.end(),
                    min(0), min(1), max(0), max(1));
   }
   else
   {
      // recursively partition the points in 3D
      HilbertSort3D(0, false, false, false,
                    points, indices.begin(), indices.end(),
                    min(0), min(1), min(2), max(0), max(1), max(2));
   }

   // return ordering in the format required by ReorderElements
   ordering.SetSize(GetNE());
   for (int i = 0; i < GetNE(); i++)
   {
      ordering[indices[i]] = i;
   }
}


void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
{
   if (NURBSext)
   {
      MFEM_WARNING("element reordering of NURBS meshes is not supported.");
      return;
   }
   if (ncmesh)
   {
      MFEM_WARNING("element reordering of non-conforming meshes is not"
                   " supported.");
      return;
   }
   MFEM_VERIFY(ordering.Size() == GetNE(), "invalid reordering array.")

   // Data members that need to be updated:

   // - elements   - reorder of the pointers and the vertex ids if reordering
   //                the vertices
   // - vertices   - if reordering the vertices
   // - boundary   - update the vertex ids, if reordering the vertices
   // - faces      - regenerate
   // - faces_info - regenerate

   // Deleted by DeleteTables():
   // - el_to_edge  - rebuild in 2D and 3D only
   // - el_to_face  - rebuild in 3D only
   // - bel_to_edge - rebuild in 3D only
   // - el_to_el    - no need to rebuild
   // - face_edge   - no need to rebuild
   // - edge_vertex - no need to rebuild
   // - geom_factors - no need to rebuild

   // - be_to_edge  - 2D only
   // - be_to_face  - 3D only

   // - Nodes

   // Save the locations of the Nodes so we can rebuild them later
   Array<Vector*> old_elem_node_vals;
   FiniteElementSpace *nodes_fes = NULL;
   if (Nodes)
   {
      old_elem_node_vals.SetSize(GetNE());
      nodes_fes = Nodes->FESpace();
      Array<int> old_dofs;
      Vector vals;
      for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
      {
         nodes_fes->GetElementVDofs(old_elid, old_dofs);
         Nodes->GetSubVector(old_dofs, vals);
         old_elem_node_vals[old_elid] = new Vector(vals);
      }
   }

   // Get the newly ordered elements
   Array<Element *> new_elements(GetNE());
   for (int old_elid = 0; old_elid < ordering.Size(); ++old_elid)
   {
      int new_elid = ordering[old_elid];
      new_elements[new_elid] = elements[old_elid];
   }
   mfem::Swap(elements, new_elements);
   new_elements.DeleteAll();

   if (reorder_vertices)
   {
      // Get the new vertex ordering permutation vectors and fill the new
      // vertices
      Array<int> vertex_ordering(GetNV());
      vertex_ordering = -1;
      Array<Vertex> new_vertices(GetNV());
      int new_vertex_ind = 0;
      for (int new_elid = 0; new_elid < GetNE(); ++new_elid)
      {
         int *elem_vert = elements[new_elid]->GetVertices();
         int nv = elements[new_elid]->GetNVertices();
         for (int vi = 0; vi < nv; ++vi)
         {
            int old_vertex_ind = elem_vert[vi];
            if (vertex_ordering[old_vertex_ind] == -1)
            {
               vertex_ordering[old_vertex_ind] = new_vertex_ind;
               new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
               new_vertex_ind++;
            }
         }
      }
      mfem::Swap(vertices, new_vertices);
      new_vertices.DeleteAll();

      // Replace the vertex ids in the elements with the reordered vertex
      // numbers
      for (int new_elid = 0; new_elid < GetNE(); ++new_elid)
      {
         int *elem_vert = elements[new_elid]->GetVertices();
         int nv = elements[new_elid]->GetNVertices();
         for (int vi = 0; vi < nv; ++vi)
         {
            elem_vert[vi] = vertex_ordering[elem_vert[vi]];
         }
      }

      // Replace the vertex ids in the boundary with reordered vertex numbers
      for (int belid = 0; belid < GetNBE(); ++belid)
      {
         int *be_vert = boundary[belid]->GetVertices();
         int nv = boundary[belid]->GetNVertices();
         for (int vi = 0; vi < nv; ++vi)
         {
            be_vert[vi] = vertex_ordering[be_vert[vi]];
         }
      }
   }

   // Destroy tables that need to be rebuild
   DeleteTables();

   if (Dim > 1)
   {
      // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
   }
   if (Dim > 2)
   {
      // generate el_to_face, be_to_face
      GetElementToFaceTable();
   }
   // Update faces and faces_info
   GenerateFaces();

   // Build the nodes from the saved locations if they were around before
   if (Nodes)
   {
      // To force FE space update, we need to increase 'sequence':
      sequence++;
      last_operation = Mesh::NONE;
      nodes_fes->Update(false); // want_transform = false
      Nodes->Update(); // just needed to update Nodes->sequence
      Array<int> new_dofs;
      for (int old_elid = 0; old_elid < GetNE(); ++old_elid)
      {
         int new_elid = ordering[old_elid];
         nodes_fes->GetElementVDofs(new_elid, new_dofs);
         Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
         delete old_elem_node_vals[old_elid];
      }
   }
}


void Mesh::MarkForRefinement()
{
   if (meshgen & 1)
   {
      if (Dim == 2)
      {
         MarkTriMeshForRefinement();
      }
      else if (Dim == 3)
      {
         DSTable v_to_v(NumOfVertices);
         GetVertexToVertexTable(v_to_v);
         MarkTetMeshForRefinement(v_to_v);
      }
   }
}

void Mesh::MarkTriMeshForRefinement()
{
   // Mark the longest triangle edge by rotating the indeces so that
   // vertex 0 - vertex 1 is the longest edge in the triangle.
   DenseMatrix pmat;
   for (int i = 0; i < NumOfElements; i++)
   {
      if (elements[i]->GetType() == Element::TRIANGLE)
      {
         GetPointMatrix(i, pmat);
         static_cast<Triangle*>(elements[i])->MarkEdge(pmat);
      }
   }
}

void Mesh::GetEdgeOrdering(DSTable &v_to_v, Array<int> &order)
{
   NumOfEdges = v_to_v.NumberOfEntries();
   order.SetSize(NumOfEdges);
   Array<Pair<double, int> > length_idx(NumOfEdges);

   for (int i = 0; i < NumOfVertices; i++)
   {
      for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
      {
         int j = it.Index();
         length_idx[j].one = GetLength(i, it.Column());
         length_idx[j].two = j;
      }
   }

   // Sort by increasing edge-length.
   length_idx.Sort();

   for (int i = 0; i < NumOfEdges; i++)
   {
      order[length_idx[i].two] = i;
   }
}

void Mesh::MarkTetMeshForRefinement(DSTable &v_to_v)
{
   // Mark the longest tetrahedral edge by rotating the indices so that
   // vertex 0 - vertex 1 is the longest edge in the element.
   Array<int> order;
   GetEdgeOrdering(v_to_v, order);

   for (int i = 0; i < NumOfElements; i++)
   {
      if (elements[i]->GetType() == Element::TETRAHEDRON)
      {
         elements[i]->MarkEdge(v_to_v, order);
      }
   }
   for (int i = 0; i < NumOfBdrElements; i++)
   {
      if (boundary[i]->GetType() == Element::TRIANGLE)
      {
         boundary[i]->MarkEdge(v_to_v, order);
      }
   }
}

void Mesh::PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
{
   if (*old_v_to_v && *old_elem_vert)
   {
      return;
   }

   FiniteElementSpace *fes = Nodes->FESpace();

   if (*old_v_to_v == NULL)
   {
      bool need_v_to_v = false;
      Array<int> dofs;
      for (int i = 0; i < GetNEdges(); i++)
      {
         // Since edge indices may change, we need to permute edge interior dofs
         // any time an edge index changes and there is at least one dof on that
         // edge.
         fes->GetEdgeInteriorDofs(i, dofs);
         if (dofs.Size() > 0)
         {
            need_v_to_v = true;
            break;
         }
      }
      if (need_v_to_v)
      {
         *old_v_to_v = new DSTable(NumOfVertices);
         GetVertexToVertexTable(*(*old_v_to_v));
      }
   }
   if (*old_elem_vert == NULL)
   {
      bool need_elem_vert = false;
      Array<int> dofs;
      for (int i = 0; i < GetNE(); i++)
      {
         // Since element indices do not change, we need to permute element
         // interior dofs only when there are at least 2 interior dofs in an
         // element (assuming the nodal dofs are non-directional).
         fes->GetElementInteriorDofs(i, dofs);
         if (dofs.Size() > 1)
         {
            need_elem_vert = true;
            break;
         }
      }
      if (need_elem_vert)
      {
         *old_elem_vert = new Table;
         (*old_elem_vert)->MakeI(GetNE());
         for (int i = 0; i < GetNE(); i++)
         {
            (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
         }
         (*old_elem_vert)->MakeJ();
         for (int i = 0; i < GetNE(); i++)
         {
            (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
                                             elements[i]->GetNVertices());
         }
         (*old_elem_vert)->ShiftUpI();
      }
   }
}

void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
{
   FiniteElementSpace *fes = Nodes->FESpace();
   const FiniteElementCollection *fec = fes->FEColl();
   Array<int> old_dofs, new_dofs;

   // assuming that all edges have the same number of dofs
   if (NumOfEdges) { fes->GetEdgeInteriorDofs(0, old_dofs); }
   const int num_edge_dofs = old_dofs.Size();

   // Save the original nodes
   const Vector onodes = *Nodes;

   // vertex dofs do not need to be moved
   fes->GetVertexDofs(0, old_dofs);
   int offset = NumOfVertices * old_dofs.Size();

   // edge dofs:
   // edge enumeration may be different but edge orientation is the same
   if (num_edge_dofs > 0)
   {
      DSTable new_v_to_v(NumOfVertices);
      GetVertexToVertexTable(new_v_to_v);

      for (int i = 0; i < NumOfVertices; i++)
      {
         for (DSTable::RowIterator it(new_v_to_v, i); !it; ++it)
         {
            const int old_i = (*old_v_to_v)(i, it.Column());
            const int new_i = it.Index();
            if (new_i == old_i) { continue; }

            old_dofs.SetSize(num_edge_dofs);
            new_dofs.SetSize(num_edge_dofs);
            for (int j = 0; j < num_edge_dofs; j++)
            {
               old_dofs[j] = offset + old_i * num_edge_dofs + j;
               new_dofs[j] = offset + new_i * num_edge_dofs + j;
            }
            fes->DofsToVDofs(old_dofs);
            fes->DofsToVDofs(new_dofs);
            for (int j = 0; j < old_dofs.Size(); j++)
            {
               (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
            }
         }
      }
      offset += NumOfEdges * num_edge_dofs;
   }

   // face dofs:
   // both enumeration and orientation of the faces may be different
   if (fes->GetNFDofs() > 0)
   {
      // generate the old face-vertex table using the unmodified 'faces'
      Table old_face_vertex;
      old_face_vertex.MakeI(NumOfFaces);
      for (int i = 0; i < NumOfFaces; i++)
      {
         old_face_vertex.AddColumnsInRow(i, faces[i]->GetNVertices());
      }
      old_face_vertex.MakeJ();
      for (int i = 0; i < NumOfFaces; i++)
         old_face_vertex.AddConnections(i, faces[i]->GetVertices(),
                                        faces[i]->GetNVertices());
      old_face_vertex.ShiftUpI();

      // update 'el_to_face', 'be_to_face', 'faces', and 'faces_info'
      STable3D *faces_tbl = GetElementToFaceTable(1);
      GenerateFaces();

      // compute the new face dof offsets
      Array<int> new_fdofs(NumOfFaces+1);
      new_fdofs[0] = 0;
      for (int i = 0; i < NumOfFaces; i++) // i = old face index
      {
         const int *old_v = old_face_vertex.GetRow(i);
         int new_i; // new face index
         switch (old_face_vertex.RowSize(i))
         {
            case 3:
               new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
               break;
            case 4:
            default:
               new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
               break;
         }
         fes->GetFaceInteriorDofs(i, old_dofs);
         new_fdofs[new_i+1] = old_dofs.Size();
      }
      new_fdofs.PartialSum();

      // loop over the old face numbers
      for (int i = 0; i < NumOfFaces; i++)
      {
         const int *old_v = old_face_vertex.GetRow(i), *new_v;
         const int *dof_ord;
         int new_i, new_or;
         switch (old_face_vertex.RowSize(i))
         {
            case 3:
               new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
               new_v = faces[new_i]->GetVertices();
               new_or = GetTriOrientation(old_v, new_v);
               dof_ord = fec->DofOrderForOrientation(Geometry::TRIANGLE, new_or);
               break;
            case 4:
            default:
               new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
               new_v = faces[new_i]->GetVertices();
               new_or = GetQuadOrientation(old_v, new_v);
               dof_ord = fec->DofOrderForOrientation(Geometry::SQUARE, new_or);
               break;
         }

         fes->GetFaceInteriorDofs(i, old_dofs);
         new_dofs.SetSize(old_dofs.Size());
         for (int j = 0; j < old_dofs.Size(); j++)
         {
            // we assume the dofs are non-directional, i.e. dof_ord[j] is >= 0
            const int old_j = dof_ord[j];
            new_dofs[old_j] = offset + new_fdofs[new_i] + j;
         }
         fes->DofsToVDofs(old_dofs);
         fes->DofsToVDofs(new_dofs);
         for (int j = 0; j < old_dofs.Size(); j++)
         {
            (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
         }
      }

      offset += fes->GetNFDofs();
      delete faces_tbl;
   }

   // element dofs:
   // element orientation may be different
   if (old_elem_vert) // have elements with 2 or more dofs
   {
      // matters when the 'fec' is
      // (this code is executed only for triangles/tets)
      // - Pk on triangles, k >= 4
      // - Qk on quads,     k >= 3
      // - Pk on tets,      k >= 5
      // - Qk on hexes,     k >= 3
      // - DG spaces
      // - ...

      // loop over all elements
      for (int i = 0; i < GetNE(); i++)
      {
         const int *old_v = old_elem_vert->GetRow(i);
         const int *new_v = elements[i]->GetVertices();
         const int *dof_ord;
         int new_or;
         const Geometry::Type geom = elements[i]->GetGeometryType();
         switch (geom)
         {
            case Geometry::SEGMENT:
               new_or = (old_v[0] == new_v[0]) ? +1 : -1;
               break;
            case Geometry::TRIANGLE:
               new_or = GetTriOrientation(old_v, new_v);
               break;
            case Geometry::SQUARE:
               new_or = GetQuadOrientation(old_v, new_v);
               break;
            case Geometry::TETRAHEDRON:
               new_or = GetTetOrientation(old_v, new_v);
               break;
            default:
               new_or = 0;
               MFEM_ABORT(Geometry::Name[geom] << " elements (" << fec->Name()
                          << " FE collection) are not supported yet!");
               break;
         }
         dof_ord = fec->DofOrderForOrientation(geom, new_or);
         MFEM_VERIFY(dof_ord != NULL,
                     "FE collection '" << fec->Name()
                     << "' does not define reordering for "
                     << Geometry::Name[geom] << " elements!");
         fes->GetElementInteriorDofs(i, old_dofs);
         new_dofs.SetSize(old_dofs.Size());
         for (int j = 0; j < new_dofs.Size(); j++)
         {
            // we assume the dofs are non-directional, i.e. dof_ord[j] is >= 0
            const int old_j = dof_ord[j];
            new_dofs[old_j] = offset + j;
         }
         offset += new_dofs.Size();
         fes->DofsToVDofs(old_dofs);
         fes->DofsToVDofs(new_dofs);
         for (int j = 0; j < old_dofs.Size(); j++)
         {
            (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
         }
      }
   }

   // Update Tables, faces, etc
   if (Dim > 2)
   {
      if (fes->GetNFDofs() == 0)
      {
         // needed for FE spaces that have face dofs, even if
         // the 'Nodes' do not have face dofs.
         GetElementToFaceTable();
         GenerateFaces();
      }
      CheckBdrElementOrientation();
   }
   if (el_to_edge)
   {
      // update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      if (Dim == 2)
      {
         // update 'faces' and 'faces_info'
         GenerateFaces();
         CheckBdrElementOrientation();
      }
   }
   // To force FE space update, we need to increase 'sequence':
   sequence++;
   last_operation = Mesh::NONE;
   fes->Update(false); // want_transform = false
   Nodes->Update(); // just needed to update Nodes->sequence
}

void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
{
   FinalizeCheck();
   CheckElementOrientation(fix_orientation);

   if (NumOfBdrElements == 0)
   {
      GetElementToFaceTable();
      GenerateFaces();
      GenerateBoundaryElements();
   }

   if (refine)
   {
      DSTable v_to_v(NumOfVertices);
      GetVertexToVertexTable(v_to_v);
      MarkTetMeshForRefinement(v_to_v);
   }

   GetElementToFaceTable();
   GenerateFaces();

   CheckBdrElementOrientation();

   if (generate_edges == 1)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
   }
   else
   {
      el_to_edge = NULL;  // Not really necessary -- InitTables was called
      bel_to_edge = NULL;
      NumOfEdges = 0;
   }

   SetAttributes();

   SetMeshGen();
}

void Mesh::FinalizeWedgeMesh(int generate_edges, int refine,
                             bool fix_orientation)
{
   FinalizeCheck();
   CheckElementOrientation(fix_orientation);

   if (NumOfBdrElements == 0)
   {
      GetElementToFaceTable();
      GenerateFaces();
      GenerateBoundaryElements();
   }

   GetElementToFaceTable();
   GenerateFaces();

   CheckBdrElementOrientation();

   if (generate_edges == 1)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
   }
   else
   {
      el_to_edge = NULL;  // Not really necessary -- InitTables was called
      bel_to_edge = NULL;
      NumOfEdges = 0;
   }

   SetAttributes();

   SetMeshGen();
}

void Mesh::FinalizeHexMesh(int generate_edges, int refine, bool fix_orientation)
{
   FinalizeCheck();
   CheckElementOrientation(fix_orientation);

   GetElementToFaceTable();
   GenerateFaces();

   if (NumOfBdrElements == 0)
   {
      GenerateBoundaryElements();
   }

   CheckBdrElementOrientation();

   if (generate_edges)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
   }
   else
   {
      NumOfEdges = 0;
   }

   SetAttributes();

   SetMeshGen();
}

void Mesh::FinalizeMesh(int refine, bool fix_orientation)
{
   FinalizeTopology();

   Finalize(refine, fix_orientation);
}

void Mesh::FinalizeTopology(bool generate_bdr)
{
   // Requirements: the following should be defined:
   //   1) Dim
   //   2) NumOfElements, elements
   //   3) NumOfBdrElements, boundary
   //   4) NumOfVertices
   // Optional:
   //   2) ncmesh may be defined
   //   3) el_to_edge may be allocated (it will be re-computed)

   FinalizeCheck();
   bool generate_edges = true;

   if (spaceDim == 0) { spaceDim = Dim; }
   if (ncmesh) { ncmesh->spaceDim = spaceDim; }

   // if the user defined any hanging nodes (see AddVertexParent),
   // we're initializing a non-conforming mesh
   if (tmp_vertex_parents.Size())
   {
      MFEM_VERIFY(ncmesh == NULL, "");
      ncmesh = new NCMesh(this);

      // we need to recreate the Mesh because NCMesh reorders the vertices
      // (see NCMesh::UpdateVertices())
      InitFromNCMesh(*ncmesh);
      ncmesh->OnMeshUpdated(this);
      GenerateNCFaceInfo();

      SetAttributes();

      tmp_vertex_parents.DeleteAll();
      return;
   }

   // set the mesh type: 'meshgen', ...
   SetMeshGen();

   // generate the faces
   if (Dim > 2)
   {
      GetElementToFaceTable();
      GenerateFaces();
      if (NumOfBdrElements == 0 && generate_bdr)
      {
         GenerateBoundaryElements();
         GetElementToFaceTable(); // update be_to_face
      }
   }
   else
   {
      NumOfFaces = 0;
   }

   // generate edges if requested
   if (Dim > 1 && generate_edges)
   {
      // el_to_edge may already be allocated (P2 VTK meshes)
      if (!el_to_edge) { el_to_edge = new Table; }
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      if (Dim == 2)
      {
         GenerateFaces(); // 'Faces' in 2D refers to the edges
         if (NumOfBdrElements == 0 && generate_bdr)
         {
            GenerateBoundaryElements();
         }
      }
   }
   else
   {
      NumOfEdges = 0;
   }

   if (Dim == 1)
   {
      GenerateFaces();
   }

   if (ncmesh)
   {
      // tell NCMesh the numbering of edges/faces
      ncmesh->OnMeshUpdated(this);

      // update faces_info with NC relations
      GenerateNCFaceInfo();
   }

   // generate the arrays 'attributes' and 'bdr_attributes'
   SetAttributes();
}

void Mesh::Finalize(bool refine, bool fix_orientation)
{
   if (NURBSext || ncmesh)
   {
      MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
      MFEM_ASSERT(CheckBdrElementOrientation() == 0, "");
      return;
   }

   // Requirements:
   //  1) FinalizeTopology() or equivalent was called
   //  2) if (Nodes == NULL), vertices must be defined
   //  3) if (Nodes != NULL), Nodes must be defined

   const bool check_orientation = true; // for regular elements, not boundary
   const bool curved = (Nodes != NULL);
   const bool may_change_topology =
      ( refine && (Dim > 1 && (meshgen & 1)) ) ||
      ( check_orientation && fix_orientation &&
        (Dim == 2 || (Dim == 3 && (meshgen & 1))) );

   DSTable *old_v_to_v = NULL;
   Table *old_elem_vert = NULL;

   if (curved && may_change_topology)
   {
      PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
   }

   if (check_orientation)
   {
      // check and optionally fix element orientation
      CheckElementOrientation(fix_orientation);
   }
   if (refine)
   {
      MarkForRefinement();   // may change topology!
   }

   if (may_change_topology)
   {
      if (curved)
      {
         DoNodeReorder(old_v_to_v, old_elem_vert); // updates the mesh topology
         delete old_elem_vert;
         delete old_v_to_v;
      }
      else
      {
         FinalizeTopology(); // Re-computes some data unnecessarily.
      }

      // TODO: maybe introduce Mesh::NODE_REORDER operation and FESpace::
      // NodeReorderMatrix and do Nodes->Update() instead of DoNodeReorder?
   }

   // check and fix boundary element orientation
   CheckBdrElementOrientation();

#ifdef MFEM_DEBUG
   // For non-orientable surfaces/manifolds, the check below will fail, so we
   // only perform it when Dim == spaceDim.
   if (Dim >= 2 && Dim == spaceDim)
   {
      const int num_faces = GetNumFaces();
      for (int i = 0; i < num_faces; i++)
      {
         MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
                     faces_info[i].Elem2Inf%2 != 0, "Invalid mesh topology."
                     " Interior face with incompatible orientations.");
      }
   }
#endif
}

void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
                  double sx, double sy, double sz, bool sfc_ordering)
{
   int x, y, z;

   int NVert, NElem, NBdrElem;

   NVert = (nx+1) * (ny+1) * (nz+1);
   NElem = nx * ny * nz;
   NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
   if (type == Element::TETRAHEDRON)
   {
      NElem *= 6;
      NBdrElem *= 2;
   }
   else if (type == Element::WEDGE)
   {
      NElem *= 2;
      NBdrElem += 2*nx*ny;
   }

   InitMesh(3, 3, NVert, NElem, NBdrElem);

   double coord[3];
   int ind[8];

   // Sets vertices and the corresponding coordinates
   for (z = 0; z <= nz; z++)
   {
      coord[2] = ((double) z / nz) * sz;
      for (y = 0; y <= ny; y++)
      {
         coord[1] = ((double) y / ny) * sy;
         for (x = 0; x <= nx; x++)
         {
            coord[0] = ((double) x / nx) * sx;
            AddVertex(coord);
         }
      }
   }

#define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))

   // Sets elements and the corresponding indices of vertices
   if (sfc_ordering && type == Element::HEXAHEDRON)
   {
      Array<int> sfc;
      NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
      MFEM_VERIFY(sfc.Size() == 3*nx*ny*nz, "");

      for (int k = 0; k < nx*ny*nz; k++)
      {
         x = sfc[3*k + 0];
         y = sfc[3*k + 1];
         z = sfc[3*k + 2];

         ind[0] = VTX(x  , y  , z  );
         ind[1] = VTX(x+1, y  , z  );
         ind[2] = VTX(x+1, y+1, z  );
         ind[3] = VTX(x  , y+1, z  );
         ind[4] = VTX(x  , y  , z+1);
         ind[5] = VTX(x+1, y  , z+1);
         ind[6] = VTX(x+1, y+1, z+1);
         ind[7] = VTX(x  , y+1, z+1);

         AddHex(ind, 1);
      }
   }
   else
   {
      for (z = 0; z < nz; z++)
      {
         for (y = 0; y < ny; y++)
         {
            for (x = 0; x < nx; x++)
            {
               ind[0] = VTX(x  , y  , z  );
               ind[1] = VTX(x+1, y  , z  );
               ind[2] = VTX(x+1, y+1, z  );
               ind[3] = VTX(x  , y+1, z  );
               ind[4] = VTX(x  , y  , z+1);
               ind[5] = VTX(x+1, y  , z+1);
               ind[6] = VTX(x+1, y+1, z+1);
               ind[7] = VTX(  x, y+1, z+1);
               if (type == Element::TETRAHEDRON)
               {
                  AddHexAsTets(ind, 1);
               }
               else if (type == Element::WEDGE)
               {
                  AddHexAsWedges(ind, 1);
               }
               else
               {
                  AddHex(ind, 1);
               }
            }
         }
      }
   }

   // Sets boundary elements and the corresponding indices of vertices
   // bottom, bdr. attribute 1
   for (y = 0; y < ny; y++)
   {
      for (x = 0; x < nx; x++)
      {
         ind[0] = VTX(x  , y  , 0);
         ind[1] = VTX(x  , y+1, 0);
         ind[2] = VTX(x+1, y+1, 0);
         ind[3] = VTX(x+1, y  , 0);
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 1);
         }
         else if (type == Element::WEDGE)
         {
            AddBdrQuadAsTriangles(ind, 1);
         }
         else
         {
            AddBdrQuad(ind, 1);
         }
      }
   }
   // top, bdr. attribute 6
   for (y = 0; y < ny; y++)
   {
      for (x = 0; x < nx; x++)
      {
         ind[0] = VTX(x  , y  , nz);
         ind[1] = VTX(x+1, y  , nz);
         ind[2] = VTX(x+1, y+1, nz);
         ind[3] = VTX(x  , y+1, nz);
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 6);
         }
         else if (type == Element::WEDGE)
         {
            AddBdrQuadAsTriangles(ind, 1);
         }
         else
         {
            AddBdrQuad(ind, 6);
         }
      }
   }
   // left, bdr. attribute 5
   for (z = 0; z < nz; z++)
   {
      for (y = 0; y < ny; y++)
      {
         ind[0] = VTX(0  , y  , z  );
         ind[1] = VTX(0  , y  , z+1);
         ind[2] = VTX(0  , y+1, z+1);
         ind[3] = VTX(0  , y+1, z  );
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 5);
         }
         else
         {
            AddBdrQuad(ind, 5);
         }
      }
   }
   // right, bdr. attribute 3
   for (z = 0; z < nz; z++)
   {
      for (y = 0; y < ny; y++)
      {
         ind[0] = VTX(nx, y  , z  );
         ind[1] = VTX(nx, y+1, z  );
         ind[2] = VTX(nx, y+1, z+1);
         ind[3] = VTX(nx, y  , z+1);
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 3);
         }
         else
         {
            AddBdrQuad(ind, 3);
         }
      }
   }
   // front, bdr. attribute 2
   for (x = 0; x < nx; x++)
   {
      for (z = 0; z < nz; z++)
      {
         ind[0] = VTX(x  , 0, z  );
         ind[1] = VTX(x+1, 0, z  );
         ind[2] = VTX(x+1, 0, z+1);
         ind[3] = VTX(x  , 0, z+1);
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 2);
         }
         else
         {
            AddBdrQuad(ind, 2);
         }
      }
   }
   // back, bdr. attribute 4
   for (x = 0; x < nx; x++)
   {
      for (z = 0; z < nz; z++)
      {
         ind[0] = VTX(x  , ny, z  );
         ind[1] = VTX(x  , ny, z+1);
         ind[2] = VTX(x+1, ny, z+1);
         ind[3] = VTX(x+1, ny, z  );
         if (type == Element::TETRAHEDRON)
         {
            AddBdrQuadAsTriangles(ind, 4);
         }
         else
         {
            AddBdrQuad(ind, 4);
         }
      }
   }

#undef VTX

#if 0
   ofstream test_stream("debug.mesh");
   Print(test_stream);
   test_stream.close();
#endif

   FinalizeTopology();

   // Finalize(...) can be called after this method, if needed
}

void Mesh::Make2D(int nx, int ny, Element::Type type,
                  double sx, double sy,
                  bool generate_edges, bool sfc_ordering)
{
   int i, j, k;

   SetEmpty();

   Dim = spaceDim = 2;

   // Creates quadrilateral mesh
   if (type == Element::QUADRILATERAL)
   {
      NumOfVertices = (nx+1) * (ny+1);
      NumOfElements = nx * ny;
      NumOfBdrElements = 2 * nx + 2 * ny;

      vertices.SetSize(NumOfVertices);
      elements.SetSize(NumOfElements);
      boundary.SetSize(NumOfBdrElements);

      double cx, cy;
      int ind[4];

      // Sets vertices and the corresponding coordinates
      k = 0;
      for (j = 0; j < ny+1; j++)
      {
         cy = ((double) j / ny) * sy;
         for (i = 0; i < nx+1; i++)
         {
            cx = ((double) i / nx) * sx;
            vertices[k](0) = cx;
            vertices[k](1) = cy;
            k++;
         }
      }

      // Sets elements and the corresponding indices of vertices
      if (sfc_ordering)
      {
         Array<int> sfc;
         NCMesh::GridSfcOrdering2D(nx, ny, sfc);
         MFEM_VERIFY(sfc.Size() == 2*nx*ny, "");

         for (k = 0; k < nx*ny; k++)
         {
            i = sfc[2*k + 0];
            j = sfc[2*k + 1];
            ind[0] = i + j*(nx+1);
            ind[1] = i + 1 +j*(nx+1);
            ind[2] = i + 1 + (j+1)*(nx+1);
            ind[3] = i + (j+1)*(nx+1);
            elements[k] = new Quadrilateral(ind);
         }
      }
      else
      {
         k = 0;
         for (j = 0; j < ny; j++)
         {
            for (i = 0; i < nx; i++)
            {
               ind[0] = i + j*(nx+1);
               ind[1] = i + 1 +j*(nx+1);
               ind[2] = i + 1 + (j+1)*(nx+1);
               ind[3] = i + (j+1)*(nx+1);
               elements[k] = new Quadrilateral(ind);
               k++;
            }
         }
      }

      // Sets boundary elements and the corresponding indices of vertices
      int m = (nx+1)*ny;
      for (i = 0; i < nx; i++)
      {
         boundary[i] = new Segment(i, i+1, 1);
         boundary[nx+i] = new Segment(m+i+1, m+i, 3);
      }
      m = nx+1;
      for (j = 0; j < ny; j++)
      {
         boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
         boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
      }
   }
   // Creates triangular mesh
   else if (type == Element::TRIANGLE)
   {
      NumOfVertices = (nx+1) * (ny+1);
      NumOfElements = 2 * nx * ny;
      NumOfBdrElements = 2 * nx + 2 * ny;

      vertices.SetSize(NumOfVertices);
      elements.SetSize(NumOfElements);
      boundary.SetSize(NumOfBdrElements);

      double cx, cy;
      int ind[3];

      // Sets vertices and the corresponding coordinates
      k = 0;
      for (j = 0; j < ny+1; j++)
      {
         cy = ((double) j / ny) * sy;
         for (i = 0; i < nx+1; i++)
         {
            cx = ((double) i / nx) * sx;
            vertices[k](0) = cx;
            vertices[k](1) = cy;
            k++;
         }
      }

      // Sets the elements and the corresponding indices of vertices
      k = 0;
      for (j = 0; j < ny; j++)
      {
         for (i = 0; i < nx; i++)
         {
            ind[0] = i + j*(nx+1);
            ind[1] = i + 1 + (j+1)*(nx+1);
            ind[2] = i + (j+1)*(nx+1);
            elements[k] = new Triangle(ind);
            k++;
            ind[1] = i + 1 + j*(nx+1);
            ind[2] = i + 1 + (j+1)*(nx+1);
            elements[k] = new Triangle(ind);
            k++;
         }
      }

      // Sets boundary elements and the corresponding indices of vertices
      int m = (nx+1)*ny;
      for (i = 0; i < nx; i++)
      {
         boundary[i] = new Segment(i, i+1, 1);
         boundary[nx+i] = new Segment(m+i+1, m+i, 3);
      }
      m = nx+1;
      for (j = 0; j < ny; j++)
      {
         boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
         boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
      }

      // MarkTriMeshForRefinement(); // done in Finalize(...)
   }
   else
   {
      MFEM_ABORT("Unsupported element type.");
   }

   SetMeshGen();
   CheckElementOrientation();

   if (generate_edges == 1)
   {
      el_to_edge = new Table;
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      GenerateFaces();
      CheckBdrElementOrientation();
   }
   else
   {
      NumOfEdges = 0;
   }

   NumOfFaces = 0;

   attributes.Append(1);
   bdr_attributes.Append(1); bdr_attributes.Append(2);
   bdr_attributes.Append(3); bdr_attributes.Append(4);

   // Finalize(...) can be called after this method, if needed
}

void Mesh::Make1D(int n, double sx)
{
   int j, ind[1];

   SetEmpty();

   Dim = 1;
   spaceDim = 1;

   NumOfVertices = n + 1;
   NumOfElements = n;
   NumOfBdrElements = 2;
   vertices.SetSize(NumOfVertices);
   elements.SetSize(NumOfElements);
   boundary.SetSize(NumOfBdrElements);

   // Sets vertices and the corresponding coordinates
   for (j = 0; j < n+1; j++)
   {
      vertices[j](0) = ((double) j / n) * sx;
   }

   // Sets elements and the corresponding indices of vertices
   for (j = 0; j < n; j++)
   {
      elements[j] = new Segment(j, j+1, 1);
   }

   // Sets the boundary elements
   ind[0] = 0;
   boundary[0] = new Point(ind, 1);
   ind[0] = n;
   boundary[1] = new Point(ind, 2);

   NumOfEdges = 0;
   NumOfFaces = 0;

   SetMeshGen();
   GenerateFaces();

   attributes.Append(1);
   bdr_attributes.Append(1); bdr_attributes.Append(2);
}

Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
{
   Dim = mesh.Dim;
   spaceDim = mesh.spaceDim;

   NumOfVertices = mesh.NumOfVertices;
   NumOfElements = mesh.NumOfElements;
   NumOfBdrElements = mesh.NumOfBdrElements;
   NumOfEdges = mesh.NumOfEdges;
   NumOfFaces = mesh.NumOfFaces;
   nbInteriorFaces = mesh.nbInteriorFaces;
   nbBoundaryFaces = mesh.nbBoundaryFaces;

   meshgen = mesh.meshgen;
   mesh_geoms = mesh.mesh_geoms;

   // Create the new Mesh instance without a record of its refinement history
   sequence = 0;
   last_operation = Mesh::NONE;

   // Duplicate the elements
   elements.SetSize(NumOfElements);
   for (int i = 0; i < NumOfElements; i++)
   {
      elements[i] = mesh.elements[i]->Duplicate(this);
   }

   // Copy the vertices
   mesh.vertices.Copy(vertices);

   // Duplicate the boundary
   boundary.SetSize(NumOfBdrElements);
   for (int i = 0; i < NumOfBdrElements; i++)
   {
      boundary[i] = mesh.boundary[i]->Duplicate(this);
   }

   // Copy the element-to-face Table, el_to_face
   el_to_face = (mesh.el_to_face) ? new Table(*mesh.el_to_face) : NULL;

   // Copy the boundary-to-face Array, be_to_face.
   mesh.be_to_face.Copy(be_to_face);

   // Copy the element-to-edge Table, el_to_edge
   el_to_edge = (mesh.el_to_edge) ? new Table(*mesh.el_to_edge) : NULL;

   // Copy the boundary-to-edge Table, bel_to_edge (3D)
   bel_to_edge = (mesh.bel_to_edge) ? new Table(*mesh.bel_to_edge) : NULL;

   // Copy the boundary-to-edge Array, be_to_edge (2D)
   mesh.be_to_edge.Copy(be_to_edge);

   // Duplicate the faces and faces_info.
   faces.SetSize(mesh.faces.Size());
   for (int i = 0; i < faces.Size(); i++)
   {
      Element *face = mesh.faces[i]; // in 1D the faces are NULL
      faces[i] = (face) ? face->Duplicate(this) : NULL;
   }
   mesh.faces_info.Copy(faces_info);
   mesh.nc_faces_info.Copy(nc_faces_info);

   // Do NOT copy the element-to-element Table, el_to_el
   el_to_el = NULL;

   // Do NOT copy the face-to-edge Table, face_edge
   face_edge = NULL;

   // Copy the edge-to-vertex Table, edge_vertex
   edge_vertex = (mesh.edge_vertex) ? new Table(*mesh.edge_vertex) : NULL;

   // Copy the attributes and bdr_attributes
   mesh.attributes.Copy(attributes);
   mesh.bdr_attributes.Copy(bdr_attributes);

   // Deep copy the NURBSExtension.
#ifdef MFEM_USE_MPI
   ParNURBSExtension *pNURBSext =
      dynamic_cast<ParNURBSExtension *>(mesh.NURBSext);
   if (pNURBSext)
   {
      NURBSext = new ParNURBSExtension(*pNURBSext);
   }
   else
#endif
   {
      NURBSext = mesh.NURBSext ? new NURBSExtension(*mesh.NURBSext) : NULL;
   }

   // Deep copy the NCMesh.
#ifdef MFEM_USE_MPI
   if (dynamic_cast<const ParMesh*>(&mesh))
   {
      ncmesh = NULL; // skip; will be done in ParMesh copy ctor
   }
   else
#endif
   {
      ncmesh = mesh.ncmesh ? new NCMesh(*mesh.ncmesh) : NULL;
   }

   // Duplicate the Nodes, including the FiniteElementCollection and the
   // FiniteElementSpace
   if (mesh.Nodes && copy_nodes)
   {
      FiniteElementSpace *fes = mesh.Nodes->FESpace();
      const FiniteElementCollection *fec = fes->FEColl();
      FiniteElementCollection *fec_copy =
         FiniteElementCollection::New(fec->Name());
      FiniteElementSpace *fes_copy =
         new FiniteElementSpace(*fes, this, fec_copy);
      Nodes = new GridFunction(fes_copy);
      Nodes->MakeOwner(fec_copy);
      *Nodes = *mesh.Nodes;
      own_nodes = 1;
   }
   else
   {
      Nodes = mesh.Nodes;
      own_nodes = 0;
   }
}

Mesh::Mesh(Mesh &&mesh) : Mesh()
{
   Swap(mesh, true);
}

Mesh& Mesh::operator=(Mesh &&mesh)
{
   Swap(mesh, true);
   return *this;
}

Mesh Mesh::LoadFromFile(const char *filename, int generate_edges, int refine,
                        bool fix_orientation)
{
   Mesh mesh;
   named_ifgzstream imesh(filename);
   if (!imesh) { MFEM_ABORT("Mesh file not found: " << filename << '\n'); }
   else { mesh.Load(imesh, generate_edges, refine, fix_orientation); }
   return mesh;
}

Mesh Mesh::MakeCartesian1D(int n, double sx)
{
   Mesh mesh;
   mesh.Make1D(n, sx);
   // mesh.Finalize(); not needed in this case
   return mesh;
}

Mesh Mesh::MakeCartesian2D(
   int nx, int ny, Element::Type type, bool generate_edges,
   double sx, double sy, bool sfc_ordering)
{
   Mesh mesh;
   mesh.Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
   mesh.Finalize(true); // refine = true
   return mesh;
}

Mesh Mesh::MakeCartesian3D(
   int nx, int ny, int nz, Element::Type type,
   double sx, double sy, double sz, bool sfc_ordering)
{
   Mesh mesh;
   mesh.Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
   mesh.Finalize(true); // refine = true
   return mesh;
}

Mesh Mesh::MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
{
   Mesh mesh;
   Array<int> ref_factors(orig_mesh.GetNE());
   ref_factors = ref_factor;
   mesh.MakeRefined_(orig_mesh, ref_factors, ref_type);
   return mesh;
}

Mesh Mesh::MakeRefined(Mesh &orig_mesh, const Array<int> &ref_factors,
                       int ref_type)
{
   Mesh mesh;
   mesh.MakeRefined_(orig_mesh, ref_factors, ref_type);
   return mesh;
}

Mesh::Mesh(const char *filename, int generate_edges, int refine,
           bool fix_orientation)
{
   // Initialization as in the default constructor
   SetEmpty();

   named_ifgzstream imesh(filename);
   if (!imesh)
   {
      // Abort with an error message.
      MFEM_ABORT("Mesh file not found: " << filename << '\n');
   }
   else
   {
      Load(imesh, generate_edges, refine, fix_orientation);
   }
}

Mesh::Mesh(std::istream &input, int generate_edges, int refine,
           bool fix_orientation)
{
   SetEmpty();
   Load(input, generate_edges, refine, fix_orientation);
}

void Mesh::ChangeVertexDataOwnership(double *vertex_data, int len_vertex_data,
                                     bool zerocopy)
{
   // A dimension of 3 is now required since we use mfem::Vertex objects as PODs
   // and these object have a hardcoded double[3] entry
   MFEM_VERIFY(len_vertex_data >= NumOfVertices * 3,
               "Not enough vertices in external array : "
               "len_vertex_data = "<< len_vertex_data << ", "
               "NumOfVertices * 3 = " << NumOfVertices * 3);
   // Allow multiple calls to this method with the same vertex_data
   if (vertex_data == (double *)(vertices.GetData()))
   {
      MFEM_ASSERT(!vertices.OwnsData(), "invalid ownership");
      return;
   }
   if (!zerocopy)
   {
      memcpy(vertex_data, vertices.GetData(),
             NumOfVertices * 3 * sizeof(double));
   }
   // Vertex is POD double[3]
   vertices.MakeRef(reinterpret_cast<Vertex*>(vertex_data), NumOfVertices);
}

Mesh::Mesh(double *vertices_, int num_vertices,
           int *element_indices, Geometry::Type element_type,
           int *element_attributes, int num_elements,
           int *boundary_indices, Geometry::Type boundary_type,
           int *boundary_attributes, int num_boundary_elements,
           int dimension, int space_dimension)
{
   if (space_dimension == -1)
   {
      space_dimension = dimension;
   }

   InitMesh(dimension, space_dimension, /*num_vertices*/ 0, num_elements,
            num_boundary_elements);

   int element_index_stride = Geometry::NumVerts[element_type];
   int boundary_index_stride = num_boundary_elements > 0 ?
                               Geometry::NumVerts[boundary_type] : 0;

   // assuming Vertex is POD
   vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
   NumOfVertices = num_vertices;

   for (int i = 0; i < num_elements; i++)
   {
      elements[i] = NewElement(element_type);
      elements[i]->SetVertices(element_indices + i * element_index_stride);
      elements[i]->SetAttribute(element_attributes[i]);
   }
   NumOfElements = num_elements;

   for (int i = 0; i < num_boundary_elements; i++)
   {
      boundary[i] = NewElement(boundary_type);
      boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
      boundary[i]->SetAttribute(boundary_attributes[i]);
   }
   NumOfBdrElements = num_boundary_elements;

   FinalizeTopology();
}

Element *Mesh::NewElement(int geom)
{
   switch (geom)
   {
      case Geometry::POINT:     return (new Point);
      case Geometry::SEGMENT:   return (new Segment);
      case Geometry::TRIANGLE:  return (new Triangle);
      case Geometry::SQUARE:    return (new Quadrilateral);
      case Geometry::TETRAHEDRON:
#ifdef MFEM_USE_MEMALLOC
         return TetMemory.Alloc();
#else
         return (new Tetrahedron);
#endif
      case Geometry::CUBE:      return (new Hexahedron);
      case Geometry::PRISM:     return (new Wedge);
      default:
         MFEM_ABORT("invalid Geometry::Type, geom = " << geom);
   }

   return NULL;
}

Element *Mesh::ReadElementWithoutAttr(std::istream &input)
{
   int geom, nv, *v;
   Element *el;

   input >> geom;
   el = NewElement(geom);
   MFEM_VERIFY(el, "Unsupported element type: " << geom);
   nv = el->GetNVertices();
   v  = el->GetVertices();
   for (int i = 0; i < nv; i++)
   {
      input >> v[i];
   }

   return el;
}

void Mesh::PrintElementWithoutAttr(const Element *el, std::ostream &out)
{
   out << el->GetGeometryType();
   const int nv = el->GetNVertices();
   const int *v = el->GetVertices();
   for (int j = 0; j < nv; j++)
   {
      out << ' ' << v[j];
   }
   out << '\n';
}

Element *Mesh::ReadElement(std::istream &input)
{
   int attr;
   Element *el;

   input >> attr;
   el = ReadElementWithoutAttr(input);
   el->SetAttribute(attr);

   return el;
}

void Mesh::PrintElement(const Element *el, std::ostream &out)
{
   out << el->GetAttribute() << ' ';
   PrintElementWithoutAttr(el, out);
}

void Mesh::SetMeshGen()
{
   meshgen = mesh_geoms = 0;
   for (int i = 0; i < NumOfElements; i++)
   {
      const Element::Type type = GetElement(i)->GetType();
      switch (type)
      {
         case Element::TETRAHEDRON:
            mesh_geoms |= (1 << Geometry::TETRAHEDRON);
         case Element::TRIANGLE:
            mesh_geoms |= (1 << Geometry::TRIANGLE);
         case Element::SEGMENT:
            mesh_geoms |= (1 << Geometry::SEGMENT);
         case Element::POINT:
            mesh_geoms |= (1 << Geometry::POINT);
            meshgen |= 1;
            break;

         case Element::HEXAHEDRON:
            mesh_geoms |= (1 << Geometry::CUBE);
         case Element::QUADRILATERAL:
            mesh_geoms |= (1 << Geometry::SQUARE);
            mesh_geoms |= (1 << Geometry::SEGMENT);
            mesh_geoms |= (1 << Geometry::POINT);
            meshgen |= 2;
            break;

         case Element::WEDGE:
            mesh_geoms |= (1 << Geometry::PRISM);
            mesh_geoms |= (1 << Geometry::SQUARE);
            mesh_geoms |= (1 << Geometry::TRIANGLE);
            mesh_geoms |= (1 << Geometry::SEGMENT);
            mesh_geoms |= (1 << Geometry::POINT);
            meshgen |= 4;
            break;

         default:
            MFEM_ABORT("invalid element type: " << type);
            break;
      }
   }
}

void Mesh::Loader(std::istream &input, int generate_edges,
                  std::string parse_tag)
{
   int curved = 0, read_gf = 1;
   bool finalize_topo = true;

   if (!input)
   {
      MFEM_ABORT("Input stream is not open");
   }

   Clear();

   string mesh_type;
   input >> ws;
   getline(input, mesh_type);
   filter_dos(mesh_type);

   // MFEM's conforming mesh formats
   int mfem_version = 0;
   if (mesh_type == "MFEM mesh v1.0") { mfem_version = 10; } // serial
   else if (mesh_type == "MFEM mesh v1.2") { mfem_version = 12; } // parallel

   // MFEM nonconforming mesh format
   // (NOTE: previous v1.1 is now under this branch for backward compatibility)
   int mfem_nc_version = 0;
   if (mesh_type == "MFEM NC mesh v1.0") { mfem_nc_version = 10; }
   else if (mesh_type == "MFEM mesh v1.1") { mfem_nc_version = 1 /*legacy*/; }

   if (mfem_version)
   {
      // Formats mfem_v12 and newer have a tag indicating the end of the mesh
      // section in the stream. A user provided parse tag can also be provided
      // via the arguments. For example, if this is called from parallel mesh
      // object, it can indicate to read until parallel mesh section begins.
      if (mfem_version == 12 && parse_tag.empty())
      {
         parse_tag = "mfem_mesh_end";
      }
      ReadMFEMMesh(input, mfem_version, curved);
   }
   else if (mfem_nc_version)
   {
      MFEM_ASSERT(ncmesh == NULL, "internal error");
      int is_nc = 1;

#ifdef MFEM_USE_MPI
      ParMesh *pmesh = dynamic_cast<ParMesh*>(this);
      if (pmesh)
      {
         MFEM_VERIFY(mfem_nc_version >= 10,
                     "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
                     "used to load a parallel nonconforming mesh, sorry.");

         ncmesh = new ParNCMesh(pmesh->GetComm(),
                                input, mfem_nc_version, curved, is_nc);
      }
      else
#endif
      {
         ncmesh = new NCMesh(input, mfem_nc_version, curved, is_nc);
      }
      InitFromNCMesh(*ncmesh);

      if (!is_nc)
      {
         // special case for backward compatibility with MFEM <=4.2:
         // if the "vertex_parents" section is missing in the v1.1 format,
         // the mesh is treated as conforming
         delete ncmesh;
         ncmesh = NULL;
      }
   }
   else if (mesh_type == "linemesh") // 1D mesh
   {
      ReadLineMesh(input);
   }
   else if (mesh_type == "areamesh2" || mesh_type == "curved_areamesh2")
   {
      if (mesh_type == "curved_areamesh2")
      {
         curved = 1;
      }
      ReadNetgen2DMesh(input, curved);
   }
   else if (mesh_type == "NETGEN" || mesh_type == "NETGEN_Neutral_Format")
   {
      ReadNetgen3DMesh(input);
   }
   else if (mesh_type == "TrueGrid")
   {
      ReadTrueGridMesh(input);
   }
   else if (mesh_type.rfind("# vtk DataFile Version") == 0)
   {
      int major_vtk_version = mesh_type[mesh_type.length()-3] - '0';
      // int minor_vtk_version = mesh_type[mesh_type.length()-1] - '0';
      MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
                  "Unsupported VTK format");
      ReadVTKMesh(input, curved, read_gf, finalize_topo);
   }
   else if (mesh_type.rfind("<VTKFile ") == 0 || mesh_type.rfind("<?xml") == 0)
   {
      ReadXML_VTKMesh(input, curved, read_gf, finalize_topo, mesh_type);
   }
   else if (mesh_type == "MFEM NURBS mesh v1.0")
   {
      ReadNURBSMesh(input, curved, read_gf);
   }
   else if (mesh_type == "MFEM INLINE mesh v1.0")
   {
      ReadInlineMesh(input, generate_edges);
      return; // done with inline mesh construction
   }
   else if (mesh_type == "$MeshFormat") // Gmsh
   {
      ReadGmshMesh(input, curved, read_gf);
   }
   else if
   ((mesh_type.size() > 2 &&
     mesh_type[0] == 'C' && mesh_type[1] == 'D' && mesh_type[2] == 'F') ||
    (mesh_type.size() > 3 &&
     mesh_type[1] == 'H' && mesh_type[2] == 'D' && mesh_type[3] == 'F'))
   {
      named_ifgzstream *mesh_input = dynamic_cast<named_ifgzstream *>(&input);
      if (mesh_input)
      {
#ifdef MFEM_USE_NETCDF
         ReadCubit(mesh_input->filename.c_str(), curved, read_gf);
#else
         MFEM_ABORT("NetCDF support requires configuration with"
                    " MFEM_USE_NETCDF=YES");
         return;
#endif
      }
      else
      {
         MFEM_ABORT("Can not determine Cubit mesh filename!"
                    " Use mfem::named_ifgzstream for input.");
         return;
      }
   }
   else
   {
      MFEM_ABORT("Unknown input mesh format: " << mesh_type);
      return;
   }

   // at this point the following should be defined:
   //  1) Dim
   //  2) NumOfElements, elements
   //  3) NumOfBdrElements, boundary
   //  4) NumOfVertices, with allocated space in vertices
   //  5) curved
   //  5a) if curved == 0, vertices must be defined
   //  5b) if curved != 0 and read_gf != 0,
   //         'input' must point to a GridFunction
   //  5c) if curved != 0 and read_gf == 0,
   //         vertices and Nodes must be defined
   // optional:
   //  1) el_to_edge may be allocated (as in the case of P2 VTK meshes)
   //  2) ncmesh may be allocated

   // FinalizeTopology() will:
   // - assume that generate_edges is true
   // - assume that refine is false
   // - does not check the orientation of regular and boundary elements
   if (finalize_topo)
   {
      // don't generate any boundary elements, especially in parallel
      bool generate_bdr = false;

      FinalizeTopology(generate_bdr);
   }

   if (curved && read_gf)
   {
      Nodes = new GridFunction(this, input);

      own_nodes = 1;
      spaceDim = Nodes->VectorDim();
      if (ncmesh) { ncmesh->spaceDim = spaceDim; }

      // Set vertex coordinates from the 'Nodes'
      SetVerticesFromNodes(Nodes);
   }

   // If a parse tag was supplied, keep reading the stream until the tag is
   // encountered.
   if (mfem_version == 12)
   {
      string line;
      do
      {
         skip_comment_lines(input, '#');
         MFEM_VERIFY(input.good(), "Required mesh-end tag not found");
         getline(input, line);
         filter_dos(line);
         // mfem v1.2 may not have parse_tag in it, e.g. if trying to read a
         // serial mfem v1.2 mesh as parallel with "mfem_serial_mesh_end" as
         // parse_tag. That's why, regardless of parse_tag, we stop reading if
         // we find "mfem_mesh_end" which is required by mfem v1.2 format.
         if (line == "mfem_mesh_end") { break; }
      }
      while (line != parse_tag);
   }
   else if (mfem_nc_version >= 10)
   {
      string ident;
      skip_comment_lines(input, '#');
      input >> ident;
      MFEM_VERIFY(ident == "mfem_mesh_end",
                  "invalid mesh: end of file tag not found");
   }

   // Finalize(...) should be called after this, if needed.
}

Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
{
   int      i, j, ie, ib, iv, *v, nv;
   Element *el;
   Mesh    *m;

   SetEmpty();

   Dim = mesh_array[0]->Dimension();
   spaceDim = mesh_array[0]->SpaceDimension();

   if (mesh_array[0]->NURBSext)
   {
      // assuming the pieces form a partition of a NURBS mesh
      NURBSext = new NURBSExtension(mesh_array, num_pieces);

      NumOfVertices = NURBSext->GetNV();
      NumOfElements = NURBSext->GetNE();

      NURBSext->GetElementTopo(elements);

      // NumOfBdrElements = NURBSext->GetNBE();
      // NURBSext->GetBdrElementTopo(boundary);

      Array<int> lvert_vert, lelem_elem;

      // Here, for visualization purposes, we copy the boundary elements from
      // the individual pieces which include the interior boundaries.  This
      // creates 'boundary' array that is different from the one generated by
      // the NURBSExtension which, in particular, makes the boundary-dof table
      // invalid. This, in turn, causes GetBdrElementTransformation to not
      // function properly.
      NumOfBdrElements = 0;
      for (i = 0; i < num_pieces; i++)
      {
         NumOfBdrElements += mesh_array[i]->GetNBE();
      }
      boundary.SetSize(NumOfBdrElements);
      vertices.SetSize(NumOfVertices);
      ib = 0;
      for (i = 0; i < num_pieces; i++)
      {
         m = mesh_array[i];
         m->NURBSext->GetVertexLocalToGlobal(lvert_vert);
         m->NURBSext->GetElementLocalToGlobal(lelem_elem);
         // copy the element attributes
         for (j = 0; j < m->GetNE(); j++)
         {
            elements[lelem_elem[j]]->SetAttribute(m->GetAttribute(j));
         }
         // copy the boundary
         for (j = 0; j < m->GetNBE(); j++)
         {
            el = m->GetBdrElement(j)->Duplicate(this);
            v  = el->GetVertices();
            nv = el->GetNVertices();
            for (int k = 0; k < nv; k++)
            {
               v[k] = lvert_vert[v[k]];
            }
            boundary[ib++] = el;
         }
         // copy the vertices
         for (j = 0; j < m->GetNV(); j++)
         {
            vertices[lvert_vert[j]].SetCoords(m->SpaceDimension(),
                                              m->GetVertex(j));
         }
      }
   }
   else // not a NURBS mesh
   {
      NumOfElements    = 0;
      NumOfBdrElements = 0;
      NumOfVertices    = 0;
      for (i = 0; i < num_pieces; i++)
      {
         m = mesh_array[i];
         NumOfElements    += m->GetNE();
         NumOfBdrElements += m->GetNBE();
         NumOfVertices    += m->GetNV();
      }
      elements.SetSize(NumOfElements);
      boundary.SetSize(NumOfBdrElements);
      vertices.SetSize(NumOfVertices);
      ie = ib = iv = 0;
      for (i = 0; i < num_pieces; i++)
      {
         m = mesh_array[i];
         // copy the elements
         for (j = 0; j < m->GetNE(); j++)
         {
            el = m->GetElement(j)->Duplicate(this);
            v  = el->GetVertices();
            nv = el->GetNVertices();
            for (int k = 0; k < nv; k++)
            {
               v[k] += iv;
            }
            elements[ie++] = el;
         }
         // copy the boundary elements
         for (j = 0; j < m->GetNBE(); j++)
         {
            el = m->GetBdrElement(j)->Duplicate(this);
            v  = el->GetVertices();
            nv = el->GetNVertices();
            for (int k = 0; k < nv; k++)
            {
               v[k] += iv;
            }
            boundary[ib++] = el;
         }
         // copy the vertices
         for (j = 0; j < m->GetNV(); j++)
         {
            vertices[iv++].SetCoords(m->SpaceDimension(), m->GetVertex(j));
         }
      }
   }

   FinalizeTopology();

   // copy the nodes (curvilinear meshes)
   GridFunction *g = mesh_array[0]->GetNodes();
   if (g)
   {
      Array<GridFunction *> gf_array(num_pieces);
      for (i = 0; i < num_pieces; i++)
      {
         gf_array[i] = mesh_array[i]->GetNodes();
      }
      Nodes = new GridFunction(this, gf_array, num_pieces);
      own_nodes = 1;
   }

#ifdef MFEM_DEBUG
   CheckElementOrientation(false);
   CheckBdrElementOrientation(false);
#endif
}

Mesh::Mesh(Mesh *orig_mesh, int ref_factor, int ref_type)
{
   Array<int> ref_factors(orig_mesh->GetNE());
   ref_factors = ref_factor;
   MakeRefined_(*orig_mesh, ref_factors, ref_type);
}

void Mesh::MakeRefined_(Mesh &orig_mesh, const Array<int> ref_factors,
                        int ref_type)
{
   SetEmpty();
   Dim = orig_mesh.Dimension();
   spaceDim = orig_mesh.SpaceDimension();

   int orig_ne = orig_mesh.GetNE();
   MFEM_VERIFY(ref_factors.Size() == orig_ne && orig_ne > 0,
               "Number of refinement factors must equal number of elements")
   MFEM_VERIFY(ref_factors.Min() >= 1, "Refinement factor must be >= 1");
   const int q_type = BasisType::GetQuadrature1D(ref_type);
   MFEM_VERIFY(Quadrature1D::CheckClosed(q_type) != Quadrature1D::Invalid,
               "Invalid refinement type. Must use closed basis type.");

   int min_ref = ref_factors.Min();
   int max_ref = ref_factors.Max();

   bool var_order = (min_ref != max_ref);

   // variable order space can only be constructed over an NC mesh
   if (var_order) { orig_mesh.EnsureNCMesh(true); }

   // Construct a scalar H1 FE space of order ref_factor and use its dofs as
   // the indices of the new, refined vertices.
   H1_FECollection rfec(min_ref, Dim, ref_type);
   FiniteElementSpace rfes(&orig_mesh, &rfec);

   if (var_order)
   {
      rfes.SetRelaxedHpConformity(false);
      for (int i = 0; i < orig_ne; i++)
      {
         rfes.SetElementOrder(i, ref_factors[i]);
      }
      rfes.Update(false);
   }

   // Set the number of vertices, set the actual coordinates later
   NumOfVertices = rfes.GetNDofs();
   vertices.SetSize(NumOfVertices);

   Array<int> rdofs;
   DenseMatrix phys_pts;

   GeometryRefiner refiner;
   refiner.SetType(q_type);

   // Add refined elements and set vertex coordinates
   for (int el = 0; el < orig_ne; el++)
   {
      Geometry::Type geom = orig_mesh.GetElementGeometry(el);
      int attrib = orig_mesh.GetAttribute(el);
      int nvert = Geometry::NumVerts[geom];
      RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[el]);

      rfes.GetElementDofs(el, rdofs);
      MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
      const FiniteElement *rfe = rfes.GetFE(el);
      orig_mesh.GetElementTransformation(el)->Transform(rfe->GetNodes(),
                                                        phys_pts);
      const int *c2h_map = rfec.GetDofMap(geom, ref_factors[el]);
      for (int i = 0; i < phys_pts.Width(); i++)
      {
         vertices[rdofs[i]].SetCoords(spaceDim, phys_pts.GetColumn(i));
      }
      for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
      {
         Element *elem = NewElement(geom);
         elem->SetAttribute(attrib);
         int *v = elem->GetVertices();
         for (int k = 0; k < nvert; k++)
         {
            int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
            v[k] = rdofs[c2h_map[cid]];
         }
         AddElement(elem);
      }
   }

   // Add refined boundary elements
   for (int el = 0; el < orig_mesh.GetNBE(); el++)
   {
      int i, info;
      orig_mesh.GetBdrElementAdjacentElement(el, i, info);
      Geometry::Type geom = orig_mesh.GetBdrElementBaseGeometry(el);
      int attrib = orig_mesh.GetBdrAttribute(el);
      int nvert = Geometry::NumVerts[geom];
      RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[i]);

      rfes.GetBdrElementDofs(el, rdofs);
      MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
      const int *c2h_map = rfec.GetDofMap(geom, ref_factors[i]);
      for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
      {
         Element *elem = NewElement(geom);
         elem->SetAttribute(attrib);
         int *v = elem->GetVertices();
         for (int k = 0; k < nvert; k++)
         {
            int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
            v[k] = rdofs[c2h_map[cid]];
         }
         AddBdrElement(elem);
      }
   }
   FinalizeTopology(false);
   sequence = orig_mesh.GetSequence() + 1;
   last_operation = Mesh::REFINE;

   // Set up the nodes of the new mesh (if the original mesh has nodes). The new
   // mesh is always straight-sided (i.e. degree 1 finite element space), but
   // the nodes are required for e.g. periodic meshes.
   if (orig_mesh.GetNodes())
   {
      bool discont = orig_mesh.GetNodalFESpace()->IsDGSpace();
      Ordering::Type dof_ordering = orig_mesh.GetNodalFESpace()->GetOrdering();
      Mesh::SetCurvature(1, discont, spaceDim, dof_ordering);
      FiniteElementSpace *nodal_fes = Nodes->FESpace();
      const FiniteElementCollection *nodal_fec = nodal_fes->FEColl();
      H1_FECollection vertex_fec(1, Dim);
      Array<int> dofs;
      int el_counter = 0;
      for (int iel = 0; iel < orig_ne; iel++)
      {
         Geometry::Type geom = orig_mesh.GetElementBaseGeometry(iel);
         int nvert = Geometry::NumVerts[geom];
         RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[iel]);
         rfes.GetElementDofs(iel, rdofs);
         const FiniteElement *rfe = rfes.GetFE(iel);
         orig_mesh.GetElementTransformation(iel)->Transform(rfe->GetNodes(),
                                                            phys_pts);
         const int *node_map = NULL;
         const H1_FECollection *h1_fec =
            dynamic_cast<const H1_FECollection *>(nodal_fec);
         if (h1_fec != NULL) { node_map = h1_fec->GetDofMap(geom); }
         const int *vertex_map = vertex_fec.GetDofMap(geom);
         const int *c2h_map = rfec.GetDofMap(geom, ref_factors[iel]);
         for (int jel = 0; jel < RG.RefGeoms.Size()/nvert; jel++)
         {
            nodal_fes->GetElementVDofs(el_counter++, dofs);
            for (int iv_lex=0; iv_lex<nvert; ++iv_lex)
            {
               // convert from lexicographic to vertex index
               int iv = vertex_map[iv_lex];
               // index of vertex of current element in phys_pts matrix
               int pt_idx = c2h_map[RG.RefGeoms[iv+nvert*jel]];
               // index of current vertex into DOF array
               int node_idx = node_map ? node_map[iv_lex] : iv_lex;
               for (int d=0; d<spaceDim; ++d)
               {
                  (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
               }
            }
         }
      }
   }

   // Setup the data for the coarse-fine refinement transformations
   CoarseFineTr.embeddings.SetSize(GetNE());
   // First, compute total number of point matrices that we need per geometry
   // and the offsets into that array
   using GeomRef = std::pair<Geometry::Type, int>;
   std::map<GeomRef, int> point_matrices_offsets;
   int n_point_matrices[Geometry::NumGeom] = {}; // initialize to zero
   for (int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
   {
      Geometry::Type geom = orig_mesh.GetElementBaseGeometry(el_coarse);
      // Have we seen this pair of (goemetry, refinement level) before?
      GeomRef id(geom, ref_factors[el_coarse]);
      if (point_matrices_offsets.find(id) == point_matrices_offsets.end())
      {
         RefinedGeometry &RG = *refiner.Refine(geom, ref_factors[el_coarse]);
         int nvert = Geometry::NumVerts[geom];
         int nref_el = RG.RefGeoms.Size()/nvert;
         // If not, then store the offset and add to the size required
         point_matrices_offsets[id] = n_point_matrices[geom];
         n_point_matrices[geom] += nref_el;
      }
   }

   // Set up the sizes
   for (int geom = 0; geom < Geometry::NumGeom; ++geom)
   {
      int nmatrices = n_point_matrices[geom];
      int nvert = Geometry::NumVerts[geom];
      CoarseFineTr.point_matrices[geom].SetSize(Dim, nvert, nmatrices);
   }

   // Compute the point matrices and embeddings
   int el_fine = 0;
   for (int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
   {
      Geometry::Type geom = orig_mesh.GetElementBaseGeometry(el_coarse);
      int ref = ref_factors[el_coarse];
      int offset = point_matrices_offsets[GeomRef(geom, ref)];
      int nvert = Geometry::NumVerts[geom];
      RefinedGeometry &RG = *refiner.Refine(geom, ref);
      for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
      {
         DenseMatrix &Pj = CoarseFineTr.point_matrices[geom](offset + j);
         for (int k = 0; k < nvert; k++)
         {
            int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
            const IntegrationPoint &ip = RG.RefPts[cid];
            ip.Get(Pj.GetColumn(k), Dim);
         }

         Embedding &emb = CoarseFineTr.embeddings[el_fine];
         emb.parent = el_coarse;
         emb.matrix = offset + j;
         ++el_fine;
      }
   }

   MFEM_ASSERT(CheckElementOrientation(false) == 0, "");

   // The check below is disabled because is fails for parallel meshes with
   // interior "boundary" element that, when such "boundary" element is between
   // two elements on different processors.
   // MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
}

Mesh Mesh::MakeSimplicial(const Mesh &orig_mesh)
{
   Mesh mesh;
   mesh.MakeSimplicial_(orig_mesh, NULL);
   return mesh;
}

void Mesh::MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
{
   MFEM_VERIFY(const_cast<Mesh&>(orig_mesh).CheckElementOrientation(false) == 0,
               "Mesh::MakeSimplicial requires a properly oriented input mesh");
   MFEM_VERIFY(orig_mesh.Conforming(),
               "Mesh::MakeSimplicial does not support non-conforming meshes.")

   int dim = orig_mesh.Dimension();
   int sdim = orig_mesh.SpaceDimension();

   if (dim == 1)
   {
      Mesh copy(orig_mesh);
      Swap(copy, true);
      return;
   }

   int nv = orig_mesh.GetNV();
   int ne = orig_mesh.GetNE();
   int nbe = orig_mesh.GetNBE();

   static int num_subdivisions[Geometry::NUM_GEOMETRIES];
   num_subdivisions[Geometry::POINT] = 1;
   num_subdivisions[Geometry::SEGMENT] = 1;
   num_subdivisions[Geometry::TRIANGLE] = 1;
   num_subdivisions[Geometry::TETRAHEDRON] = 1;
   num_subdivisions[Geometry::SQUARE] = 2;
   num_subdivisions[Geometry::PRISM] = 3;
   num_subdivisions[Geometry::CUBE] = 6;
   // NOTE: some hexes may be subdivided into only 5 tets, so this is an
   // estimate only. The actual number of created tets may be less, so the
   // elements array will need to be shrunk after mesh creation.
   int new_ne = 0, new_nbe = 0;
   for (int i=0; i<ne; ++i)
   {
      new_ne += num_subdivisions[orig_mesh.GetElementBaseGeometry(i)];
   }
   for (int i=0; i<nbe; ++i)
   {
      new_nbe += num_subdivisions[orig_mesh.GetBdrElementBaseGeometry(i)];
   }

   InitMesh(dim, sdim, nv, new_ne, new_nbe);

   // Vertices of the new mesh are same as the original mesh
   NumOfVertices = nv;
   for (int i=0; i<nv; ++i)
   {
      vertices[i].SetCoords(dim, orig_mesh.vertices[i]());
   }

   // We need a global vertex numbering to identify which diagonals to split
   // (quad faces are split using the diagonal originating from the smallest
   // global vertex number). Use the supplied global numbering, if it is
   // non-NULL, otherwise use the local numbering.
   Array<int> vglobal_id;
   if (vglobal == NULL)
   {
      vglobal_id.SetSize(nv);
      for (int i=0; i<nv; ++i) { vglobal_id[i] = i; }
      vglobal = vglobal_id.GetData();
   }

   constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
   constexpr int quad_ntris = 2, prism_ntets = 3;
   static const int quad_trimap[2][nv_tri*quad_ntris] =
   {
      {
         0, 0,
         1, 2,
         2, 3
      },{
         0, 1,
         1, 2,
         3, 3
      }
   };
   static const int prism_rot[nv_prism*nv_prism] =
   {
      0, 1, 2, 3, 4, 5,
      1, 2, 0, 4, 5, 3,
      2, 0, 1, 5, 3, 4,
      3, 5, 4, 0, 2, 1,
      4, 3, 5, 1, 0, 2,
      5, 4, 3, 2, 1, 0
   };
   static const int prism_f[nv_quad] = {1, 2, 5, 4};
   static const int prism_tetmaps[2][nv_prism*prism_ntets] =
   {
      {
         0, 0, 0,
         1, 1, 4,
         2, 5, 5,
         5, 4, 3
      },{
         0, 0, 0,
         1, 4, 4,
         2, 2, 5,
         4, 5, 3
      }
   };
   static const int hex_rot[nv_hex*nv_hex] =
   {
      0, 1, 2, 3, 4, 5, 6, 7,
      1, 0, 4, 5, 2, 3, 7, 6,
      2, 1, 5, 6, 3, 0, 4, 7,
      3, 0, 1, 2, 7, 4, 5, 6,
      4, 0, 3, 7, 5, 1, 2, 6,
      5, 1, 0, 4, 6, 2, 3, 7,
      6, 2, 1, 5, 7, 3, 0, 4,
      7, 3, 2, 6, 4, 0, 1, 5
   };
   static const int hex_f0[nv_quad] = {1, 2, 6, 5};
   static const int hex_f1[nv_quad] = {2, 3, 7, 6};
   static const int hex_f2[nv_quad] = {4, 5, 6, 7};
   static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
   static const int hex_tetmap0[nv_tet*5] =
   {
      0, 0, 0, 0, 2,
      1, 2, 2, 5, 7,
      2, 7, 3, 7, 5,
      5, 5, 7, 4, 6
   };
   static const int hex_tetmap1[nv_tet*6] =
   {
      0, 0, 1, 0, 0, 1,
      5, 1, 6, 7, 7, 7,
      7, 7, 7, 2, 1, 6,
      4, 5, 5, 3, 2, 2
   };
   static const int hex_tetmap2[nv_tet*6] =
   {
      0, 0, 0, 0, 0, 0,
      4, 3, 7, 1, 3, 6,
      5, 7, 4, 2, 6, 5,
      6, 6, 6, 5, 2, 2
   };
   static const int hex_tetmap3[nv_tet*6] =
   {
      0, 0, 0, 0, 1, 1,
      2, 3, 7, 5, 5, 6,
      3, 7, 4, 6, 6, 2,
      6, 6, 6, 4, 0, 0
   };
   static const int *hex_tetmaps[4] =
   {
      hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
   };

   auto find_min = [](const int*a, int n) { return std::min_element(a,a+n)-a; };

   for (int i=0; i<ne; ++i)
   {
      const int *v = orig_mesh.elements[i]->GetVertices();
      const int attrib = orig_mesh.GetAttribute(i);
      const Geometry::Type orig_geom = orig_mesh.GetElementBaseGeometry(i);

      if (num_subdivisions[orig_geom] == 1)
      {
         // (num_subdivisions[orig_geom] == 1) implies that the element does
         // not need to be further split (it is either a segment, triangle,
         // or tetrahedron), and so it is left unchanged.
         Element *e = NewElement(orig_geom);
         e->SetAttribute(attrib);
         e->SetVertices(v);
         AddElement(e);
      }
      else if (orig_geom == Geometry::SQUARE)
      {
         for (int itri=0; itri<quad_ntris; ++itri)
         {
            Element *e = NewElement(Geometry::TRIANGLE);
            e->SetAttribute(attrib);
            int *v2 = e->GetVertices();
            for (int iv=0; iv<nv_tri; ++iv)
            {
               v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
            }
            AddElement(e);
         }
      }
      else if (orig_geom == Geometry::PRISM)
      {
         int vg[nv_prism];
         for (int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
         // Rotate the vertices of the prism so that the smallest vertex index
         // is in the first place
         int irot = find_min(vg, nv_prism);
         for (int iv=0; iv<nv_prism; ++iv)
         {
            int jv = prism_rot[iv + irot*nv_prism];
            vg[iv] = v[jv];
         }
         // Two cases according to which diagonal splits third quad face
         int q[nv_quad];
         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
         int j = find_min(q, nv_quad);
         const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
         for (int itet=0; itet<prism_ntets; ++itet)
         {
            Element *e = NewElement(Geometry::TETRAHEDRON);
            e->SetAttribute(attrib);
            int *v2 = e->GetVertices();
            for (int iv=0; iv<nv_tet; ++iv)
            {
               v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
            }
            AddElement(e);
         }
      }
      else if (orig_geom == Geometry::CUBE)
      {
         int vg[nv_hex];
         for (int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }

         // Rotate the vertices of the hex so that the smallest vertex index is
         // in the first place
         int irot = find_min(vg, nv_hex);
         for (int iv=0; iv<nv_hex; ++iv)
         {
            int jv = hex_rot[iv + irot*nv_hex];
            vg[iv] = v[jv];
         }

         int q[nv_quad];
         // Bitmask is three binary digits, each digit is 1 if the diagonal of
         // the corresponding face goes through the 7th vertex, and 0 if not.
         int bitmask = 0;
         int j;
         // First quad face
         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
         j = find_min(q, nv_quad);
         if (j == 0 || j == 2) { bitmask += 4; }
         // Second quad face
         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
         j = find_min(q, nv_quad);
         if (j == 1 || j == 3) { bitmask += 2; }
         // Third quad face
         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
         j = find_min(q, nv_quad);
         if (j == 0 || j == 2) { bitmask += 1; }

         // Apply rotations
         int nrot = num_rot[bitmask];
         for (int irot=0; irot<nrot; ++irot)
         {
            int vtemp;
            vtemp = vg[1];
            vg[1] = vg[4];
            vg[4] = vg[3];
            vg[3] = vtemp;
            vtemp = vg[5];
            vg[5] = vg[7];
            vg[7] = vg[2];
            vg[2] = vtemp;
         }

         // Sum up nonzero bits in bitmask
         int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
         int ntets = (ndiags == 0) ? 5 : 6;
         const int *tetmap = hex_tetmaps[ndiags];
         for (int itet=0; itet<ntets; ++itet)
         {
            Element *e = NewElement(Geometry::TETRAHEDRON);
            e->SetAttribute(attrib);
            int *v2 = e->GetVertices();
            for (int iv=0; iv<nv_tet; ++iv)
            {
               v2[iv] = vg[tetmap[itet + iv*ntets]];
            }
            AddElement(e);
         }
      }
   }
   // In 3D, shrink the element array because some hexes have only 5 tets
   if (dim == 3) { elements.SetSize(NumOfElements); }

   for (int i=0; i<nbe; ++i)
   {
      const int *v = orig_mesh.boundary[i]->GetVertices();
      const int attrib = orig_mesh.GetBdrAttribute(i);
      const Geometry::Type orig_geom = orig_mesh.GetBdrElementBaseGeometry(i);
      if (num_subdivisions[orig_geom] == 1)
      {
         Element *be = NewElement(orig_geom);
         be->SetAttribute(attrib);
         be->SetVertices(v);
         AddBdrElement(be);
      }
      else if (orig_geom == Geometry::SQUARE)
      {
         int vg[nv_quad];
         for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
         // Split quad according the smallest (global) vertex
         int iv_min = find_min(vg, nv_quad);
         int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
         for (int itri=0; itri<quad_ntris; ++itri)
         {
            Element *be = NewElement(Geometry::TRIANGLE);
            be->SetAttribute(attrib);
            int *v2 = be->GetVertices();
            for (int iv=0; iv<nv_tri; ++iv)
            {
               v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
            }
            AddBdrElement(be);
         }
      }
      else
      {
         MFEM_ABORT("Unreachable");
      }
   }

   FinalizeTopology(false);
   sequence = orig_mesh.GetSequence();
   last_operation = orig_mesh.last_operation;

   MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
   MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
}

Mesh Mesh::MakePeriodic(const Mesh &orig_mesh, const std::vector<int> &v2v)
{
   Mesh periodic_mesh(orig_mesh, true); // Make a copy of the original mesh
   const FiniteElementSpace *nodal_fes = orig_mesh.GetNodalFESpace();
   int nodal_order = nodal_fes ? nodal_fes->GetMaxElementOrder() : 1;
   periodic_mesh.SetCurvature(nodal_order, true);

   // renumber element vertices
   for (int i = 0; i < periodic_mesh.GetNE(); i++)
   {
      Element *el = periodic_mesh.GetElement(i);
      int *v = el->GetVertices();
      int nv = el->GetNVertices();
      for (int j = 0; j < nv; j++)
      {
         v[j] = v2v[v[j]];
      }
   }
   // renumber boundary element vertices
   for (int i = 0; i < periodic_mesh.GetNBE(); i++)
   {
      Element *el = periodic_mesh.GetBdrElement(i);
      int *v = el->GetVertices();
      int nv = el->GetNVertices();
      for (int j = 0; j < nv; j++)
      {
         v[j] = v2v[v[j]];
      }
   }

   periodic_mesh.RemoveUnusedVertices();
   return periodic_mesh;
}

std::vector<int> Mesh::CreatePeriodicVertexMapping(
   const std::vector<Vector> &translations, double tol) const
{
   int sdim = SpaceDimension();

   Vector coord(sdim), at(sdim), dx(sdim);
   Vector xMax(sdim), xMin(sdim), xDiff(sdim);
   xMax = xMin = xDiff = 0.0;

   // Get a list of all vertices on the boundary
   set<int> bdr_v;
   for (int be = 0; be < GetNBE(); be++)
   {
      Array<int> dofs;
      GetBdrElementVertices(be,dofs);

      for (int i = 0; i < dofs.Size(); i++)
      {
         bdr_v.insert(dofs[i]);

         coord = GetVertex(dofs[i]);
         for (int j = 0; j < sdim; j++)
         {
            xMax[j] = max(xMax[j], coord[j]);
            xMin[j] = min(xMin[j], coord[j]);
         }
      }
   }
   add(xMax, -1.0, xMin, xDiff);
   double dia = xDiff.Norml2(); // compute mesh diameter

   // We now identify coincident vertices. Several originally distinct vertices
   // may become coincident under the periodic mapping. One of these vertices
   // will be identified as the "primary" vertex, and all other coincident
   // vertices will be considered as "replicas".

   // replica2primary[v] is the index of the primary vertex of replica `v`
   map<int, int> replica2primary;
   // primary2replicas[v] is a set of indices of replicas of primary vertex `v`
   map<int, set<int>> primary2replicas;

   // We begin with the assumption that all vertices are primary, and that there
   // are no replicas.
   for (int v : bdr_v) { primary2replicas[v]; }

   // Make `r` and all of `r`'s replicas be replicas of `p`. Delete `r` from the
   // list of primary vertices.
   auto make_replica = [&replica2primary, &primary2replicas](int r, int p)
   {
      if (r == p) { return; }
      primary2replicas[p].insert(r);
      replica2primary[r] = p;
      for (int s : primary2replicas[r])
      {
         primary2replicas[p].insert(s);
         replica2primary[s] = p;
      }
      primary2replicas.erase(r);
   };

   for (unsigned int i = 0; i < translations.size(); i++)
   {
      for (int vi : bdr_v)
      {
         coord = GetVertex(vi);
         add(coord, translations[i], at);

         for (int vj : bdr_v)
         {
            coord = GetVertex(vj);
            add(at, -1.0, coord, dx);
            if (dx.Norml2() > dia*tol) { continue; }

            // The two vertices vi and vj are coincident.

            // Are vertices `vi` and `vj` already primary?
            bool pi = primary2replicas.find(vi) != primary2replicas.end();
            bool pj = primary2replicas.find(vj) != primary2replicas.end();

            if (pi && pj)
            {
               // Both vertices are currently primary
               // Demote `vj` to be a replica of `vi`
               make_replica(vj, vi);
            }
            else if (pi && !pj)
            {
               // `vi` is primary and `vj` is a replica
               int owner_of_vj = replica2primary[vj];
               // Make `vi` and its replicas be replicas of `vj`'s owner
               make_replica(vi, owner_of_vj);
            }
            else if (!pi && pj)
            {
               // `vi` is currently a replica and `vj` is currently primary
               // Make `vj` and its replicas be replicas of `vi`'s owner
               int owner_of_vi = replica2primary[vi];
               make_replica(vj, owner_of_vi);
            }
            else
            {
               // Both vertices are currently replicas
               // Make `vj`'s owner and all of its owner's replicas be replicas
               // of `vi`'s owner
               int owner_of_vi = replica2primary[vi];
               int owner_of_vj = replica2primary[vj];
               make_replica(owner_of_vj, owner_of_vi);
            }
            break;
         }
      }
   }

   std::vector<int> v2v(GetNV());
   for (size_t i = 0; i < v2v.size(); i++)
   {
      v2v[i] = i;
   }
   for (auto &&r2p : replica2primary)
   {
      v2v[r2p.first] = r2p.second;
   }
   return v2v;
}

void Mesh::KnotInsert(Array<KnotVector *> &kv)
{
   if (NURBSext == NULL)
   {
      mfem_error("Mesh::KnotInsert : Not a NURBS mesh!");
   }

   if (kv.Size() != NURBSext->GetNKV())
   {
      mfem_error("Mesh::KnotInsert : KnotVector array size mismatch!");
   }

   NURBSext->ConvertToPatches(*Nodes);

   NURBSext->KnotInsert(kv);

   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
   sequence++;

   UpdateNURBS();
}

void Mesh::KnotInsert(Array<Vector *> &kv)
{
   if (NURBSext == NULL)
   {
      mfem_error("Mesh::KnotInsert : Not a NURBS mesh!");
   }

   if (kv.Size() != NURBSext->GetNKV())
   {
      mfem_error("Mesh::KnotInsert : KnotVector array size mismatch!");
   }

   NURBSext->ConvertToPatches(*Nodes);

   NURBSext->KnotInsert(kv);

   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
   sequence++;

   UpdateNURBS();
}

void Mesh::NURBSUniformRefinement()
{
   // do not check for NURBSext since this method is protected
   NURBSext->ConvertToPatches(*Nodes);

   NURBSext->UniformRefinement();

   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
   sequence++;

   UpdateNURBS();
}

void Mesh::DegreeElevate(int rel_degree, int degree)
{
   if (NURBSext == NULL)
   {
      mfem_error("Mesh::DegreeElevate : Not a NURBS mesh!");
   }

   NURBSext->ConvertToPatches(*Nodes);

   NURBSext->DegreeElevate(rel_degree, degree);

   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
   sequence++;

   UpdateNURBS();
}

void Mesh::UpdateNURBS()
{
   ResetLazyData();

   NURBSext->SetKnotsFromPatches();

   Dim = NURBSext->Dimension();
   spaceDim = Dim;

   if (NumOfElements != NURBSext->GetNE())
   {
      for (int i = 0; i < elements.Size(); i++)
      {
         FreeElement(elements[i]);
      }
      NumOfElements = NURBSext->GetNE();
      NURBSext->GetElementTopo(elements);
   }

   if (NumOfBdrElements != NURBSext->GetNBE())
   {
      for (int i = 0; i < boundary.Size(); i++)
      {
         FreeElement(boundary[i]);
      }
      NumOfBdrElements = NURBSext->GetNBE();
      NURBSext->GetBdrElementTopo(boundary);
   }

   Nodes->FESpace()->Update();
   Nodes->Update();
   NURBSext->SetCoordsFromPatches(*Nodes);

   if (NumOfVertices != NURBSext->GetNV())
   {
      NumOfVertices = NURBSext->GetNV();
      vertices.SetSize(NumOfVertices);
      int vd = Nodes->VectorDim();
      for (int i = 0; i < vd; i++)
      {
         Vector vert_val;
         Nodes->GetNodalValues(vert_val, i+1);
         for (int j = 0; j < NumOfVertices; j++)
         {
            vertices[j](i) = vert_val(j);
         }
      }
   }

   if (el_to_edge)
   {
      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      if (Dim == 2)
      {
         GenerateFaces();
      }
   }

   if (el_to_face)
   {
      GetElementToFaceTable();
      GenerateFaces();
   }
}

void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_knot)
{
   SetEmpty();

   // Read MFEM NURBS mesh v1.0 format
   string ident;

   skip_comment_lines(input, '#');

   input >> ident; // 'dimension'
   input >> Dim;
   spaceDim = Dim;

   skip_comment_lines(input, '#');

   input >> ident; // 'elements'
   input >> NumOfElements;
   elements.SetSize(NumOfElements);
   for (int j = 0; j < NumOfElements; j++)
   {
      elements[j] = ReadElement(input);
   }

   skip_comment_lines(input, '#');

   input >> ident; // 'boundary'
   input >> NumOfBdrElements;
   boundary.SetSize(NumOfBdrElements);
   for (int j = 0; j < NumOfBdrElements; j++)
   {
      boundary[j] = ReadElement(input);
   }

   skip_comment_lines(input, '#');

   input >> ident; // 'edges'
   input >> NumOfEdges;
   edge_vertex = new Table(NumOfEdges, 2);
   edge_to_knot.SetSize(NumOfEdges);
   for (int j = 0; j < NumOfEdges; j++)
   {
      int *v = edge_vertex->GetRow(j);
      input >> edge_to_knot[j] >> v[0] >> v[1];
      if (v[0] > v[1])
      {
         edge_to_knot[j] = -1 - edge_to_knot[j];
      }
   }

   skip_comment_lines(input, '#');

   input >> ident; // 'vertices'
   input >> NumOfVertices;
   vertices.SetSize(0);

   FinalizeTopology();
   CheckBdrElementOrientation(); // check and fix boundary element orientation
}

void XYZ_VectorFunction(const Vector &p, Vector &v)
{
   if (p.Size() >= v.Size())
   {
      for (int d = 0; d < v.Size(); d++)
      {
         v(d) = p(d);
      }
   }
   else
   {
      int d;
      for (d = 0; d < p.Size(); d++)
      {
         v(d) = p(d);
      }
      for ( ; d < v.Size(); d++)
      {
         v(d) = 0.0;
      }
   }
}

void Mesh::GetNodes(GridFunction &nodes) const
{
   if (Nodes == NULL || Nodes->FESpace() != nodes.FESpace())
   {
      const int newSpaceDim = nodes.FESpace()->GetVDim();
      VectorFunctionCoefficient xyz(newSpaceDim, XYZ_VectorFunction);
      nodes.ProjectCoefficient(xyz);
   }
   else
   {
      nodes = *Nodes;
   }
}

void Mesh::SetNodalFESpace(FiniteElementSpace *nfes)
{
   GridFunction *nodes = new GridFunction(nfes);
   SetNodalGridFunction(nodes, true);
}

void Mesh::EnsureNodes()
{
   if (Nodes)
   {
      const FiniteElementCollection *fec = GetNodalFESpace()->FEColl();
      if (dynamic_cast<const H1_FECollection*>(fec)
          || dynamic_cast<const L2_FECollection*>(fec))
      {
         return;
      }
      else // Mesh using a legacy FE_Collection
      {
         const int order = GetNodalFESpace()->GetElementOrder(0);
         SetCurvature(order, false, -1, Ordering::byVDIM);
      }
   }
   else // First order H1 mesh
   {
      SetCurvature(1, false, -1, Ordering::byVDIM);
   }
}

void Mesh::SetNodalGridFunction(GridFunction *nodes, bool make_owner)
{
   GetNodes(*nodes);
   NewNodes(*nodes, make_owner);
}

const FiniteElementSpace *Mesh::GetNodalFESpace() const
{
   return ((Nodes) ? Nodes->FESpace() : NULL);
}

void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
{
   space_dim = (space_dim == -1) ? spaceDim : space_dim;
   FiniteElementCollection* nfec;
   if (discont)
   {
      const int type = 1; // Gauss-Lobatto points
      nfec = new L2_FECollection(order, Dim, type);
   }
   else
   {
      nfec = new H1_FECollection(order, Dim);
   }
   FiniteElementSpace* nfes = new FiniteElementSpace(this, nfec, space_dim,
                                                     ordering);
   SetNodalFESpace(nfes);
   Nodes->MakeOwner(nfec);
}

void Mesh::SetVerticesFromNodes(const GridFunction *nodes)
{
   MFEM_ASSERT(nodes != NULL, "");
   for (int i = 0; i < spaceDim; i++)
   {
      Vector vert_val;
      nodes->GetNodalValues(vert_val, i+1);
      for (int j = 0; j < NumOfVertices; j++)
      {
         vertices[j](i) = vert_val(j);
      }
   }
}

int Mesh::GetNumFaces() const
{
   switch (Dim)
   {
      case 1: return GetNV();
      case 2: return GetNEdges();
      case 3: return GetNFaces();
   }
   return 0;
}

static int CountFacesByType(const Mesh &mesh, const FaceType type)
{
   int e1, e2;
   int inf1, inf2;
   int nf = 0;
   for (int f = 0; f < mesh.GetNumFaces(); ++f)
   {
      mesh.GetFaceElements(f, &e1, &e2);
      mesh.GetFaceInfos(f, &inf1, &inf2);
      if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
          (type==FaceType::Boundary && e2<0 && inf2<0) ) { nf++; }
   }
   return nf;
}

int Mesh::GetNFbyType(FaceType type) const
{
   const bool isInt = type==FaceType::Interior;
   int &nf = isInt ? nbInteriorFaces : nbBoundaryFaces;
   if (nf<0) { nf = CountFacesByType(*this, type); }
   return nf;
}

#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
#endif

int Mesh::CheckElementOrientation(bool fix_it)
{
   int i, j, k, wo = 0, fo = 0;
   double *v[4];

   if (Dim == 2 && spaceDim == 2)
   {
      DenseMatrix J(2, 2);

      for (i = 0; i < NumOfElements; i++)
      {
         int *vi = elements[i]->GetVertices();
         if (Nodes == NULL)
         {
            for (j = 0; j < 3; j++)
            {
               v[j] = vertices[vi[j]]();
            }
            for (j = 0; j < 2; j++)
               for (k = 0; k < 2; k++)
               {
                  J(j, k) = v[j+1][k] - v[0][k];
               }
         }
         else
         {
            // only check the Jacobian at the center of the element
            GetElementJacobian(i, J);
         }
         if (J.Det() < 0.0)
         {
            if (fix_it)
            {
               switch (GetElementType(i))
               {
                  case Element::TRIANGLE:
                     mfem::Swap(vi[0], vi[1]);
                     break;
                  case Element::QUADRILATERAL:
                     mfem::Swap(vi[1], vi[3]);
                     break;
                  default:
                     MFEM_ABORT("Invalid 2D element type \""
                                << GetElementType(i) << "\"");
                     break;
               }
               fo++;
            }
            wo++;
         }
      }
   }

   if (Dim == 3)
   {
      DenseMatrix J(3, 3);

      for (i = 0; i < NumOfElements; i++)
      {
         int *vi = elements[i]->GetVertices();
         switch (GetElementType(i))
         {
            case Element::TETRAHEDRON:
               if (Nodes == NULL)
               {
                  for (j = 0; j < 4; j++)
                  {
                     v[j] = vertices[vi[j]]();
                  }
                  for (j = 0; j < 3; j++)
                     for (k = 0; k < 3; k++)
                     {
                        J(j, k) = v[j+1][k] - v[0][k];
                     }
               }
               else
               {
                  // only check the Jacobian at the center of the element
                  GetElementJacobian(i, J);
               }
               if (J.Det() < 0.0)
               {
                  wo++;
                  if (fix_it)
                  {
                     mfem::Swap(vi[0], vi[1]);
                     fo++;
                  }
               }
               break;

            case Element::WEDGE:
               // only check the Jacobian at the center of the element
               GetElementJacobian(i, J);
               if (J.Det() < 0.0)
               {
                  wo++;
                  if (fix_it)
                  {
                     // how?
                  }
               }
               break;

            case Element::HEXAHEDRON:
               // only check the Jacobian at the center of the element
               GetElementJacobian(i, J);
               if (J.Det() < 0.0)
               {
                  wo++;
                  if (fix_it)
                  {
                     // how?
                  }
               }
               break;

            default:
               MFEM_ABORT("Invalid 3D element type \""
                          << GetElementType(i) << "\"");
               break;
         }
      }
   }
#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
   if (wo > 0)
   {
      mfem::out << "Elements with wrong orientation: " << wo << " / "
                << NumOfElements << " (" << fixed_or_not[(wo == fo) ? 0 : 1]
                << ")" << endl;
   }
#endif
   return wo;
}

int Mesh::GetTriOrientation(const int *base, const int *test)
{
   // Static method.
   // This function computes the index 'j' of the permutation that transforms
   // test into base: test[tri_orientation[j][i]]=base[i].
   // tri_orientation = Geometry::Constants<Geometry::TRIANGLE>::Orient
   int orient;

   if (test[0] == base[0])
      if (test[1] == base[1])
      {
         orient = 0;   //  (0, 1, 2)
      }
      else
      {
         orient = 5;   //  (0, 2, 1)
      }
   else if (test[0] == base[1])
      if (test[1] == base[0])
      {
         orient = 1;   //  (1, 0, 2)
      }
      else
      {
         orient = 2;   //  (1, 2, 0)
      }
   else // test[0] == base[2]
      if (test[1] == base[0])
      {
         orient = 4;   //  (2, 0, 1)
      }
      else
      {
         orient = 3;   //  (2, 1, 0)
      }

#ifdef MFEM_DEBUG
   const int *aor = tri_t::Orient[orient];
   for (int j = 0; j < 3; j++)
      if (test[aor[j]] != base[j])
      {
         mfem_error("Mesh::GetTriOrientation(...)");
      }
#endif

   return orient;
}

int Mesh::GetQuadOrientation(const int *base, const int *test)
{
   int i;

   for (i = 0; i < 4; i++)
      if (test[i] == base[0])
      {
         break;
      }

#ifdef MFEM_DEBUG
   int orient;
   if (test[(i+1)%4] == base[1])
   {
      orient = 2*i;
   }
   else
   {
      orient = 2*i+1;
   }
   const int *aor = quad_t::Orient[orient];
   for (int j = 0; j < 4; j++)
      if (test[aor[j]] != base[j])
      {
         mfem::err << "Mesh::GetQuadOrientation(...)" << endl;
         mfem::err << " base = [";
         for (int k = 0; k < 4; k++)
         {
            mfem::err << " " << base[k];
         }
         mfem::err << " ]\n test = [";
         for (int k = 0; k < 4; k++)
         {
            mfem::err << " " << test[k];
         }
         mfem::err << " ]" << endl;
         mfem_error();
      }
#endif

   if (test[(i+1)%4] == base[1])
   {
      return 2*i;
   }

   return 2*i+1;
}

int Mesh::GetTetOrientation(const int *base, const int *test)
{
   // Static method.
   // This function computes the index 'j' of the permutation that transforms
   // test into base: test[tet_orientation[j][i]]=base[i].
   // tet_orientation = Geometry::Constants<Geometry::TETRAHEDRON>::Orient
   int orient;

   if (test[0] == base[0])
      if (test[1] == base[1])
         if (test[2] == base[2])
         {
            orient = 0;   //  (0, 1, 2, 3)
         }
         else
         {
            orient = 1;   //  (0, 1, 3, 2)
         }
      else if (test[2] == base[1])
         if (test[3] == base[2])
         {
            orient = 2;   //  (0, 2, 3, 1)
         }
         else
         {
            orient = 3;   //  (0, 2, 1, 3)
         }
      else // test[3] == base[1]
         if (test[1] == base[2])
         {
            orient = 4;   //  (0, 3, 1, 2)
         }
         else
         {
            orient = 5;   //  (0, 3, 2, 1)
         }
   else if (test[1] == base[0])
      if (test[2] == base[1])
         if (test[0] == base[2])
         {
            orient = 6;   //  (1, 2, 0, 3)
         }
         else
         {
            orient = 7;   //  (1, 2, 3, 0)
         }
      else if (test[3] == base[1])
         if (test[2] == base[2])
         {
            orient = 8;   //  (1, 3, 2, 0)
         }
         else
         {
            orient = 9;   //  (1, 3, 0, 2)
         }
      else // test[0] == base[1]
         if (test[3] == base[2])
         {
            orient = 10;   //  (1, 0, 3, 2)
         }
         else
         {
            orient = 11;   //  (1, 0, 2, 3)
         }
   else if (test[2] == base[0])
      if (test[3] == base[1])
         if (test[0] == base[2])
         {
            orient = 12;   //  (2, 3, 0, 1)
         }
         else
         {
            orient = 13;   //  (2, 3, 1, 0)
         }
      else if (test[0] == base[1])
         if (test[1] == base[2])
         {
            orient = 14;   //  (2, 0, 1, 3)
         }
         else
         {
            orient = 15;   //  (2, 0, 3, 1)
         }
      else // test[1] == base[1]
         if (test[3] == base[2])
         {
            orient = 16;   //  (2, 1, 3, 0)
         }
         else
         {
            orient = 17;   //  (2, 1, 0, 3)
         }
   else // (test[3] == base[0])
      if (test[0] == base[1])
         if (test[2] == base[2])
         {
            orient = 18;   //  (3, 0, 2, 1)
         }
         else
         {
            orient = 19;   //  (3, 0, 1, 2)
         }
      else if (test[1] == base[1])
         if (test[0] == base[2])
         {
            orient = 20;   //  (3, 1, 0, 2)
         }
         else
         {
            orient = 21;   //  (3, 1, 2, 0)
         }
      else // test[2] == base[1]
         if (test[1] == base[2])
         {
            orient = 22;   //  (3, 2, 1, 0)
         }
         else
         {
            orient = 23;   //  (3, 2, 0, 1)
         }

#ifdef MFEM_DEBUG
   const int *aor = tet_t::Orient[orient];
   for (int j = 0; j < 4; j++)
      if (test[aor[j]] != base[j])
      {
         mfem_error("Mesh::GetTetOrientation(...)");
      }
#endif

   return orient;
}

int Mesh::CheckBdrElementOrientation(bool fix_it)
{
   int wo = 0; // count wrong orientations

   if (Dim == 2)
   {
      if (el_to_edge == NULL) // edges were not generated
      {
         el_to_edge = new Table;
         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
         GenerateFaces(); // 'Faces' in 2D refers to the edges
      }
      for (int i = 0; i < NumOfBdrElements; i++)
      {
         if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
         {
            int *bv = boundary[i]->GetVertices();
            int *fv = faces[be_to_edge[i]]->GetVertices();
            if (bv[0] != fv[0])
            {
               if (fix_it)
               {
                  mfem::Swap<int>(bv[0], bv[1]);
               }
               wo++;
            }
         }
      }
   }

   if (Dim == 3)
   {
      for (int i = 0; i < NumOfBdrElements; i++)
      {
         const int fi = be_to_face[i];

         if (faces_info[fi].Elem2No >= 0) { continue; }

         // boundary face
         int *bv = boundary[i]->GetVertices();
         // Make sure the 'faces' are generated:
         MFEM_ASSERT(fi < faces.Size(), "internal error");
         const int *fv = faces[fi]->GetVertices();
         int orientation; // orientation of the bdr. elem. w.r.t. the
         // corresponding face element (that's the base)
         const Element::Type bdr_type = GetBdrElementType(i);
         switch (bdr_type)
         {
            case Element::TRIANGLE:
            {
               orientation = GetTriOrientation(fv, bv);
               break;
            }
            case Element::QUADRILATERAL:
            {
               orientation = GetQuadOrientation(fv, bv);
               break;
            }
            default:
               MFEM_ABORT("Invalid 2D boundary element type \""
                          << bdr_type << "\"");
               orientation = 0; // suppress a warning
               break;
         }

         if (orientation % 2 == 0) { continue; }
         wo++;
         if (!fix_it) { continue; }

         switch (bdr_type)
         {
            case Element::TRIANGLE:
            {
               // swap vertices 0 and 1 so that we don't change the marked edge:
               // (0,1,2) -> (1,0,2)
               mfem::Swap<int>(bv[0], bv[1]);
               if (bel_to_edge)
               {
                  int *be = bel_to_edge->GetRow(i);
                  mfem::Swap<int>(be[1], be[2]);
               }
               break;
            }
            case Element::QUADRILATERAL:
            {
               mfem::Swap<int>(bv[0], bv[2]);
               if (bel_to_edge)
               {
                  int *be = bel_to_edge->GetRow(i);
                  mfem::Swap<int>(be[0], be[1]);
                  mfem::Swap<int>(be[2], be[3]);
               }
               break;
            }
            default: // unreachable
               break;
         }
      }
   }
   // #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
#ifdef MFEM_DEBUG
   if (wo > 0)
   {
      mfem::out << "Boundary elements with wrong orientation: " << wo << " / "
                << NumOfBdrElements << " (" << fixed_or_not[fix_it ? 0 : 1]
                << ")" << endl;
   }
#endif
   return wo;
}

int Mesh::GetNumGeometries(int dim) const
{
   MFEM_ASSERT(0 <= dim && dim <= Dim, "invalid dim: " << dim);
   int num_geoms = 0;
   for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
   {
      if (HasGeometry(Geometry::Type(g))) { num_geoms++; }
   }
   return num_geoms;
}

void Mesh::GetGeometries(int dim, Array<Geometry::Type> &el_geoms) const
{
   MFEM_ASSERT(0 <= dim && dim <= Dim, "invalid dim: " << dim);
   el_geoms.SetSize(0);
   for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
   {
      if (HasGeometry(Geometry::Type(g)))
      {
         el_geoms.Append(Geometry::Type(g));
      }
   }
}

void Mesh::GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const
{
   if (el_to_edge)
   {
      el_to_edge->GetRow(i, edges);
   }
   else
   {
      mfem_error("Mesh::GetElementEdges(...) element to edge table "
                 "is not generated.");
   }

   const int *v = elements[i]->GetVertices();
   const int ne = elements[i]->GetNEdges();
   cor.SetSize(ne);
   for (int j = 0; j < ne; j++)
   {
      const int *e = elements[i]->GetEdgeVertices(j);
      cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
   }
}

void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const
{
   if (Dim == 2)
   {
      edges.SetSize(1);
      cor.SetSize(1);
      edges[0] = be_to_edge[i];
      const int *v = boundary[i]->GetVertices();
      cor[0] = (v[0] < v[1]) ? (1) : (-1);
   }
   else if (Dim == 3)
   {
      if (bel_to_edge)
      {
         bel_to_edge->GetRow(i, edges);
      }
      else
      {
         mfem_error("Mesh::GetBdrElementEdges(...)");
      }

      const int *v = boundary[i]->GetVertices();
      const int ne = boundary[i]->GetNEdges();
      cor.SetSize(ne);
      for (int j = 0; j < ne; j++)
      {
         const int *e = boundary[i]->GetEdgeVertices(j);
         cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
      }
   }
}

void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
{
   if (Dim == 2)
   {
      edges.SetSize(1);
      edges[0] = i;
      o.SetSize(1);
      const int *v = faces[i]->GetVertices();
      o[0] = (v[0] < v[1]) ? (1) : (-1);
   }

   if (Dim != 3)
   {
      return;
   }

   GetFaceEdgeTable(); // generate face_edge Table (if not generated)

   face_edge->GetRow(i, edges);

   const int *v = faces[i]->GetVertices();
   const int ne = faces[i]->GetNEdges();
   o.SetSize(ne);
   for (int j = 0; j < ne; j++)
   {
      const int *e = faces[i]->GetEdgeVertices(j);
      o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
   }
}

void Mesh::GetEdgeVertices(int i, Array<int> &vert) const
{
   // the two vertices are sorted: vert[0] < vert[1]
   // this is consistent with the global edge orientation
   // generate edge_vertex Table (if not generated)
   if (!edge_vertex) { GetEdgeVertexTable(); }
   edge_vertex->GetRow(i, vert);
}

Table *Mesh::GetFaceEdgeTable() const
{
   if (face_edge)
   {
      return face_edge;
   }

   if (Dim != 3)
   {
      return NULL;
   }

#ifdef MFEM_DEBUG
   if (faces.Size() != NumOfFaces)
   {
      mfem_error("Mesh::GetFaceEdgeTable : faces were not generated!");
   }
#endif

   DSTable v_to_v(NumOfVertices);
   GetVertexToVertexTable(v_to_v);

   face_edge = new Table;
   GetElementArrayEdgeTable(faces, v_to_v, *face_edge);

   return (face_edge);
}

Table *Mesh::GetEdgeVertexTable() const
{
   if (edge_vertex)
   {
      return edge_vertex;
   }

   DSTable v_to_v(NumOfVertices);
   GetVertexToVertexTable(v_to_v);

   int nedges = v_to_v.NumberOfEntries();
   edge_vertex = new Table(nedges, 2);
   for (int i = 0; i < NumOfVertices; i++)
   {
      for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
      {
         int j = it.Index();
         edge_vertex->Push(j, i);
         edge_vertex->Push(j, it.Column());
      }
   }
   edge_vertex->Finalize();

   return edge_vertex;
}

Table *Mesh::GetVertexToElementTable()
{
   int i, j, nv, *v;

   Table *vert_elem = new Table;

   vert_elem->MakeI(NumOfVertices);

   for (i = 0; i < NumOfElements; i++)
   {
      nv = elements[i]->GetNVertices();
      v  = elements[i]->GetVertices();
      for (j = 0; j < nv; j++)
      {
         vert_elem->AddAColumnInRow(v[j]);
      }
   }

   vert_elem->MakeJ();

   for (i = 0; i < NumOfElements; i++)
   {
      nv = elements[i]->GetNVertices();
      v  = elements[i]->GetVertices();
      for (j = 0; j < nv; j++)
      {
         vert_elem->AddConnection(v[j], i);
      }
   }

   vert_elem->ShiftUpI();

   return vert_elem;
}

Table *Mesh::GetFaceToElementTable() const
{
   Table *face_elem = new Table;

   face_elem->MakeI(faces_info.Size());

   for (int i = 0; i < faces_info.Size(); i++)
   {
      if (faces_info[i].Elem2No >= 0)
      {
         face_elem->AddColumnsInRow(i, 2);
      }
      else
      {
         face_elem->AddAColumnInRow(i);
      }
   }

   face_elem->MakeJ();

   for (int i = 0; i < faces_info.Size(); i++)
   {
      face_elem->AddConnection(i, faces_info[i].Elem1No);
      if (faces_info[i].Elem2No >= 0)
      {
         face_elem->AddConnection(i, faces_info[i].Elem2No);
      }
   }

   face_elem->ShiftUpI();

   return face_elem;
}

void Mesh::GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const
{
   MFEM_VERIFY(el_to_face != NULL, "el_to_face not generated");

   el_to_face->GetRow(i, faces);

   int n = faces.Size();
   ori.SetSize(n);

   for (int j = 0; j < n; j++)
   {
      if (faces_info[faces[j]].Elem1No == i)
      {
         ori[j] = faces_info[faces[j]].Elem1Inf % 64;
      }
      else
      {
         MFEM_ASSERT(faces_info[faces[j]].Elem2No == i, "internal error");
         ori[j] = faces_info[faces[j]].Elem2Inf % 64;
      }
   }
}

void Mesh::GetBdrElementFace(int i, int *f, int *o) const
{
   const int *bv, *fv;

   *f = be_to_face[i];
   bv = boundary[i]->GetVertices();
   fv = faces[be_to_face[i]]->GetVertices();

   // find the orientation of the bdr. elem. w.r.t.
   // the corresponding face element (that's the base)
   switch (GetBdrElementType(i))
   {
      case Element::TRIANGLE:
         *o = GetTriOrientation(fv, bv);
         break;
      case Element::QUADRILATERAL:
         *o = GetQuadOrientation(fv, bv);
         break;
      default:
         MFEM_ABORT("invalid geometry");
   }
}

int Mesh::GetBdrElementEdgeIndex(int i) const
{
   switch (Dim)
   {
      case 1: return boundary[i]->GetVertices()[0];
      case 2: return be_to_edge[i];
      case 3: return be_to_face[i];
      default: MFEM_ABORT("invalid dimension!");
   }
   return -1;
}

void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
{
   int fid = GetBdrElementEdgeIndex(bdr_el);

   const FaceInfo &fi = faces_info[fid];
   MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0

   const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
   const int *bv = boundary[bdr_el]->GetVertices();
   int ori;
   switch (GetBdrElementGeometry(bdr_el))
   {
      case Geometry::POINT:    ori = 0; break;
      case Geometry::SEGMENT:  ori = (fv[0] == bv[0]) ? 0 : 1; break;
      case Geometry::TRIANGLE: ori = GetTriOrientation(fv, bv); break;
      case Geometry::SQUARE:   ori = GetQuadOrientation(fv, bv); break;
      default: MFEM_ABORT("boundary element type not implemented"); ori = 0;
   }
   el   = fi.Elem1No;
   info = fi.Elem1Inf + ori;
}

Element::Type Mesh::GetElementType(int i) const
{
   return elements[i]->GetType();
}

Element::Type Mesh::GetBdrElementType(int i) const
{
   return boundary[i]->GetType();
}

void Mesh::GetPointMatrix(int i, DenseMatrix &pointmat) const
{
   int k, j, nv;
   const int *v;

   v  = elements[i]->GetVertices();
   nv = elements[i]->GetNVertices();

   pointmat.SetSize(spaceDim, nv);
   for (k = 0; k < spaceDim; k++)
   {
      for (j = 0; j < nv; j++)
      {
         pointmat(k, j) = vertices[v[j]](k);
      }
   }
}

void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
{
   int k, j, nv;
   const int *v;

   v  = boundary[i]->GetVertices();
   nv = boundary[i]->GetNVertices();

   pointmat.SetSize(spaceDim, nv);
   for (k = 0; k < spaceDim; k++)
      for (j = 0; j < nv; j++)
      {
         pointmat(k, j) = vertices[v[j]](k);
      }
}

double Mesh::GetLength(int i, int j) const
{
   const double *vi = vertices[i]();
   const double *vj = vertices[j]();
   double length = 0.;

   for (int k = 0; k < spaceDim; k++)
   {
      length += (vi[k]-vj[k])*(vi[k]-vj[k]);
   }

   return sqrt(length);
}

// static method
void Mesh::GetElementArrayEdgeTable(const Array<Element*> &elem_array,
                                    const DSTable &v_to_v, Table &el_to_edge)
{
   el_to_edge.MakeI(elem_array.Size());
   for (int i = 0; i < elem_array.Size(); i++)
   {
      el_to_edge.AddColumnsInRow(i, elem_array[i]->GetNEdges());
   }
   el_to_edge.MakeJ();
   for (int i = 0; i < elem_array.Size(); i++)
   {
      const int *v = elem_array[i]->GetVertices();
      const int ne = elem_array[i]->GetNEdges();
      for (int j = 0; j < ne; j++)
      {
         const int *e = elem_array[i]->GetEdgeVertices(j);
         el_to_edge.AddConnection(i, v_to_v(v[e[0]], v[e[1]]));
      }
   }
   el_to_edge.ShiftUpI();
}

void Mesh::GetVertexToVertexTable(DSTable &v_to_v) const
{
   if (edge_vertex)
   {
      for (int i = 0; i < edge_vertex->Size(); i++)
      {
         const int *v = edge_vertex->GetRow(i);
         v_to_v.Push(v[0], v[1]);
      }
   }
   else
   {
      for (int i = 0; i < NumOfElements; i++)
      {
         const int *v = elements[i]->GetVertices();
         const int ne = elements[i]->GetNEdges();
         for (int j = 0; j < ne; j++)
         {
            const int *e = elements[i]->GetEdgeVertices(j);
            v_to_v.Push(v[e[0]], v[e[1]]);
         }
      }
   }
}

int Mesh::GetElementToEdgeTable(Table & e_to_f, Array<int> &be_to_f)
{
   int i, NumberOfEdges;

   DSTable v_to_v(NumOfVertices);
   GetVertexToVertexTable(v_to_v);

   NumberOfEdges = v_to_v.NumberOfEntries();

   // Fill the element to edge table
   GetElementArrayEdgeTable(elements, v_to_v, e_to_f);

   if (Dim == 2)
   {
      // Initialize the indices for the boundary elements.
      be_to_f.SetSize(NumOfBdrElements);
      for (i = 0; i < NumOfBdrElements; i++)
      {
         const int *v = boundary[i]->GetVertices();
         be_to_f[i] = v_to_v(v[0], v[1]);
      }
   }
   else if (Dim == 3)
   {
      if (bel_to_edge == NULL)
      {
         bel_to_edge = new Table;
      }
      GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
   }
   else
   {
      mfem_error("1D GetElementToEdgeTable is not yet implemented.");
   }

   // Return the number of edges
   return NumberOfEdges;
}

const Table & Mesh::ElementToElementTable()
{
   if (el_to_el)
   {
      return *el_to_el;
   }

   // Note that, for ParNCMeshes, faces_info will contain also the ghost faces
   MFEM_ASSERT(faces_info.Size() >= GetNumFaces(), "faces were not generated!");

   Array<Connection> conn;
   conn.Reserve(2*faces_info.Size());

   for (int i = 0; i < faces_info.Size(); i++)
   {
      const FaceInfo &fi = faces_info[i];
      if (fi.Elem2No >= 0)
      {
         conn.Append(Connection(fi.Elem1No, fi.Elem2No));
         conn.Append(Connection(fi.Elem2No, fi.Elem1No));
      }
      else if (fi.Elem2Inf >= 0)
      {
         int nbr_elem_idx = NumOfElements - 1 - fi.Elem2No;
         conn.Append(Connection(fi.Elem1No, nbr_elem_idx));
         conn.Append(Connection(nbr_elem_idx, fi.Elem1No));
      }
   }

   conn.Sort();
   conn.Unique();
   el_to_el = new Table(NumOfElements, conn);

   return *el_to_el;
}

const Table & Mesh::ElementToFaceTable() const
{
   if (el_to_face == NULL)
   {
      mfem_error("Mesh::ElementToFaceTable()");
   }
   return *el_to_face;
}

const Table & Mesh::ElementToEdgeTable() const
{
   if (el_to_edge == NULL)
   {
      mfem_error("Mesh::ElementToEdgeTable()");
   }
   return *el_to_edge;
}

void Mesh::AddPointFaceElement(int lf, int gf, int el)
{
   if (faces_info[gf].Elem1No == -1)  // this will be elem1
   {
      // faces[gf] = new Point(&gf);
      faces_info[gf].Elem1No  = el;
      faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
      faces_info[gf].Elem2No  = -1; // in case there's no other side
      faces_info[gf].Elem2Inf = -1; // face is not shared
   }
   else  //  this will be elem2
   {
      /* WARNING: Without the following check the mesh faces_info data structure
         may contain unreliable data. Normally, the order in which elements are
         processed could swap which elements appear as Elem1No and Elem2No. In
         branched meshes, where more than two elements can meet at a given node,
         the indices stored in Elem1No and Elem2No will be the first and last,
         respectively, elements found which touch a given node. This can lead to
         inconsistencies in any algorithms which rely on this data structure. To
         properly support branched meshes this data structure should be extended
         to support multiple elements per face. */
      /*
      MFEM_VERIFY(faces_info[gf].Elem2No < 0, "Invalid mesh topology. "
                  "Interior point found connecting 1D elements "
                  << faces_info[gf].Elem1No << ", " << faces_info[gf].Elem2No
                  << " and " << el << ".");
      */
      faces_info[gf].Elem2No  = el;
      faces_info[gf].Elem2Inf = 64 * lf + 1;
   }
}

void Mesh::AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
{
   if (faces[gf] == NULL)  // this will be elem1
   {
      faces[gf] = new Segment(v0, v1);
      faces_info[gf].Elem1No  = el;
      faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
      faces_info[gf].Elem2No  = -1; // in case there's no other side
      faces_info[gf].Elem2Inf = -1; // face is not shared
   }
   else  //  this will be elem2
   {
      MFEM_VERIFY(faces_info[gf].Elem2No < 0, "Invalid mesh topology.  "
                  "Interior edge found between 2D elements "
                  << faces_info[gf].Elem1No << ", " << faces_info[gf].Elem2No
                  << " and " << el << ".");
      int *v = faces[gf]->GetVertices();
      faces_info[gf].Elem2No  = el;
      if ( v[1] == v0 && v[0] == v1 )
      {
         faces_info[gf].Elem2Inf = 64 * lf + 1;
      }
      else if ( v[0] == v0 && v[1] == v1 )
      {
         // Temporarily allow even edge orientations: see the remark in
         // AddTriangleFaceElement().
         // Also, in a non-orientable surface mesh, the orientation will be even
         // for edges that connect elements with opposite orientations.
         faces_info[gf].Elem2Inf = 64 * lf;
      }
      else
      {
         MFEM_ABORT("internal error");
      }
   }
}

void Mesh::AddTriangleFaceElement(int lf, int gf, int el,
                                  int v0, int v1, int v2)
{
   if (faces[gf] == NULL)  // this will be elem1
   {
      faces[gf] = new Triangle(v0, v1, v2);
      faces_info[gf].Elem1No  = el;
      faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
      faces_info[gf].Elem2No  = -1; // in case there's no other side
      faces_info[gf].Elem2Inf = -1; // face is not shared
   }
   else  //  this will be elem2
   {
      MFEM_VERIFY(faces_info[gf].Elem2No < 0, "Invalid mesh topology.  "
                  "Interior triangular face found connecting elements "
                  << faces_info[gf].Elem1No << ", " << faces_info[gf].Elem2No
                  << " and " << el << ".");
      int orientation, vv[3] = { v0, v1, v2 };
      orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
      // In a valid mesh, we should have (orientation % 2 != 0), however, if
      // one of the adjacent elements has wrong orientation, both face
      // orientations can be even, until the element orientations are fixed.
      // MFEM_ASSERT(orientation % 2 != 0, "");
      faces_info[gf].Elem2No  = el;
      faces_info[gf].Elem2Inf = 64 * lf + orientation;
   }
}

void Mesh::AddQuadFaceElement(int lf, int gf, int el,
                              int v0, int v1, int v2, int v3)
{
   if (faces_info[gf].Elem1No < 0)  // this will be elem1
   {
      faces[gf] = new Quadrilateral(v0, v1, v2, v3);
      faces_info[gf].Elem1No  = el;
      faces_info[gf].Elem1Inf = 64 * lf; // face lf with orientation 0
      faces_info[gf].Elem2No  = -1; // in case there's no other side
      faces_info[gf].Elem2Inf = -1; // face is not shared
   }
   else  //  this will be elem2
   {
      MFEM_VERIFY(faces_info[gf].Elem2No < 0, "Invalid mesh topology.  "
                  "Interior quadrilateral face found connecting elements "
                  << faces_info[gf].Elem1No << ", " << faces_info[gf].Elem2No
                  << " and " << el << ".");
      int vv[4] = { v0, v1, v2, v3 };
      int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
      // Temporarily allow even face orientations: see the remark in
      // AddTriangleFaceElement().
      // MFEM_ASSERT(oo % 2 != 0, "");
      faces_info[gf].Elem2No  = el;
      faces_info[gf].Elem2Inf = 64 * lf + oo;
   }
}

void Mesh::GenerateFaces()
{
   int i, nfaces = GetNumFaces();

   for (i = 0; i < faces.Size(); i++)
   {
      FreeElement(faces[i]);
   }

   // (re)generate the interior faces and the info for them
   faces.SetSize(nfaces);
   faces_info.SetSize(nfaces);
   for (i = 0; i < nfaces; i++)
   {
      faces[i] = NULL;
      faces_info[i].Elem1No = -1;
      faces_info[i].NCFace = -1;
   }
   for (i = 0; i < NumOfElements; i++)
   {
      const int *v = elements[i]->GetVertices();
      const int *ef;
      if (Dim == 1)
      {
         AddPointFaceElement(0, v[0], i);
         AddPointFaceElement(1, v[1], i);
      }
      else if (Dim == 2)
      {
         ef = el_to_edge->GetRow(i);
         const int ne = elements[i]->GetNEdges();
         for (int j = 0; j < ne; j++)
         {
            const int *e = elements[i]->GetEdgeVertices(j);
            AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
         }
      }
      else
      {
         ef = el_to_face->GetRow(i);
         switch (GetElementType(i))
         {
            case Element::TETRAHEDRON:
            {
               for (int j = 0; j < 4; j++)
               {
                  const int *fv = tet_t::FaceVert[j];
                  AddTriangleFaceElement(j, ef[j], i,
                                         v[fv[0]], v[fv[1]], v[fv[2]]);
               }
               break;
            }
            case Element::WEDGE:
            {
               for (int j = 0; j < 2; j++)
               {
                  const int *fv = pri_t::FaceVert[j];
                  AddTriangleFaceElement(j, ef[j], i,
                                         v[fv[0]], v[fv[1]], v[fv[2]]);
               }
               for (int j = 2; j < 5; j++)
               {
                  const int *fv = pri_t::FaceVert[j];
                  AddQuadFaceElement(j, ef[j], i,
                                     v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
               }
               break;
            }
            case Element::HEXAHEDRON:
            {
               for (int j = 0; j < 6; j++)
               {
                  const int *fv = hex_t::FaceVert[j];
                  AddQuadFaceElement(j, ef[j], i,
                                     v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
               }
               break;
            }
            default:
               MFEM_ABORT("Unexpected type of Element.");
         }
      }
   }
}

void Mesh::GenerateNCFaceInfo()
{
   MFEM_VERIFY(ncmesh, "missing NCMesh.");

   for (int i = 0; i < faces_info.Size(); i++)
   {
      faces_info[i].NCFace = -1;
   }

   const NCMesh::NCList &list =
      (Dim == 2) ? ncmesh->GetEdgeList() : ncmesh->GetFaceList();

   nc_faces_info.SetSize(0);
   nc_faces_info.Reserve(list.masters.Size() + list.slaves.Size());

   int nfaces = GetNumFaces();

   // add records for master faces
   for (int i = 0; i < list.masters.Size(); i++)
   {
      const NCMesh::Master &master = list.masters[i];
      if (master.index >= nfaces) { continue; }

      faces_info[master.index].NCFace = nc_faces_info.Size();
      nc_faces_info.Append(NCFaceInfo(false, master.local, NULL));
      // NOTE: one of the unused members stores local face no. to be used below
   }

   // add records for slave faces
   for (int i = 0; i < list.slaves.Size(); i++)
   {
      const NCMesh::Slave &slave = list.slaves[i];

      if (slave.index < 0 || // degenerate slave face
          slave.index >= nfaces || // ghost slave
          slave.master >= nfaces) // has ghost master
      {
         continue;
      }

      FaceInfo &slave_fi = faces_info[slave.index];
      FaceInfo &master_fi = faces_info[slave.master];
      NCFaceInfo &master_nc = nc_faces_info[master_fi.NCFace];

      slave_fi.NCFace = nc_faces_info.Size();
      slave_fi.Elem2No = master_fi.Elem1No;
      slave_fi.Elem2Inf = 64 * master_nc.MasterFace; // get lf no. stored above
      // NOTE: In 3D, the orientation part of Elem2Inf is encoded in the point
      //       matrix. In 2D, the point matrix has the orientation of the parent
      //       edge, so its columns need to be flipped when applying it, see
      //       ApplyLocalSlaveTransformation.

      nc_faces_info.Append(
         NCFaceInfo(true, slave.master,
                    list.point_matrices[slave.geom][slave.matrix]));
   }
}

STable3D *Mesh::GetFacesTable()
{
   STable3D *faces_tbl = new STable3D(NumOfVertices);
   for (int i = 0; i < NumOfElements; i++)
   {
      const int *v = elements[i]->GetVertices();
      switch (GetElementType(i))
      {
         case Element::TETRAHEDRON:
         {
            for (int j = 0; j < 4; j++)
            {
               const int *fv = tet_t::FaceVert[j];
               faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
            }
            break;
         }
         case Element::WEDGE:
         {
            for (int j = 0; j < 2; j++)
            {
               const int *fv = pri_t::FaceVert[j];
               faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
            }
            for (int j = 2; j < 5; j++)
            {
               const int *fv = pri_t::FaceVert[j];
               faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
            }
            break;
         }
         case Element::HEXAHEDRON:
         {
            // find the face by the vertices with the smallest 3 numbers
            // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
            for (int j = 0; j < 6; j++)
            {
               const int *fv = hex_t::FaceVert[j];
               faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
            }
            break;
         }
         default:
            MFEM_ABORT("Unexpected type of Element.");
      }
   }
   return faces_tbl;
}

STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
{
   int i, *v;
   STable3D *faces_tbl;

   if (el_to_face != NULL)
   {
      delete el_to_face;
   }
   el_to_face = new Table(NumOfElements, 6);  // must be 6 for hexahedra
   faces_tbl = new STable3D(NumOfVertices);
   for (i = 0; i < NumOfElements; i++)
   {
      v = elements[i]->GetVertices();
      switch (GetElementType(i))
      {
         case Element::TETRAHEDRON:
         {
            for (int j = 0; j < 4; j++)
            {
               const int *fv = tet_t::FaceVert[j];
               el_to_face->Push(
                  i, faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]));
            }
            break;
         }
         case Element::WEDGE:
         {
            for (int j = 0; j < 2; j++)
            {
               const int *fv = pri_t::FaceVert[j];
               el_to_face->Push(
                  i, faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]));
            }
            for (int j = 2; j < 5; j++)
            {
               const int *fv = pri_t::FaceVert[j];
               el_to_face->Push(
                  i, faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
            }
            break;
         }
         case Element::HEXAHEDRON:
         {
            // find the face by the vertices with the smallest 3 numbers
            // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
            for (int j = 0; j < 6; j++)
            {
               const int *fv = hex_t::FaceVert[j];
               el_to_face->Push(
                  i, faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
            }
            break;
         }
         default:
            MFEM_ABORT("Unexpected type of Element.");
      }
   }
   el_to_face->Finalize();
   NumOfFaces = faces_tbl->NumberOfElements();
   be_to_face.SetSize(NumOfBdrElements);
   for (i = 0; i < NumOfBdrElements; i++)
   {
      v = boundary[i]->GetVertices();
      switch (GetBdrElementType(i))
      {
         case Element::TRIANGLE:
         {
            be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
            break;
         }
         case Element::QUADRILATERAL:
         {
            be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
            break;
         }
         default:
            MFEM_ABORT("Unexpected type of boundary Element.");
      }
   }

   if (ret_ftbl)
   {
      return faces_tbl;
   }
   delete faces_tbl;
   return NULL;
}

// shift cyclically 3 integers so that the smallest is first
static inline
void Rotate3(int &a, int &b, int &c)
{
   if (a < b)
   {
      if (a > c)
      {
         ShiftRight(a, b, c);
      }
   }
   else
   {
      if (b < c)
      {
         ShiftRight(c, b, a);
      }
      else
      {
         ShiftRight(a, b, c);
      }
   }
}

void Mesh::ReorientTetMesh()
{
   if (Dim != 3 || !(meshgen & 1))
   {
      return;
   }

   ResetLazyData();

   DSTable *old_v_to_v = NULL;
   Table *old_elem_vert = NULL;

   if (Nodes)
   {
      PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
   }

   for (int i = 0; i < NumOfElements; i++)
   {
      if (GetElementType(i) == Element::TETRAHEDRON)
      {
         int *v = elements[i]->GetVertices();

         Rotate3(v[0], v[1], v[2]);
         if (v[0] < v[3])
         {
            Rotate3(v[1], v[2], v[3]);
         }
         else
         {
            ShiftRight(v[0], v[1], v[3]);
         }
      }
   }

   for (int i = 0; i < NumOfBdrElements; i++)
   {
      if (GetBdrElementType(i) == Element::TRIANGLE)
      {
         int *v = boundary[i]->GetVertices();

         Rotate3(v[0], v[1], v[2]);
      }
   }

   if (!Nodes)
   {
      GetElementToFaceTable();
      GenerateFaces();
      if (el_to_edge)
      {
         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
      }
   }
   else
   {
      DoNodeReorder(old_v_to_v, old_elem_vert);
      delete old_elem_vert;
      delete old_v_to_v;
   }
}

int *Mesh::CartesianPartitioning(int nxyz[])
{
   int *partitioning;
   double pmin[3] = { infinity(), infinity(), infinity() };
   double pmax[3] = { -infinity(), -infinity(), -infinity() };
   // find a bounding box using the vertices
   for (int vi = 0; vi < NumOfVertices; vi++)
   {
      const double *p = vertices[vi]();
      for (int i = 0; i < spaceDim; i++)
      {
         if (p[i] < pmin[i]) { pmin[i] = p[i]; }
         if (p[i] > pmax[i]) { pmax[i] = p[i]; }
      }
   }

   partitioning = new int[NumOfElements];

   // determine the partitioning using the centers of the elements
   double ppt[3];
   Vector pt(ppt, spaceDim);
   for (int el = 0; el < NumOfElements; el++)
   {
      GetElementTransformation(el)->Transform(
         Geometries.GetCenter(GetElementBaseGeometry(el)), pt);
      int part = 0;
      for (int i = spaceDim-1; i >= 0; i--)
      {
         int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
         if (idx < 0) { idx = 0; }