Merge branch 'master' into meshing-additions

as suggested in GitHub issue #4612.

fem/bilininteg.hpp

@@ -2858,8 +2858,8 @@ public:
 };
 
 /** Integrator for $(Q u, v)$, where $Q$ is an optional coefficient (of type scalar,
-    vector (diagonal matrix), or matrix), trial function $u$ is in $H(curl$ or
-    $H(div)$, and test function $v$ is in $H(curl$, $H(div)$, or $v=(v_1,\dots,v_n)$, where
+    vector (diagonal matrix), or matrix), trial function $u$ is in $H(curl)$ or
+    $H(div)$, and test function $v$ is in $H(curl)$, $H(div)$, or $v=(v_1,\dots,v_n)$, where
     $v_i$ are in $H^1$. */
 class VectorFEMassIntegrator: public BilinearFormIntegrator
 {
@@ -3665,7 +3665,7 @@ public:
 /** Integrator for the form: $\langle v, w \times n \rangle$ over a face (the interface)
  *  In 3D the trial variable $v$ is defined on the interface ($H^{-1/2}$(curl), trace of $H(curl$)
  *  In 2D it's defined on the interface ($H^{1/2}$, trace of $H^1$)
- *  The test variable $w$ is in an $H(curl$-conforming space. */
+ *  The test variable $w$ is in an $H(curl)$-conforming space. */
 class TangentTraceIntegrator : public BilinearFormIntegrator
 {
 private:

fem/coefficient.cpp

@@ -262,6 +262,84 @@ void RestrictedCoefficient::SetTime(real_t t)
    this->Coefficient::SetTime(t);
 }
 
+real_t JacobianFunctionCoefficient::Eval(
+   ElementTransformation &T, const IntegrationPoint &ip)
+{
+   MFEM_VERIFY(&T.GetIntPoint() == &ip, "invalid input");
+   const DenseMatrix &JT = T.Jacobian();
+   if (!use_perf_J)
+   {
+      return JFunction(JT);
+   }
+   real_t JP_data[9];
+   DenseMatrix JP(JP_data, JT.Height(), JT.Width());
+   Geometries.JacToPerfJac(T.GetGeometryType(), JT, JP);
+   return JFunction(JP);
+}
+
+real_t MeshFunctionCoefficient::Eval(
+   ElementTransformation &T, const IntegrationPoint &ip)
+{
+   MFEM_VERIFY(&T.GetIntPoint() == &ip, "invalid input");
+   const DenseMatrix &JT = T.Jacobian();
+   real_t x_data[3];
+   Vector x(x_data, JT.Height());
+   T.Transform(ip, x);
+   if (!use_perf_J)
+   {
+      return XJFunction(x, JT);
+   }
+   real_t JP_data[9];
+   DenseMatrix JP(JP_data, JT.Height(), JT.Width());
+   Geometries.JacToPerfJac(T.GetGeometryType(), JT, JP);
+   return XJFunction(x, JP);
+}
+
+real_t JacobianDeterminantCoefficient::Eval(
+   ElementTransformation &T, const IntegrationPoint &ip)
+{
+   MFEM_VERIFY(&T.GetIntPoint() == &ip, "invalid input");
+   const real_t w = T.Weight();
+   if (!use_perf_J) { return w; }
+   const DenseMatrix *Jp = Geometries.GetPerfGeomToGeomJac(T.GetGeometryType());
+   return Jp ? w*Jp->Det() : w;
+}
+
+real_t MeshSizeCoefficient::Eval(
+   ElementTransformation &T, const IntegrationPoint &ip)
+{
+   MFEM_VERIFY(&T.GetIntPoint() == &ip, "invalid input");
+   const DenseMatrix &JT = T.Jacobian();
+   const int dim = JT.Width();
+   real_t JP_data[9];
+   DenseMatrix JP(JP_data, JT.Height(), dim);
+   Geometries.JacToPerfJac(T.GetGeometryType(), JT, JP);
+   if (dim == 1)
+   {
+      return JP.Weight();
+   }
+   if (type == 0)
+   {
+      const real_t w = JP.Weight();
+      return (w >= 0) ? std::pow(w, 1_r/dim) : -std::pow(-w, 1_r/dim);
+   }
+   if (dim == JP.Height())
+   {
+      if (type == 1) { return JP.CalcSingularvalue(dim-1); } // h_min
+      return JP.CalcSingularvalue(0); // h_max
+   }
+   if (dim == 2 && JP.Height() > 2) // suface mesh
+   {
+      real_t JtJ_data[4], e_val[2], e_vec[4];
+      DenseMatrix JtJ(JtJ_data, 2, 2);
+      MultAtB(JP, JP, JtJ);
+      JtJ.CalcEigenvalues(e_val, e_vec); // e_val are in increasing order
+      if (type == 1) { return std::sqrt(e_val[0]); } // h_min
+      return std::sqrt(e_val[1]); // h_max
+   }
+   MFEM_ABORT("unexpected Jacobian size: " << JP.Height() << " x " << dim);
+}
+
 void VectorCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                              const IntegrationRule &ir)
 {

fem/coefficient.hpp

@@ -562,6 +562,91 @@ public:
    { return active_attr[T.Attribute-1] ? c->Eval(T, ip, GetTime()) : 0.0; }
 };
 
+/** @brief Coefficient that evaluates a user specified scalar function of the
+    Jacobian of the mapping from the reference element (or the perfect element,
+    see Geometry::JacToPerfJac()) to the element described by the
+    ElementTransformation given to the Eval() method as input. */
+class JacobianFunctionCoefficient : public Coefficient
+{
+protected:
+   std::function<real_t(const DenseMatrix &)> JFunction;
+   bool use_perf_J;
+
+public:
+   /** @brief Construct a JacobianFunctionCoefficient that uses the perfect
+      Jacobian, when @a use_perf_J = true (default), or the reference Jacobian,
+      otherwise. */
+   JacobianFunctionCoefficient(std::function<real_t(const DenseMatrix &)> JF,
+                               bool use_perf_J = true)
+      : JFunction(std::move(JF)), use_perf_J(use_perf_J) { }
+
+   /// Evaluate the coefficient at the given point @a ip.
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+};
+
+/** @brief Coefficient that evaluates a user specified scalar function of the
+    mesh position and the Jacobian of the mapping from the reference element (or
+    the perfect element, see Geometry::JacToPerfJac()) to the element described
+    by the ElementTransformation given to the Eval() method as input. */
+class MeshFunctionCoefficient : public Coefficient
+{
+protected:
+   std::function<real_t(const Vector &, const DenseMatrix &)> XJFunction;
+   bool use_perf_J;
+
+public:
+   /** @brief Construct a MeshFunctionCoefficient that uses the perfect
+      Jacobian, when @a use_perf_J = true (default), or the reference Jacobian,
+      otherwise. */
+   MeshFunctionCoefficient(
+      std::function<real_t(const Vector &, const DenseMatrix &)> XJF,
+      bool use_perf_J = true)
+      : XJFunction(std::move(XJF)), use_perf_J(use_perf_J) { }
+
+   /// Evaluate the coefficient at the given point @a ip.
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+};
+
+/** @brief Coefficient that returns the Jacobian determinant (or integration
+    weight if the Jacobian is not square) of the mapping from the reference
+    element (or the perfect element, see Geometry::JacToPerfJac()) to the
+    element described by the ElementTransformation given to the Eval() method
+    as input. */
+class JacobianDeterminantCoefficient : public Coefficient
+{
+protected:
+   bool use_perf_J;
+
+public:
+   /** @brief Construct a JacobianDeterminantCoefficient that uses the perfect
+      Jacobian, when @a use_perf_J = true (default), or the reference Jacobian,
+      otherwise. */
+   JacobianDeterminantCoefficient(bool use_perf_J = true)
+      : use_perf_J(use_perf_J) { }
+
+   /// Evaluate the coefficient at the given point @a ip.
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+};
+
+/** @brief Coefficient that evaluates the pointwise mesh size as a function of
+    the Jacobian of the transformation from the perfect (unit) element. */
+class MeshSizeCoefficient : public Coefficient
+{
+protected:
+   int type;
+
+public:
+   /** @brief Construct a MeshSizeCoefficient of the given @a type. */
+   /** @param[in] type is one of:
+           - 0 - Jacobian determinant to the power 1/dim,
+           - 1 - h_min = minimal singular value of the Jacobian,
+           - 2 - h_max = maximal singular value of the Jacobian. */
+   MeshSizeCoefficient(int type = 0) : type(type) { }
+
+   /// Evaluate the coefficient at the given point @a ip.
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+};
+
 /// Base class for vector Coefficients that optionally depend on time and space.
 class VectorCoefficient
 {

fem/fe/fe_base.hpp

@@ -268,8 +268,8 @@ public:
       $ u(x) $ on a general physical element in following ways:
        - $ x = T(\hat x) $ is the image of the reference point $ \hat x $
        - $ J = J(\hat x) $ is the Jacobian matrix of the transformation T
-       - $ w = w(\hat x) = det(J) $ is the transformation weight factor for square J
-       - $ w = w(\hat x) = det(J^t J)^{1/2} $ is the transformation weight factor in general
+       - $ w = w(\hat x) = \det(J) $ is the transformation weight factor for square J
+       - $ w = w(\hat x) = \det(J^t J)^{1/2} $ is the transformation weight factor in general
    */
    enum MapType
    {

fem/fe_coll.cpp

@@ -858,6 +858,16 @@ const int *CubicFECollection::DofOrderForOrientation(Geometry::Type GeomType,
       };
       return sq_ind[Or];
    }
+   else if (GeomType == Geometry::CUBE)
+   {
+      // Note: only orientations 0 and 1 are supported!
+      static int hex_ord[2][8] =
+      {
+         {0, 1, 2, 3, 4, 5, 6, 7}, // orientation 0 = identity
+         {0, 3, 2, 1, 4, 7, 6, 5}  // orientation 1 = (x,y,z) -> (y,x,z)
+      };
+      return (0 <= Or && Or < 2) ? hex_ord[Or] : NULL;
+   }
 
    return NULL;
 }
@@ -1733,6 +1743,14 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
    {
       TetDofOrd[i] = NULL;
    }
+   for (int i = 0; i < 48; i++)
+   {
+      HexDofOrd[i] = NULL;
+   }
+   for (int i = 0; i < 12; i++)
+   {
+      PriDofOrd[i] = NULL;
+   }
 
    H1_dof[Geometry::POINT] = 1;
    H1_Elements[Geometry::POINT] = new PointFiniteElement;
@@ -1979,6 +1997,57 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
                }
             }
          }
+
+         const int HexDof = H1_dof[Geometry::CUBE];
+         if (HexDof > 0)
+         {
+            // Only orientations 0 and 1 will be defined!
+            HexDofOrd[0] = new int[2*HexDof];
+            for (int i = 1; i < 2; i++)
+            {
+               HexDofOrd[i] = HexDofOrd[i-1] + HexDof;
+            }
+            for (int iz = 0; iz < pm1; iz++)
+            {
+               for (int iy = 0; iy < pm1; iy++)
+               {
+                  for (int ix = 0; ix < pm1; ix++)
+                  {
+                     const int o0 = ix + pm1*(iy + pm1*iz);
+                     const int o1 = iy + pm1*(ix + pm1*iz);
+                     HexDofOrd[0][o0] = o0;
+                     HexDofOrd[1][o0] = o1;
+                  }
+               }
+            }
+         }
+
+         const int PriDof = H1_dof[Geometry::PRISM];
+         if (PriDof > 0)
+         {
+            // Only orientations 0 and 1 will be defined!
+            PriDofOrd[0] = new int[2*PriDof];
+            for (int i = 1; i < 2; i++)
+            {
+               PriDofOrd[i] = PriDofOrd[i-1] + PriDof;
+            }
+            for (int iz = 0; iz < pm1; iz++)
+            {
+               for (int iy = 0; iy < pm2; iy++)
+               {
+                  for (int ix = 0; ix + iy < pm2; ix++)
+                  {
+                     int t0 = TriDof - ((pm1 - iy)*(pm2 - iy))/2 + ix;
+                     int k = pm3 - iy - ix;
+                     int t1 = TriDof - ((pm1-iy)*(pm2-iy))/2 + k;  // (1,0,2)
+                     int o0 = t0 + TriDof*iz;
+                     int o1 = t1 + TriDof*iz;
+                     PriDofOrd[0][o0] = o0;
+                     PriDofOrd[1][o0] = o1;
+                  }
+               }
+            }
+         }
       }
    }
 }
@@ -2018,6 +2087,14 @@ const int *H1_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
    {
       return TetDofOrd[Or%24];
    }
+   else if (GeomType == Geometry::CUBE)
+   {
+      return HexDofOrd[Or%48];
+   }
+   else if (GeomType == Geometry::PRISM)
+   {
+      return PriDofOrd[Or%12];
+   }
    return NULL;
 }
 
@@ -2067,10 +2144,12 @@ const int *H1_FECollection::GetDofMap(Geometry::Type GeomType, int p) const
 
 H1_FECollection::~H1_FECollection()
 {
-   delete [] SegDofOrd[0];
-   delete [] TriDofOrd[0];
-   delete [] QuadDofOrd[0];
+   delete [] PriDofOrd[0];
+   delete [] HexDofOrd[0];
    delete [] TetDofOrd[0];
+   delete [] QuadDofOrd[0];
+   delete [] TriDofOrd[0];
+   delete [] SegDofOrd[0];
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
       delete H1_Elements[g];

fem/fe_coll.hpp

@@ -275,6 +275,7 @@ protected:
    FiniteElement *H1_Elements[Geometry::NumGeom];
    int H1_dof[Geometry::NumGeom];
    int *SegDofOrd[2], *TriDofOrd[6], *QuadDofOrd[8], *TetDofOrd[24];
+   int *HexDofOrd[48], *PriDofOrd[12];
 
 public:
    explicit H1_FECollection(const int p, const int dim = 3,

fem/geom.cpp

@@ -193,9 +193,9 @@ Geometry::Geometry()
    GeomCenter[PRISM].y = 1.0 / 3.0;
    GeomCenter[PRISM].z = 0.5;
 
-   GeomCenter[PYRAMID].x = 0.375;
-   GeomCenter[PYRAMID].y = 0.375;
-   GeomCenter[PYRAMID].z = 0.25;
+   GeomCenter[PYRAMID].x = 0.4;
+   GeomCenter[PYRAMID].y = 0.4;
+   GeomCenter[PYRAMID].z = 0.2;
 
    GeomToPerfGeomJac[POINT]       = NULL;
    GeomToPerfGeomJac[SEGMENT]     = new DenseMatrix(1);

fem/gridfunc.cpp

@@ -75,17 +75,18 @@ GridFunction::GridFunction(Mesh *m, std::istream &input)
    fes_sequence = fes->GetSequence();
 }
 
-GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
+GridFunction::GridFunction(Mesh *m, const GridFunction * const gf_array[],
+                           int num_pieces)
 {
    UseDevice(true);
 
    // all GridFunctions must have the same FE collection, vdim, ordering
    int vdim, ordering;
 
-   fes = gf_array[0]->FESpace();
-   fec_owned = FiniteElementCollection::New(fes->FEColl()->Name());
-   vdim = fes->GetVDim();
-   ordering = fes->GetOrdering();
+   const FiniteElementSpace *base_fes = gf_array[0]->FESpace();
+   fec_owned = FiniteElementCollection::New(base_fes->FEColl()->Name());
+   vdim = base_fes->GetVDim();
+   ordering = base_fes->GetOrdering();
    fes = new FiniteElementSpace(m, fec_owned, vdim, ordering);
    SetSize(fes->GetVSize());
 
@@ -104,7 +105,7 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
    vi = ei = fi = di = 0;
    for (int i = 0; i < num_pieces; i++)
    {
-      FiniteElementSpace *l_fes = gf_array[i]->FESpace();
+      const FiniteElementSpace *l_fes = gf_array[i]->FESpace();
       int l_ndofs  = l_fes->GetNDofs();
       int l_nvdofs = l_fes->GetNVDofs();
       int l_nedofs = l_fes->GetNEDofs();

fem/gridfunc.hpp

@@ -105,7 +105,7 @@ public:
        are owned by the GridFunction. */
    GridFunction(Mesh *m, std::istream &input);
 
-   GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces);
+   GridFunction(Mesh *m, const GridFunction * const gf_array[], int num_pieces);
 
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that

fem/tmop.cpp

@@ -3915,6 +3915,7 @@ real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
    for (int i = 0; i < ir.GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir.IntPoint(i);
+      if (Tpr) { Tpr->SetIntPoint(&ip); }
 
       metric->SetTargetJacobian(Jtr(i));
       CalcInverse(Jtr(i), Jrt);
@@ -4066,11 +4067,15 @@ real_t TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
          Mult(Jpr, Jrt, Jpt);
 
          real_t val = metric_normal * h_metric->EvalW(Jpt);
-         if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
+         if (metric_coeff)
+         {
+            Tpr->SetIntPoint(&ip);
+            val *= metric_coeff->Eval(*Tpr, ip);
+         }
 
          el_energy += weight * val;
-         delete Tpr;
       }
+      delete Tpr;
       energy += el_energy;
    }
    energy /= NEsplit;
@@ -4125,7 +4130,11 @@ real_t TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
       Mult(Jpr, Jrt, Jpt);
 
       real_t val = metric_normal * h_metric->EvalW(Jpt);
-      if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
+      if (metric_coeff)
+      {
+         Tpr->SetIntPoint(&ip);
+         val *= metric_coeff->Eval(*Tpr, ip);
+      }
 
       energy += weight * val;
    }
@@ -4236,6 +4245,7 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
    for (int q = 0; q < nqp; q++)
    {
       const IntegrationPoint &ip = ir.IntPoint(q);
+      if (Tpr) { Tpr->SetIntPoint(&ip); }
       metric->SetTargetJacobian(Jtr(q));
       CalcInverse(Jtr(q), Jrt);
       weights(q) = (integ_over_target) ? ip.weight * Jtr(q).Det() : ip.weight;

general/kdtree.hpp

@@ -287,13 +287,16 @@ public:
    void FindNeighborPoints(const PointND &pt,Tfloat R, std::vector<Tindex> & res,
                            std::vector<Tfloat> & dist)
    {
+      res.clear();
+      dist.clear();
       FindNeighborPoints(pt,R,data.begin(),data.end(),0,res,dist);
    }
 
    /// Finds all points within a distance R from point pt. The indices are
-   /// returned in the vector res and the correponding distances in vector dist.
+   /// returned in the vector res.
    void FindNeighborPoints(const PointND &pt,Tfloat R, std::vector<Tindex> & res)
    {
+      res.clear();
       FindNeighborPoints(pt,R,data.begin(),data.end(),0,res);
    }
 
@@ -302,6 +305,8 @@ public:
                                std::vector<Tindex> &res,
                                std::vector<Tfloat> &dist)
    {
+      res.clear();
+      dist.clear();
       Tfloat dd;
       for (auto iti=data.begin(); iti!=data.end(); iti++)
       {
@@ -318,6 +323,7 @@ public:
    void FindNeighborPointsSlow(const PointND &pt,Tfloat R,
                                std::vector<Tindex> &res)
    {
+      res.clear();
       Tfloat dd;
       for (auto iti=data.begin(); iti!=data.end(); iti++)
       {
@@ -537,7 +543,7 @@ private:
    }
 
    /// Finds the set of indices of points within a distance R of a point pt.
-   void FindNeighborPoints(PointND& pt, Tfloat R,
+   void FindNeighborPoints(const PointND& pt, Tfloat R,
                            typename std::vector<NodeND>::iterator itb,
                            typename std::vector<NodeND>::iterator ite,
                            size_t level,
@@ -586,7 +592,7 @@ private:
    }
 
    /// Finds the set of indices of points within a distance R of a point pt.
-   void FindNeighborPoints(PointND& pt, Tfloat R,
+   void FindNeighborPoints(const PointND& pt, Tfloat R,
                            typename std::vector<NodeND>::iterator itb,
                            typename std::vector<NodeND>::iterator ite,
                            size_t level,

linalg/densemat.cpp

@@ -1329,6 +1329,20 @@ void DenseMatrix::CalcEigenvalues(real_t *lambda, real_t *vec) const
    }
 }
 
+real_t DenseMatrix::CalcConditionNumber() const
+{
+   if (width == 1) { return 1_r; }
+   if (width == 2 && height > 2)
+   {
+      real_t AtA_data[4], e_val[2], e_vec[4];
+      DenseMatrix AtA(AtA_data, 2, 2);
+      MultAtB(*this, *this, AtA);
+      AtA.CalcEigenvalues(e_val, e_vec);
+      return std::sqrt(e_val[1]/e_val[0]);
+   }
+   return CalcSingularvalue(0)/CalcSingularvalue(width-1);
+}
+
 void DenseMatrix::GetRow(int r, Vector &row) const
 {
    int m = Height();

linalg/densemat.hpp

@@ -307,10 +307,15 @@ public:
    /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
    real_t CalcSingularvalue(const int i) const;
 
-   /** Return the eigenvalues (in increasing order) and eigenvectors of a
+   /** @brief Return the eigenvalues (in increasing order) and eigenvectors of a
        2x2 or 3x3 symmetric matrix. */
    void CalcEigenvalues(real_t *lambda, real_t *vec) const;
 
+   /** @brief Return the Euclidean, $\ell_2$, condition number of the matrix,
+       $\kappa = \sigma_{max}/\sigma_{min}$. Only small matrices are supported,
+       MxN with 1 <= N <= M <= 3. */
+   real_t CalcConditionNumber() const;
+
    void GetRow(int r, Vector &row) const;
    void GetColumn(int c, Vector &col) const;
    real_t *GetColumn(int col) { return data + col*height; }

linalg/kernels.hpp

@@ -543,7 +543,7 @@ MFEM_HOST_DEVICE static inline
 void GetScalingFactor(const real_t &d_max, real_t &mult)
 {
    int d_exp;
-   if (d_max > 0.)
+   if (d_max > 0_r)
    {
       mult = frexp(d_max, &d_exp);
       if (d_exp == std::numeric_limits<real_t>::max_exponent)
@@ -554,7 +554,7 @@ void GetScalingFactor(const real_t &d_max, real_t &mult)
    }
    else
    {
-      mult = 1.;
+      mult = 1_r;
    }
    // mult = 2^d_exp is such that d_max/mult is in [0.5,1) or in other words
    // d_max is in the interval [0.5,1)*mult

linalg/sparsemat.cpp

@@ -54,9 +54,9 @@ using namespace std;
 #ifdef MFEM_USE_CUDA_OR_HIP
 int SparseMatrix::SparseMatrixCount = 0;
 // doxygen doesn't like the macro-assisted typename so let's skip parsing it:
-// \cond false
+/// \cond Suppress_Doxygen_warnings
 MFEM_cu_or_hip(sparseHandle_t) SparseMatrix::handle = nullptr;
-// \endcond
+/// \endcond
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
 #endif // MFEM_USE_CUDA_OR_HIP

mesh/mesh.cpp

@@ -221,8 +221,7 @@ void Mesh::GetCharacteristics(real_t &h_min, real_t &h_max,
    {
       GetElementJacobian(i, J);
       h = pow(fabs(J.Weight()), 1.0/real_t(dim));
-      kappa = (dim == sdim) ?
-              J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1) : -1.0;
+      kappa = J.CalcConditionNumber();
       if (Vh) { (*Vh)(i) = h; }
       if (Vk) { (*Vk)(i) = kappa; }
 
@@ -330,6 +329,15 @@ void Mesh::PrintCharacteristics(Vector *Vh, Vector *Vk, std::ostream &os)
          << "kappa_min          : " << kappa_min << '\n'
          << "kappa_max          : " << kappa_max << '\n';
    }
+   os << "Mesh curvature     : ";
+   if (GetNodalFESpace() != NULL)
+   {
+      os << GetNodalFESpace()->FEColl()->Name() << '\n';
+   }
+   else
+   {
+      os << "NONE\n";
+   }
    os << '\n' << std::flush;
 }
 
@@ -3121,7 +3129,6 @@ void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
    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
@@ -3151,6 +3158,12 @@ void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
             case Geometry::TETRAHEDRON:
                new_or = GetTetOrientation(old_v, new_v);
                break;
+            case Geometry::CUBE:
+               new_or = GetHexOrientation(old_v, new_v);
+               break;
+            case Geometry::PRISM:
+               new_or = GetPrismOrientation(old_v, new_v);
+               break;
             default:
                new_or = 0;
                MFEM_ABORT(Geometry::Name[geom] << " elements (" << fec->Name()
@@ -3160,8 +3173,9 @@ void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
          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!");
+                     << "' does not define reordering for orientation "
+                     << new_or << " on " << Geometry::Name[geom]
+                     << " elements!");
          fes->GetElementInteriorDofs(i, old_dofs);
          new_dofs.SetSize(old_dofs.Size());
          for (int j = 0; j < new_dofs.Size(); j++)
@@ -3483,8 +3497,7 @@ void Mesh::Finalize(bool refine, bool fix_orientation)
    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))) );
+      ( check_orientation && fix_orientation );
 
    DSTable *old_v_to_v = NULL;
    Table *old_elem_vert = NULL;
@@ -5017,12 +5030,12 @@ void Mesh::Loader(std::istream &input, int generate_edges,
    // Finalize(...) should be called after this, if needed.
 }
 
-Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
+Mesh::Mesh(const Mesh * const mesh_array[], int num_pieces)
  : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
 {
    int      i, j, ie, ib, iv, *v, nv;
    Element *el;
-   Mesh    *m;
+   const Mesh *m;
 
    SetEmpty();
 
@@ -5142,10 +5155,10 @@ Mesh::Mesh(Mesh *mesh_array[], int num_pieces)
    FinalizeTopology();
 
    // copy the nodes (curvilinear meshes)
-   GridFunction *g = mesh_array[0]->GetNodes();
+   const GridFunction *g = mesh_array[0]->GetNodes();
    if (g)
    {
-      Array<GridFunction *> gf_array(num_pieces);
+      Array<const GridFunction *> gf_array(num_pieces);
       for (i = 0; i < num_pieces; i++)
       {
          gf_array[i] = mesh_array[i]->GetNodes();
@@ -6662,7 +6675,11 @@ int Mesh::CheckElementOrientation(bool fix_it)
                   wo++;
                   if (fix_it)
                   {
-                     // how?
+                     // on the reference element the following is a linear
+                     // transformation that has det(J) = -1
+                     mfem::Swap(vi[0], vi[1]);
+                     mfem::Swap(vi[3], vi[4]);
+                     fo++;
                   }
                }
                break;
@@ -6675,7 +6692,10 @@ int Mesh::CheckElementOrientation(bool fix_it)
                   wo++;
                   if (fix_it)
                   {
-                     // how?
+                     // on the reference element the following is a linear
+                     // transformation that has det(J) = -1
+                     mfem::Swap(vi[1], vi[3]);
+                     fo++;
                   }
                }
                break;
@@ -6688,7 +6708,11 @@ int Mesh::CheckElementOrientation(bool fix_it)
                   wo++;
                   if (fix_it)
                   {
-                     // how?
+                     // on the reference element the following is a linear
+                     // transformation that has det(J) = -1
+                     mfem::Swap(vi[1], vi[3]);
+                     mfem::Swap(vi[5], vi[7]);
+                     fo++;
                   }
                }
                break;
@@ -7014,6 +7038,76 @@ int Mesh::GetTetOrientation(const int *base, const int *test)
    return orient;
 }
 
+int Mesh::GetHexOrientation(const int *base, const int *test)
+{
+   // static method
+
+   // The hex orientations can be defined by permuting (x,y,z), which gives 6
+   // choices, followed by optional reflections about the planes x,y,z=1/2,
+   // which give 8 choices, for a total of 48 permutations.
+   const int hex_perm[2][8] =
+   {
+      { 0, 1, 2, 3, 4, 5, 6, 7 }, // 0: identity
+      { 0, 3, 2, 1, 4, 7, 6, 5 }  // 1: (x,y,z) -> (y,x,z), inverse of itself
+   };
+
+   for (int o = 0; o < 2; o++)
+   {
+      // check if test[hex_perm[o][j]] == base[j], for all j = 1,...,8
+      bool match = true;
+      for (int j = 0; j < 8; j++)
+      {
+         if (test[hex_perm[o][j]] != base[j])
+         {
+            match = false;
+            break;
+         }
+      }
+      if (match) { return o; }
+   }
+
+   MFEM_ABORT("unsupported orientation:\n\tbase: " << base[0] << ' ' << base[1]
+              << ' ' << base[2] << ' ' << base[3] << ' ' << base[4] << ' '
+              << base[5] << ' ' << base[6] << ' ' << base[7] << "\n\ttest: "
+              << test[0] << ' ' << test[1] << ' ' << test[2] << ' ' << test[3]
+              << ' ' << test[4] << ' ' << test[5] << ' ' << test[6] << ' '
+              << test[7]);
+}
+
+int Mesh::GetPrismOrientation(const int *base, const int *test)
+{
+   // static method
+
+   // The prism orientations can be defined by the 6 triangle orientations in
+   // the xy-plane, followed by an optional reflection about the plane z=1/2,
+   // which gives 2 choices, for a total of 12 permutations.
+   const int pri_perm[2][6] =
+   {
+      { 0, 1, 2, 3, 4, 5 }, // 0: identity
+      { 1, 0, 2, 4, 3, 5 }  // 1: (x,y,z) -> (1-x-y,y,z), inverse of itself
+   };
+
+   for (int o = 0; o < 2; o++)
+   {
+      // check if test[pri_perm[o][j]] == base[j], for all j = 1,...,6
+      bool match = true;
+      for (int j = 0; j < 6; j++)
+      {
+         if (test[pri_perm[o][j]] != base[j])
+         {
+            match = false;
+            break;
+         }
+      }
+      if (match) { return o; }
+   }
+
+   MFEM_ABORT("unsupported orientation:\n\tbase: " << base[0] << ' ' << base[1]
+              << ' ' << base[2] << ' ' << base[3] << ' ' << base[4] << ' '
+              << base[5] << "\n\ttest: " << test[0] << ' ' << test[1] << ' '
+              << test[2] << ' ' << test[3] << ' ' << test[4] << ' ' << test[5]);
+}
+
 int Mesh::CheckBdrElementOrientation(bool fix_it)
 {
    int wo = 0; // count wrong orientations
@@ -13129,7 +13223,7 @@ void Mesh::ScaleElements(real_t sf)
    delete [] vn;
 }
 
-void Mesh::Transform(void (*f)(const Vector&, Vector&))
+void Mesh::Transform(std::function<void(const Vector&, Vector&)> f)
 {
    // TODO: support for different new spaceDim.
    if (Nodes == NULL)
@@ -13142,13 +13236,13 @@ void Mesh::Transform(void (*f)(const Vector&, Vector&))
             vold(j) = vertices[i](j);
          }
          vnew.SetData(vertices[i]());
-         (*f)(vold, vnew);
+         f(vold, vnew);
       }
    }
    else
    {
       GridFunction xnew(Nodes->FESpace());
-      VectorFunctionCoefficient f_pert(spaceDim, f);
+      VectorFunctionCoefficient f_pert(spaceDim, std::move(f));
       xnew.ProjectCoefficient(f_pert);
       *Nodes = xnew;
    }
@@ -13364,6 +13458,512 @@ void Mesh::RemoveInternalBoundaries()
    attribs.Copy(bdr_attributes);
 }
 
+int Mesh::FindDuplicateVertices(int action,
+                                Array<int> &v2v,
+                                int &num_new_vertices,
+                                real_t abs_tol,
+                                real_t tol_mult) const
+{
+   const int verbosity_level = 0; // 0 or 1
+
+   if (abs_tol < 0) { return 1; }
+   if (tol_mult < 2) { return 2; }
+   if (action < 0 || action > 3) { return 3; }
+
+   const real_t real_tol = [&]()
+   {
+      if (abs_tol == 0_r) { return 0_r; }
+      int expon;
+      std::frexp(abs_tol, &expon);
+      return tol_mult*std::ldexp(1_r, expon);
+   }();
+
+   const real_t real_shift = abs_tol;
+   // const real_t real_shift = real_tol/tol_mult; // alternative; not tested
+
+   auto my_floor = [&](real_t val) -> real_t
+   {
+      if (real_tol == 0_r) { return val; }
+      real_t flr = real_tol*floor(val/real_tol);
+      if (!IsFinite(flr)) { flr = val; }
+      return flr;
+   };
+
+   auto hash_combine = [](std::size_t &seed, std::size_t h) -> void
+   {
+      // According to
+      // https://www.boost.org/doc/libs/1_85_0/libs/container_hash/doc/html/hash.html#notes_hash_combine
+      // the initial boost::hash_combine implementation used:
+      seed ^= h + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+   };
+
+   auto vertex_hash = [&](const Vertex &v) -> std::size_t
+   {
+      constexpr std::hash<real_t> hash_real;
+      std::size_t seed = 0;
+      hash_combine(seed, hash_real(v(0)));
+      hash_combine(seed, hash_real(v(1)));
+      hash_combine(seed, hash_real(v(2)));
+      return seed;
+   };
+
+   auto vertex_equal = [](const Vertex &a, const Vertex &b) -> bool
+   {
+      return a(0) == b(0) && a(1) == b(1) && a(2) == b(2);
+   };
+
+   // v2v is input when action == 3, so it is resized later.
+   if (action != 3)
+   {
+      v2v.SetSize(GetNV());
+   }
+
+   // construct vertex to active vertex map, v2av, and its inverse, av2v
+   Array<int> &v2av = v2v;
+   Array<int> av2v;
+   if (action == 1)
+   {
+      // the vertex subset is defined by the vertices of boundary elements
+      v2av = -1;
+      Array<int> bev;
+      for (int i = 0; i < GetNBE(); i++)
+      {
+         GetBdrElement(i)->GetVertices(bev);
+         for (int j = 0; j < bev.Size(); j++)
+         {
+            const int vj = bev[j];
+            if (v2av[vj] == -1)
+            {
+               v2av[vj] = av2v.Size();
+               av2v.Append(vj);
+            }
+         }
+      }
+   }
+   else if (action == 2)
+   {
+      // the vertex subset is defined by the vertices of boundary face elements
+      v2av = -1;
+      Array<int> fv;
+      const int num_faces = GetNumFaces(); // (dim-1) - dimensional entries
+      for (int i = 0; i < num_faces; i++)
+      {
+         // skip interior faces
+         if (faces_info[i].Elem2No >= 0) { continue; }
+         GetFace(i)->GetVertices(fv);
+         for (int j = 0; j < fv.Size(); j++)
+         {
+            const int vj = fv[j];
+            if (v2av[vj] == -1)
+            {
+               v2av[vj] = av2v.Size();
+               av2v.Append(vj);
+            }
+         }
+      }
+   }
+   else if (action == 3)
+   {
+      // v2v defines av2v. copy v2v to av2v while checking for invalid input:
+      av2v.SetSize(v2v.Size());
+      for (int ai = 0; ai < av2v.Size(); ai++)
+      {
+         const int i = v2v[ai]; // vertex 'i' is active vertex 'ai'
+         if (i < 0 || i >= GetNV()) { return 4; } // error: index out of range
+         av2v[ai] = i;
+      }
+      // define v2av:
+      v2v.SetSize(GetNV());
+      v2av = -1;
+      for (int ai = 0; ai < av2v.Size(); ai++)
+      {
+         const int i = av2v[ai]; // vertex 'i' is active vertex 'ai'
+         if (v2av[i] >= 0) { return 5; } // error: duplicate indices found
+         v2av[i] = ai;
+      }
+   }
+   const int num_active_vert = (action == 0) ? v2v.Size() : av2v.Size();
+   if (verbosity_level >= 1)
+   {
+      mfem::out << "number of all vertices:    " << v2v.Size() << endl;
+      mfem::out << "number of active vertices: " << num_active_vert << endl;
+   }
+
+   std::unordered_map<Vertex,int,decltype(vertex_hash),
+       decltype(vertex_equal)> av_map(0, vertex_hash, vertex_equal);
+   // av_map.max_load_factor(2.0f);  // slows down things a bit
+   Array<int> av2cube(num_active_vert);
+   const int sdim = SpaceDimension();
+   for (int i = 0; i < num_active_vert; i++)
+   {
+      Vertex v(0_r, 0_r, 0_r);
+      const int vi = (action == 0) ? i : av2v[i];
+      v.SetCoords(sdim, GetVertex(vi));
+      for (int j = 0; j < sdim; j++)
+      {
+         v(j) = my_floor(v(j));
+      }
+      auto res = av_map.emplace(v, (int)av_map.size());
+      if (!res.second)
+      {
+         av2cube[i] = res.first->second;
+      }
+      else
+      {
+         av2cube[i] = (int)av_map.size()-1;
+      }
+   }
+   if (verbosity_level >= 1)
+   {
+      auto num_buckets = av_map.bucket_count();
+      mfem::out << "av_map.size()         = " << av_map.size() << endl;
+      mfem::out << "av_map.bucket_count() = " << num_buckets << endl;
+      mfem::out << "av_map.load_factor()  = " << av_map.load_factor() << endl;
+      std::size_t max_bucket_size = 0, num_buckets_other_sizes = 0;
+      const int N = 9;
+      std::size_t num_buckets_by_size[N] = {};
+      for (std::size_t i = 0; i < num_buckets; i++)
+      {
+         auto bucket_size = av_map.bucket_size(i);
+         if (bucket_size < N) { num_buckets_by_size[bucket_size]++; }
+         else { num_buckets_other_sizes++; }
+         max_bucket_size = std::max(max_bucket_size, bucket_size);
+      }
+      mfem::out << "max bucket size   = " << max_bucket_size << endl;
+      mfem::out << "num buckets by size:" << endl;
+      for (int s = 0; s < N; s++)
+      {
+         auto nb = num_buckets_by_size[s];
+         mfem::out << "size " << s << " = " << nb
+                   << '\t' << real_t(nb)/num_buckets*100 << '%' << endl;
+      }
+      mfem::out << "rest   = " << num_buckets_other_sizes << '\t'
+                << real_t(num_buckets_other_sizes)/num_buckets*100 << '%'
+                << endl;
+   }
+   Table cube2av;
+   // Note: the column (av) indices in each row (cube) of cube2av are sorted.
+   Transpose(av2cube, cube2av, (int)av_map.size());
+   if (verbosity_level >= 1)
+   {
+      // print 'cube2av' stats
+      mfem::out << "number of cubes: " << cube2av.Size() << endl;
+      mfem::out << "average number of active vertices / cube: "
+                << real_t(num_active_vert)/cube2av.Size() << endl;
+      int max_cube_size = 0;
+      for (int c = 0; c < cube2av.Size(); c++)
+      {
+         max_cube_size = std::max(max_cube_size, cube2av.RowSize(c));
+      }
+      mfem::out << "maximum number of active vertices / cube: "
+                << max_cube_size << endl;
+   }
+
+   // some statistics variabes
+   int num_dist_evals = 0;
+   real_t max_dist_sq = 0_r;
+
+   // store info about duplicates in av2av: av2av[i] != i means active
+   // vertex i is duplicate of active vertex av2av[i]
+   Array<int> av2av(num_active_vert);
+   // initialize av2av with identity
+   for (int i = 0; i < av2av.Size(); i++)
+   {
+      av2av[i] = i;
+   }
+   // loop over the active vertices
+   Vector vi_p(sdim);
+   for (int i = 0; i < num_active_vert; i++)
+   {
+      const int vi = (action == 0) ? i : av2v[i];
+      vi_p = GetVertex(vi);
+      // 1. find the min index i_p such that the distance from active
+      //    vertex i_p to active vertex i is <= abs_tol
+      int i_p = i;
+      // 1.1. check the cube containg vi_p, i.e. active vertex i
+      int cube = av2cube[i];
+      int cube_size = cube2av.RowSize(cube);
+      const int *cube_avs = cube2av.GetRow(cube);
+      for (int j = 0; true; j++)
+      {
+         const int av_j = cube_avs[j];
+         if (av_j >= i_p) { break; }
+         const int vj = (action == 0) ? av_j : av2v[av_j];
+         const real_t dist_sq = vi_p.DistanceSquaredTo(GetVertex(vj));
+         num_dist_evals++;
+         if (dist_sq <= abs_tol*abs_tol)
+         {
+            max_dist_sq = std::fmax(max_dist_sq, dist_sq);
+            i_p = av_j;
+            break;
+         }
+      }
+      // 1.2. check other adjacent cubes, if needed
+      int shift[3] = { 0, 0, 0 };
+      Vertex vi_v(0_r, 0_r, 0_r);
+      for (int j = 0; j < sdim; j++)
+      {
+         vi_v(j) = my_floor(vi_p(j));
+         if (vi_v(j) != my_floor(vi_p(j)-real_shift))
+         {
+            shift[j] = -1;
+         }
+         if (vi_v(j) != my_floor(vi_p(j)+real_shift))
+         {
+            MFEM_VERIFY(shift[j] == 0, "internal error");
+            shift[j] = +1;
+         }
+      }
+      // loop over the adjacent cubes -- there are up to 7 of those that
+      // may need to be checked
+      for (int cb = 1; cb < 8; cb++)
+      {
+         const int c[3] = { cb%2, (cb/2)%2, cb/4 };
+         int sh[3];
+         bool skip = false;
+         for (int j = 0; j < 3; j++)
+         {
+            if (c[j] != 0 && shift[j] == 0) { skip = true; break; }
+            sh[j] = (c[j] == 0) ? 0 : shift[j];
+         }
+         if (skip) { continue; }
+         for (int j = 0; j < sdim; j++)
+         {
+            vi_v(j) = my_floor(vi_p(j)+sh[j]*real_shift);
+         }
+         auto av_map_iter = av_map.find(vi_v);
+         if (av_map_iter == av_map.end()) { continue; }
+         cube = av_map_iter->second;
+         cube_size = cube2av.RowSize(cube);
+         cube_avs = cube2av.GetRow(cube);
+         for (int j = 0; j < cube_size; j++)
+         {
+            const int av_j = cube_avs[j];
+            if (av_j >= i_p) { break; }
+            const int vj = (action == 0) ? av_j : av2v[av_j];
+            const real_t dist_sq = vi_p.DistanceSquaredTo(GetVertex(vj));
+            num_dist_evals++;
+            if (dist_sq <= abs_tol*abs_tol)
+            {
+               max_dist_sq = std::fmax(max_dist_sq, dist_sq);
+               i_p = av_j;
+               break;
+            }
+         }
+      }
+      // 2. i_p is set; if i_p != i, i.e. i_p < i, update av2av to say that i
+      //    is duplicate of i_p
+      if (i_p != i)
+      {
+         av2av[i] = i_p;
+      }
+   }
+   int ret_code = 0;
+   int num_unresolved = 0;
+   for (int i = 0; i < av2av.Size(); i++)
+   {
+      const int p = av2av[i];
+      if (i != p && p != av2av[p]) { num_unresolved++; }
+   }
+   if (num_unresolved != 0)
+   {
+      ret_code = 10; // warning: transitive duplicates found
+   }
+   // using v2av (which is a ref to v2v), av2av and av2v, construct v2v
+   int num_new_vert = 0;
+   for (int i = 0; i < v2v.Size(); i++)
+   {
+      if (action != 0 && v2v[i] < 0) // 'i' is not an active vertex
+      {
+         v2v[i] = num_new_vert++;
+      }
+      else
+      {
+         // vertex 'i' is active vertex 'ai'
+         const int ai = (action == 0) ? i : v2v[i];
+         if (av2av[ai] == ai) // 'ai' is not a duplicate
+         {
+            v2v[i] = num_new_vert++;
+         }
+         else
+         {
+            // vertex 'i' is active vertex 'ai' which is a duplicate; we
+            // will assign it a new vertex id later when all non-duplicate
+            // vertices have been assigned a new vertex id
+            v2v[i] = -1;
+         }
+      }
+   }
+   // assign new vertex ids to duplicate active vertices
+   for (int ai = 0; ai < num_active_vert; ai++)
+   {
+      const int ai_p = av2av[ai]; // by construction, we have: ai_p <= ai
+      if (ai_p != ai)
+      {
+         // active vertex 'ai' is a duplicate of active vertex 'ai_p'
+         // active vertex 'ai' is vertex 'i'
+         const int i = (action == 0) ? ai : av2v[ai];
+         // active vertex 'ai_p' is vertex 'i_p'
+         const int i_p = (action == 0) ? ai_p : av2v[ai_p];
+         // since ai_p < ai, i_p has already been assigned a new vertex id in
+         // v2v -- we assign the same new vertex id to 'i':
+         v2v[i] = v2v[i_p];
+      }
+   }
+
+   num_new_vertices = num_new_vert;
+
+   if (verbosity_level >= 1)
+   {
+      mfem::out << "num removed duplicate vertices: "
+                << GetNV()-num_new_vert << endl;
+      mfem::out << "num distance evals: " << num_dist_evals << endl;
+      mfem::out << "max measured distance between duplicate vertices: "
+                << std::sqrt(max_dist_sq) << endl;
+   }
+
+   return ret_code;
+}
+
+int Mesh::ApplyVertexMap(const Array<int> &v2v, int num_new_vertices)
+{
+   Array<int> ev;
+   auto have_repetitions = [&](const int *v, int nv) -> bool
+   {
+      ev.SetSize(nv);
+      ev.CopyFrom(v);
+      ev.Sort();
+      // return true if there are repetitions
+      return std::adjacent_find(ev.begin(), ev.end()) != ev.end();
+   };
+
+   Mesh mesh_new(Dimension(), num_new_vertices, GetNE(), GetNBE(),
+                 SpaceDimension());
+   // Add elements with reassigned vertex indices
+   for (int i = 0; i < GetNE(); i++)
+   {
+      Element *el = GetElement(i)->Duplicate(&mesh_new);
+      int *v = el->GetVertices();
+      int nv = el->GetNVertices();
+      for (int j = 0; j < nv; j++)
+      {
+         v[j] = v2v[v[j]];
+      }
+      mesh_new.AddElement(el);
+      if (have_repetitions(v, nv))
+      {
+         return 1;
+      }
+   }
+   // Add boundary elements with reassigned vertex indices
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      Element *el = GetBdrElement(i)->Duplicate(&mesh_new);
+      int *v = el->GetVertices();
+      int nv = el->GetNVertices();
+      for (int j = 0; j < nv; j++)
+      {
+         v[j] = v2v[v[j]];
+      }
+      mesh_new.AddBdrElement(el);
+      if (have_repetitions(v, nv))
+      {
+         return 1;
+      }
+   }
+   // New vertex coordinates are averages of the original ones.
+   for (int i = 0; i < num_new_vertices; i++)
+   {
+      mesh_new.AddVertex(0_r, 0_r, 0_r);
+   }
+   const int sdim = SpaceDimension();
+   Array<int> new_vertex_counter(num_new_vertices);
+   new_vertex_counter = 0;
+   for (int i = 0; i < GetNV(); i++)
+   {
+      const real_t *vtx_orig = GetVertex(i);
+      real_t *vtx_new = mesh_new.GetVertex(v2v[i]);
+      for (int j = 0; j < sdim; j++)
+      {
+         vtx_new[j] += vtx_orig[j];
+      }
+      new_vertex_counter[v2v[i]]++;
+   }
+   for (int i = 0; i < num_new_vertices; i++)
+   {
+      int vtx_count = new_vertex_counter[i];
+      if (vtx_count == 0) { continue; }
+      real_t *vtx_new = mesh_new.GetVertex(i);
+      for (int j = 0; j < sdim; j++)
+      {
+         vtx_new[j] /= vtx_count;
+      }
+   }
+   new_vertex_counter.DeleteAll();
+   mesh_new.FinalizeTopology(/* generate_bdr: */ false);
+   if (Nodes)
+   {
+      FiniteElementCollection *nodal_fec_new =
+         FiniteElementCollection::New(Nodes->FESpace()->FEColl()->Name());
+      FiniteElementSpace *nodal_fes_new =
+         new FiniteElementSpace(&mesh_new, nodal_fec_new, sdim,
+                                Nodes->FESpace()->GetOrdering());
+      if (Nodes->FESpace()->IsVariableOrder())
+      {
+         for (int i = 0; i < GetNE(); i++)
+         {
+            auto p = Nodes->FESpace()->GetElementOrder(i);
+            nodal_fes_new->SetElementOrder(i, p);
+         }
+         nodal_fes_new->Update(false);
+      }
+      GridFunction *nodes_new = new GridFunction(nodal_fes_new);
+      // give ownership of nodal_fec_new and nodal_fes_new to nodes_new:
+      nodes_new->MakeOwner(nodal_fec_new);
+      // the new mesh_nodes are averages of the old ones
+      Array<int> new_vdof_counter(nodes_new->Size());
+      new_vdof_counter = 0;
+      *nodes_new = 0_r;
+      Array<int> vdofs;
+      Vector loc_nodes;
+      for (int i = 0; i < GetNE(); i++)
+      {
+         Nodes->FESpace()->GetElementVDofs(i, vdofs);
+         Nodes->GetSubVector(vdofs, loc_nodes);
+         nodal_fes_new->GetElementVDofs(i, vdofs);
+         // Verify that the number of vdofs on this element in the original and
+         // in the new meshes is the same.
+         MFEM_VERIFY(loc_nodes.Size() == vdofs.Size(), "internal error");
+         nodes_new->AddElementVector(vdofs, loc_nodes);
+         for (int j = 0; j < vdofs.Size(); j++)
+         {
+            new_vdof_counter[vdofs[j]]++;
+         }
+      }
+      for (int i = 0; i < new_vdof_counter.Size(); i++)
+      {
+         const int vd_count = new_vdof_counter[i];
+         if (vd_count > 0)
+         {
+            (*nodes_new)[i] /= vd_count;
+         }
+      }
+      new_vdof_counter.DeleteAll();
+      const Operator *Pconf = nodal_fes_new->GetProlongationMatrix();
+      if (Pconf)
+      {
+         const SparseMatrix *Rconf = nodal_fes_new->GetRestrictionMatrix();
+         Vector t_nodes_new(Rconf->Height());
+         Rconf->Mult(*nodes_new, t_nodes_new);
+         Pconf->Mult(t_nodes_new, *nodes_new);
+      }
+      mesh_new.NewNodes(*nodes_new, true);
+   }
+   Swap(mesh_new, true);
+
+   return 0;
+}
+
 void Mesh::FreeElement(Element *E)
 {
 #ifdef MFEM_USE_MEMALLOC

mesh/mesh.hpp

@@ -544,6 +544,16 @@ protected:
    /// Returns the orientation of "test" relative to "base"
    static int GetTetOrientation(const int *base, const int *test);
 
+   /// Returns the orientation of "test" relative to "base"
+   /** @warning For now, only a minimal set of orientations is supported - the
+       method with generate an MFEM error for unsupported orientations. */
+   static int GetHexOrientation(const int *base, const int *test);
+
+   /// Returns the orientation of "test" relative to "base"
+   /** @warning For now, only a minimal set of orientations is supported - the
+       method with generate an MFEM error for unsupported orientations. */
+   static int GetPrismOrientation(const int *base, const int *test);
+
    static void GetElementArrayEdgeTable(const Array<Element*> &elem_array,
                                         const DSTable &v_to_v,
                                         Table &el_to_edge);
@@ -693,7 +703,10 @@ public:
 
    /// Construct a Mesh from the given primary data.
    /** The array @a vertices is used as external data, i.e. the Mesh does not
-       copy the data and will not delete the pointer.
+       copy the data and will not delete the pointer. Note that this array needs
+       to have three components (x,y,z) for the vertex coordinates regardless
+       of @a dimension and @a space_dimension; the extra components (e.g. the z
+       component for a mesh in 2D) should be initialized to zero.
 
        The data from the other arrays is copied into the internal Mesh data
        structures.
@@ -738,7 +751,7 @@ public:
    /// Create a disjoint mesh from the given mesh array
    ///
    /// @note Data is copied from the meshes in @a mesh_array.
-   Mesh(Mesh *mesh_array[], int num_pieces);
+   Mesh(const Mesh * const mesh_array[], int num_pieces);
 
    /** This is similar to the mesh constructor with the same arguments, but here
        the current mesh is destroyed and another one created based on the data
@@ -893,6 +906,15 @@ public:
        SetCurvature() for further details. */
    static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector<int> &v2v);
 
+   /// Create a disjoint mesh from the given meshes
+   static Mesh MakeUnion(const Mesh * const meshes[], int num_meshes)
+   { return Mesh(meshes, num_meshes); }
+
+   /// Create a disjoint mesh from the given meshes
+   template <int num_meshes>
+   static Mesh MakeUnion(const Mesh * const (&meshes)[num_meshes])
+   { return Mesh(meshes, num_meshes); }
+
    ///@}
 
    /// Construct a Mesh from a NURBSExtension
@@ -1080,8 +1102,10 @@ public:
    virtual void SetAttributes();
 
    /// Check (and optionally attempt to fix) the orientation of the elements
-   /** @param[in] fix_it  If `true`, attempt to fix the orientations of some
-                          elements: triangles, quads, and tets.
+   /** @param[in] fix_it  If `true`, attempt to fix the orientations of the
+                          elements: this operation will succeed only when the
+                          determinant of the Jacobian does not change its sign
+                          throughout the element.
        @return The number of elements with wrong orientation.
 
        @note For meshes with nodes (e.g. high-order or periodic meshes), fixing
@@ -1116,6 +1140,101 @@ public:
        have two adjacent faces in 3D, or edges in 2D. */
    void RemoveInternalBoundaries();
 
+   /// Helper method to set the internal vertex coordinates from the mesh nodes.
+   void SetVerticesFromNodes() { if (Nodes) { SetVerticesFromNodes(Nodes); } }
+
+   /** @brief Find duplicate vertices, in a subset of all vertices, based on
+       their physical location. */
+   /** Here, we call the vertices in the subset that will be searched for
+       duplicates, active vertices.
+       @param[in] action
+          The following action codes are valid:
+            - 0: Find duplicates among all vertices in the mesh.
+            - 1: Find duplicates among the boundary vertices in the mesh. For
+              this action, vertices are considered to be boundary if they are
+              used by at least one boundary element.
+            - 2: Find duplicates among the boundary vertices in the mesh. For
+              this action, vertices are considered to be boundary if they are
+              used by at least one boundary-face element, i.e. face that has
+              exactly one adjacent mesh element.
+            - 3: The input parameter @a v2v defines the set of active vertices,
+              i.e. the set of vertices that will be searched for duplicates. See
+              the description of @a v2v for details.
+       @param[in,out] v2v
+          This parameter is input parameter only when @a action is 3. In this
+          case, its size must be the number of active vertices and `v2v[ai]`
+          must be the vertex index for active vertex `ai`. The entries contained
+          in @a v2v must be in the range [0,GetNV()) and there can be no
+          repeated entries. At exit, this parameter defines a map from the
+          current vertex indices to new vertex indices where duplicates have
+          been removed. In case of an error, this parameter may remain
+          unmodified, see the description of the return values.
+       @param[out] num_new_vertices
+          Set to the number of new vertices, i.e. the largest entity in @a v2v
+          plus 1. In case of an error, this parameter may remain unmodified, see
+          the description of the return values.
+       @param[in] abs_tol
+          Absolute distance below which active vertices will be considered
+          duplicates, must be >= 0.
+       @param[in] tol_mult
+          Tolerance multiplier, must be >= 2; larger values may speed up the
+          algorithm, so the default value is set to 16; however, for large
+          @a abs_tol (close to the minimum distance between non-duplicate
+          vertices) smaller values of @a tol_mult may be faster.
+       @returns
+          One of the following status codes is returned:
+            - 0: Success: no warnings or errors.
+            - 1: Error: invalid (negative) value for @a abs_tol; @a v2v and
+              @a num_new_vertices are not modified.
+            - 2: Error: invalid (less than 2) value for @a tol_mult; @a v2v and
+              @a num_new_vertices are not modified.
+            - 3: Error: invalid @a action parameter; @a v2v and
+              @a num_new_vertices are not modified.
+            - 4: Error: @a action is 3 and @a v2v contains an index outside the
+              valid range [0,GetNV()); @a v2v and @a num_new_vertices are not
+              modified.
+            - 5: Error: @a action is 3 and @a v2v contains repeated indices;
+              @a v2v is modified and @a num_new_vertices is not modified.
+            - 10: Warning: transitive duplicates were found which means that
+              there are vertices a, b, and c such that dist(a,b) <= abs_tol,
+              dist(b,c) <= abs_tol, however dist(a,c) > abs_tol. In such cases,
+              vertices a and c are marked as duplicates even though they do not
+              meet the @a abs_tol requirement. Typically, this warning can be
+              resolved by increasing @a abs_tol. Both output parameters, @a v2v
+              and @a num_new_vertices are set as in the case of success.
+
+       On success, the output parameters @a v2v and @a num_new_vertices can be
+       used directly as input for the method ApplyVertexMap() to remove the
+       duplicate vertices from the mesh connecting topologically disconnected
+       pieces.
+
+       @note This method uses the internal vertex coordinates to search for
+       duplicates, so for meshes with nodes, it maybe necessary to call the
+       method SetVerticesFromNodes() before calling this method. */
+   int FindDuplicateVertices(int action,
+                             Array<int> &v2v,
+                             int &num_new_vertices,
+                             real_t abs_tol,
+                             real_t tol_mult = 16) const;
+
+   /// Apply the @a v2v map to re-enumerate the mesh vertices.
+   /** @param[in] v2v
+          The size of this array must be GetNV(). The entries `v2v[i]` must be
+          in the range [0, @a num_new_verices). The map can add new (unused)
+          vertices, combine multiple existing vertices into a single vertex
+          (e.g. removing duplicate vertices or introducing periodicity), or
+          simply permute the existing vertices.
+       @param[in] num_new_vertices
+          The number of vertices in the new mesh. Entries in @a v2v are expected
+          to be in the range [0, @a num_new_vertices).
+       @returns
+          One of the following status codes is returned:
+            - 0: Success.
+            - 1: The mapping @a v2v is invalid: if the mapping is applied as
+              given it will result in mesh elements or boundary elements with
+              repeated vertex indices which is not allowed. */
+   int ApplyVertexMap(const Array<int> &v2v, int num_new_vertices);
+
    /**
     * @brief Clear the boundary element to edge map.
     */
@@ -1173,7 +1292,7 @@ public:
    /// @ref mfem_Mesh_named_ctors "Named mesh constructors".
    /// @{
 
-   /// Deprecated: see @a MakeCartesian3D.
+   /// Deprecated: see @a MakeCartesian3D().
    MFEM_DEPRECATED
    Mesh(int nx, int ny, int nz, Element::Type type, bool generate_edges = false,
         real_t sx = 1.0, real_t sy = 1.0, real_t sz = 1.0,
@@ -1184,7 +1303,7 @@ public:
       Finalize(true); // refine = true
    }
 
-   /// Deprecated: see @a MakeCartesian2D.
+   /// Deprecated: see @a MakeCartesian2D().
    MFEM_DEPRECATED
    Mesh(int nx, int ny, Element::Type type, bool generate_edges = false,
         real_t sx = 1.0, real_t sy = 1.0, bool sfc_ordering = true)
@@ -1194,7 +1313,7 @@ public:
       Finalize(true); // refine = true
    }
 
-   /// Deprecated: see @a MakeCartesian1D.
+   /// Deprecated: see @a MakeCartesian1D().
    MFEM_DEPRECATED
    explicit Mesh(int n, real_t sx = 1.0)
       : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
@@ -1203,7 +1322,7 @@ public:
       // Finalize(); // reminder: not needed
    }
 
-   /// Deprecated: see @a MakeRefined.
+   /// Deprecated: see @a MakeRefined().
    MFEM_DEPRECATED
    Mesh(Mesh *orig_mesh, int ref_factor, int ref_type);
 
@@ -1415,7 +1534,7 @@ public:
    /// Returns the type of boundary element i.
    Element::Type GetBdrElementType(int i) const;
 
-   /// Deprecated in favor of Mesh::GetFaceGeometry
+   /// Deprecated in favor of Mesh::GetFaceGeometry()
    MFEM_DEPRECATED Geometry::Type GetFaceGeometryType(int Face) const
    { return GetFaceGeometry(Face); }
 
@@ -1448,7 +1567,7 @@ public:
       return boundary[i]->GetGeometryType();
    }
 
-   /// Deprecated in favor of Mesh::GetFaceGeometry
+   /// Deprecated in favor of Mesh::GetFaceGeometry()
    MFEM_DEPRECATED Geometry::Type GetFaceBaseGeometry(int i) const
    { return GetFaceGeometry(i); }
 
@@ -1575,9 +1694,9 @@ public:
        @warning This only differs from GetBdrElementAdjacentElement by returning
        the face info with inverted orientation. It does @b not return
        information corresponding to a second adjacent face. This function is
-       deprecated, use Geometry::GetInverseOrientation, Mesh::EncodeFaceInfo,
-       Mesh::DecodeFaceInfoOrientation, and Mesh::DecodeFaceInfoLocalIndex
-       instead.
+       deprecated, use Geometry::GetInverseOrientation(),
+       Mesh::EncodeFaceInfo(), Mesh::DecodeFaceInfoOrientation(), and
+       Mesh::DecodeFaceInfoLocalIndex() instead.
 
        @sa GetBdrElementAdjacentElement() */
    MFEM_DEPRECATED
@@ -2197,7 +2316,7 @@ public:
    void ScaleSubdomains (real_t sf);
    void ScaleElements (real_t sf);
 
-   void Transform(void (*f)(const Vector&, Vector&));
+   void Transform(std::function<void(const Vector&, Vector&)> f);
    void Transform(VectorCoefficient &deformation);
 
    /** @brief This function should be called after the mesh node coordinates

mesh/nurbs.cpp

@@ -2331,7 +2331,7 @@ NURBSExtension::NURBSExtension(NURBSExtension *parent,
    ConnectBoundaries();
 }
 
-NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
+NURBSExtension::NURBSExtension(const Mesh * const mesh_array[], int num_pieces)
 {
    NURBSExtension *parent = mesh_array[0]->NURBSext;
 
@@ -2873,7 +2873,8 @@ void NURBSExtension::GenerateActiveBdrElems()
 }
 
 
-void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
+void NURBSExtension::MergeWeights(const Mesh * const mesh_array[],
+                                  int num_pieces)
 {
    Array<int> lelem_elem;
 
@@ -2899,7 +2900,7 @@ void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
 }
 
 void NURBSExtension::MergeGridFunctions(
-   GridFunction *gf_array[], int num_pieces, GridFunction &merged)
+   const GridFunction * const gf_array[], int num_pieces, GridFunction &merged)
 {
    FiniteElementSpace *gfes = merged.FESpace();
    Array<int> lelem_elem, dofs;
@@ -2907,7 +2908,7 @@ void NURBSExtension::MergeGridFunctions(
 
    for (int i = 0; i < num_pieces; i++)
    {
-      FiniteElementSpace *lfes = gf_array[i]->FESpace();
+      const FiniteElementSpace *lfes = gf_array[i]->FESpace();
       NURBSExtension *lext = lfes->GetMesh()->NURBSext;
 
       lext->GetElementLocalToGlobal(lelem_elem);

mesh/nurbs.hpp

@@ -666,7 +666,7 @@ protected:
 
    /** @brief Set the weights in this object to values from active elements in
        @a num_pieces meshes in @a mesh_array. */
-   void MergeWeights(Mesh *mesh_array[], int num_pieces);
+   void MergeWeights(const Mesh * const mesh_array[], int num_pieces);
 
    /// Set @a patch_to_el.
    void SetPatchToElements();
@@ -696,7 +696,7 @@ public:
                   Mode mode = Mode::H_1);
    /// Construct a NURBSExtension by merging a partitioned NURBS mesh.
 
-   NURBSExtension(Mesh *mesh_array[], int num_pieces);
+   NURBSExtension(const Mesh * const mesh_array[], int num_pieces);
 
    NURBSExtension(const Mesh *patch_topology, const Array<const NURBSPatch*> p);
 
@@ -712,8 +712,8 @@ public:
 
    /** @brief Set the DOFs of @a merged to values from active elements in
        @a num_pieces of Gridfunctions @a gf_array. */
-   void MergeGridFunctions(GridFunction *gf_array[], int num_pieces,
-                           GridFunction &merged);
+   void MergeGridFunctions(const GridFunction * const gf_array[],
+                           int num_pieces, GridFunction &merged);
 
    /// Destroy a NURBSExtension.
    virtual ~NURBSExtension();

miniapps/meshing/mesh-explorer.cpp

@@ -356,15 +356,6 @@ int main (int argc, char *argv[])
             cout << ' ' << mesh->attributes[i];
          }
          cout << endl;
-         cout << "mesh curvature     : ";
-         if (mesh->GetNodalFESpace() != NULL)
-         {
-            cout << mesh->GetNodalFESpace()->FEColl()->Name() << endl;
-         }
-         else
-         {
-            cout << "NONE" << endl;
-         }
       }
       print_char = 0;
       cout << endl;
@@ -380,8 +371,10 @@ int main (int argc, char *argv[])
            "b) View boundary\n"
            "B) View boundary partitioning\n"
            "e) View elements\n"
-           "h) View element sizes, h\n"
-           "k) View element ratios, kappa\n"
+           "h) View element sizes, h, computed at element centers\n"
+           "H) View element sizes, h, computed as a high-order field\n"
+           "k) View element ratios, kappa, computed at element centers\n"
+           "K) View element ratios, kappa, computed as high-order field\n"
            "J) View scaled Jacobian\n"
            "l) Plot a function\n"
            "x) Print sub-element stats\n"
@@ -705,9 +698,8 @@ int main (int argc, char *argv[])
                T->SetIntPoint(&ir.IntPoint(j));
                Geometries.JacToPerfJac(geom, T->Jacobian(), J);
 
-               real_t det_J = J.Det();
-               real_t kappa =
-                  J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1);
+               real_t det_J = J.Weight();
+               real_t kappa = J.CalcConditionNumber();
 
                min_det_J_z = std::min(min_det_J_z, det_J);
                max_det_J_z = std::max(max_det_J_z, det_J);
@@ -837,13 +829,15 @@ int main (int argc, char *argv[])
 
       // These are most of the cases that open a new GLVis window
       if (mk == 'm' || mk == 'b' || mk == 'e' || mk == 'v' || mk == 'h' ||
-          mk == 'k' || mk == 'J' || mk == 'p' || mk == 'B' || mk == 'P')
+          mk == 'k' || mk == 'J' || mk == 'p' || mk == 'B' || mk == 'P' ||
+          mk == 'H' || mk == 'K')
       {
          FiniteElementSpace *bdr_attr_fespace = NULL;
          FiniteElementSpace *attr_fespace =
             new FiniteElementSpace(mesh, attr_fec);
          GridFunction bdr_attr;
          GridFunction attr(attr_fespace);
+         GridFunction ho_func; // high-order function to display, if non-empty
 
          if (mk == 'm')
          {
@@ -926,45 +920,56 @@ int main (int argc, char *argv[])
                  << "- F9/F10     - 3D: cycle through visible elements\n";
          }
 
-         if (mk == 'h')
+         if (mk == 'h' || mk == 'H')
          {
-            DenseMatrix J(dim);
-            real_t h_min, h_max;
-            h_min = infinity();
-            h_max = -h_min;
-            for (int i = 0; i < mesh->GetNE(); i++)
+            int h_type = 0;
+            int ho_func_order = 0;
+            if (mk == 'H')
             {
-               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
-               ElementTransformation *T = mesh->GetElementTransformation(i);
-               T->SetIntPoint(&Geometries.GetCenter(geom));
-               Geometries.JacToPerfJac(geom, T->Jacobian(), J);
-
-               attr(i) = J.Det();
-               if (attr(i) < 0.0)
-               {
-                  attr(i) = -pow(-attr(i), 1.0/real_t(dim));
-               }
-               else
-               {
-                  attr(i) = pow(attr(i), 1.0/real_t(dim));
-               }
-               h_min = min(h_min, attr(i));
-               h_max = max(h_max, attr(i));
+               cout <<
+                    "enter mesh size type:\n"
+                    "0) det(J)^(1/dim) or weight(J)^(1/dim)\n"
+                    "1) h_min = minimal singular value of J\n"
+                    "2) h_max = maximal singular value of J\n"
+                    "--> " << flush;
+               cin >> h_type;
+               cout << "enter FE order for mesh size approximation --> "
+                    << flush;
+               cin >> ho_func_order;
+               ho_func_order = std::max(ho_func_order, 0);
             }
-            cout << "h_min = " << h_min << ", h_max = " << h_max << endl;
+            auto *ho_fec = new L2_FECollection(
+               ho_func_order, mesh->Dimension(), BasisType::GaussLobatto);
+            auto *ho_fes = new FiniteElementSpace(mesh, ho_fec);
+            ho_func.SetSpace(ho_fes);
+            ho_func.MakeOwner(ho_fec);
+            MeshSizeCoefficient h_coeff(h_type);
+            ho_func.ProjectCoefficient(h_coeff);
+            cout << "h range : " << ho_func.Min() << ' ' << ho_func.Max()
+                 << endl;
          }
 
-         if (mk == 'k')
+         if (mk == 'k' || mk == 'K')
          {
-            DenseMatrix J(dim);
-            for (int i = 0; i < mesh->GetNE(); i++)
+            int ho_func_order = 0;
+            if (mk == 'K')
             {
-               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
-               ElementTransformation *T = mesh->GetElementTransformation(i);
-               T->SetIntPoint(&Geometries.GetCenter(geom));
-               Geometries.JacToPerfJac(geom, T->Jacobian(), J);
-               attr(i) = J.CalcSingularvalue(0) / J.CalcSingularvalue(dim-1);
+               cout << "enter FE order for aspect ratio approximation --> "
+                    << flush;
+               cin >> ho_func_order;
+               ho_func_order = std::max(ho_func_order, 0);
             }
+            auto *ho_fec = new L2_FECollection(
+               ho_func_order, mesh->Dimension(), BasisType::GaussLobatto);
+            auto *ho_fes = new FiniteElementSpace(mesh, ho_fec);
+            ho_func.SetSpace(ho_fes);
+            ho_func.MakeOwner(ho_fec);
+            auto kappa_eval = [&](const DenseMatrix &J) -> real_t
+            {
+               return J.CalcConditionNumber();
+            };
+            JacobianFunctionCoefficient kappa_coeff(kappa_eval, true);
+            ho_func.ProjectCoefficient(kappa_coeff);
          }
 
          if (mk == 'J')
@@ -1158,7 +1163,7 @@ int main (int argc, char *argv[])
                      }
                   }
                }
-               attr.Save(sol_sock);
+               sol_sock << (ho_func.Size() ? ho_func : attr);
                sol_sock << "RjlmAb***********";
                if (mk == 'v')
                {
@@ -1209,7 +1214,7 @@ int main (int argc, char *argv[])
                }
                if (mk != 'b' && mk != 'B')
                {
-                  attr.Save(sol_sock);
+                  sol_sock << (ho_func.Size() ? ho_func : attr);
                   sol_sock << "maaA";
                   if (mk == 'v')
                   {
@@ -1237,17 +1242,33 @@ int main (int argc, char *argv[])
          // Project and plot the function 'f'
          int p;
          FiniteElementCollection *fec = NULL;
-         cout << "Enter projection space order: " << flush;
+         cout << "Enter projection space order (<= 0 for DG): " << flush;
          cin >> p;
          if (p >= 1)
          {
+            cout << "Using H1 space of order " << p << " (Gauss-Lobatto nodes)"
+                 << endl;
             fec = new H1_FECollection(p, mesh->Dimension(),
                                       BasisType::GaussLobatto);
          }
          else
          {
-            fec = new DG_FECollection(-p, mesh->Dimension(),
-                                      BasisType::GaussLegendre);
+            cout << "Using DG space of order " << (-p) << endl;
+            cout <<
+                 "Enter DG Basis type:\n"
+                 "0) Gauss-Legendre\n"
+                 "1) Gauss-Lobatto\n"
+                 "2) Positive (Bernstein)\n"
+                 "--> " << flush;
+            int bt;
+            cin >> bt;
+            switch (bt)
+            {
+               case 1:  bt = BasisType::GaussLobatto; break;
+               case 2:  bt = BasisType::Positive; break;
+               default: bt = BasisType::GaussLegendre;
+            }
+            fec = new DG_FECollection(-p, mesh->Dimension(), bt);
          }
          FiniteElementSpace fes(mesh, fec);
          GridFunction level(&fes);