Adding transfer operator support for anisotropic triangles

data/periodic-hexagon-mixed.mesh

@@ -0,0 +1,118 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+2
+
+elements
+16
+1 2 0 1 5 
+1 2 2 3 1
+1 2 4 5 3
+1 2 1 3 5
+2 3 9 6 1 0
+2 3 6 10 2 1
+3 2 10 8 2 
+3 2 8 9 2
+2 3 9 7 3 2
+2 3 7 10 4 3
+3 2 10 6 4
+3 2 6 9 4
+2 3 9 8 5 4
+2 3 8 10 0 5
+3 2 10 7 0
+3 2 7 9 0
+
+boundary
+0
+
+vertices
+11
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+-0.5  -0.2886751345948129
+ 0.0  -0.2886751345948129
+-0.25  0.14433756729740646
+
+ 0.5  -0.2886751345948129
+ 0.25  0.14433756729740646
+ 0.0  -0.2886751345948129
+
+ 0.0   0.5773502691896258
+-0.25  0.14433756729740646
+ 0.25  0.14433756729740646
+
+ 0.0  -0.2886751345948129
+ 0.25  0.14433756729740646
+-0.25  0.14433756729740646
+ 
+-0.5 -0.8660254037844386
+ 0.0 -0.8660254037844386
+-0.5  -0.2886751345948129
+ 0.0  -0.2886751345948129
+
+ 0.0 -0.8660254037844386
+ 0.5 -0.8660254037844386
+ 0.0  -0.2886751345948129
+ 0.5  -0.2886751345948129
+
+ 0.5 -0.8660254037844386
+ 0.75 -0.4330127018922193
+ 0.5  -0.2886751345948129
+
+ 0.75 -0.4330127018922193
+ 1.0   0.0
+ 0.5  -0.2886751345948129
+
+ 1.0   0.0
+ 0.75  0.4330127018922193
+ 0.5  -0.2886751345948129
+ 0.25  0.14433756729740646
+
+ 0.75  0.4330127018922193
+ 0.5   0.8660254037844386
+ 0.25  0.14433756729740646
+ 0.0   0.5773502691896258
+
+ 0.5   0.8660254037844386
+ 0.0   0.8660254037844386
+ 0.0   0.5773502691896258
+  
+ 0.0   0.8660254037844386
+-0.5   0.8660254037844386
+ 0.0   0.5773502691896258
+  
+-0.5   0.8660254037844386
+-0.75  0.4330127018922193
+-0.0   0.5773502691896258
+-0.25  0.14433756729740646
+
+-0.75  0.4330127018922193
+-1.0   0.0
+-0.25  0.14433756729740646
+-0.5  -0.2886751345948129
+
+-1.0   0.0
+-0.75 -0.4330127018922193
+-0.5  -0.2886751345948129
+
+-0.75 -0.4330127018922193
+-0.5  -0.8660254037844386
+-0.5  -0.2886751345948129

fem/derefmat_op.cpp

@@ -166,7 +166,7 @@ DerefineMatrixOp::DerefineMatrixOp(FiniteElementSpace &fespace_, int old_ndofs,
          if (fespace->IsVariableOrder())
          {
             const FiniteElement *fe = fespace->GetFE(emb.parent);
-            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const DenseMatrixStack &pmats = dtrans.point_matrices[geom];
             const int ldof = fe->GetDof();
 
             IsoparametricTransformation isotr;

fem/fespace.cpp

@@ -1642,7 +1642,7 @@ const FaceQuadratureInterpolator
 
 SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
    const int coarse_ndofs, const Table &coarse_elem_dof,
-   const Table *coarse_elem_fos, const DenseTensor localP[]) const
+   const Table *coarse_elem_fos, const DenseMatrixStack localP[]) const
 {
    /// TODO: Implement DofTransformation support
 
@@ -1656,7 +1656,7 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
    SparseMatrix *P;
    if (elem_geoms.Size() == 1)
    {
-      const int coarse_ldof = localP[elem_geoms[0]].SizeJ();
+      const int coarse_ldof = localP[elem_geoms[0]].SizeJ(0);
       P = new SparseMatrix(GetVSize(), coarse_ndofs*vdim, coarse_ldof);
    }
    else
@@ -1674,7 +1674,7 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
       const Embedding &emb = rtrans.embeddings[k];
       const Geometry::Type geom = mesh->GetElementBaseGeometry(k);
       const DenseMatrix &lP = localP[geom](emb.matrix);
-      const int fine_ldof = localP[geom].SizeI();
+      const int fine_ldof = localP[geom].SizeI(emb.matrix);
 
       elem_dof->GetRow(k, dofs);
       coarse_elem_dof.GetRow(emb.parent, coarse_dofs);
@@ -1731,7 +1731,7 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
       isotr.SetIdentityTransformation(geom);
       const int ldof = fe->GetDof();
       lP.SetSize(ldof, ldof);
-      const DenseTensor &pmats = rtrans.point_matrices[geom];
+      const DenseMatrixStack &pmats = rtrans.point_matrices[geom];
       isotr.SetPointMat(pmats(emb.matrix));
       fe->GetLocalInterpolation(isotr, lP);
 
@@ -1766,25 +1766,66 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
 }
 
 void FiniteElementSpace::GetLocalRefinementMatrices(
-   Geometry::Type geom, DenseTensor &localP) const
+   Geometry::Type geom, DenseMatrixStack &localP) const
 {
    const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
+   const int dim = fe->GetDim();
 
    const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
-   const DenseTensor &pmats = rtrans.point_matrices[geom];
+   const DenseMatrixStack &pmats = rtrans.point_matrices[geom];
 
    int nmat = pmats.SizeK();
    int ldof = fe->GetDof();
 
-   IsoparametricTransformation isotr;
-   isotr.SetIdentityTransformation(geom);
-
-   // calculate local interpolation matrices for all refinement types
-   localP.SetSize(ldof, ldof, nmat);
+   // Determine the number of DoFs in each fine element
+   Array<int> fine_ndofs(nmat);
+   Array<Geometry::Type> fine_geom(nmat);
    for (int i = 0; i < nmat; i++)
    {
+      const int num_fine_pts = pmats(i).Width();
+      if (dim == 1)
+      {
+         fine_geom[i] = Geometry::SEGMENT;
+      }
+      else if (dim == 2)
+      {
+         fine_geom[i] =
+            num_fine_pts == 3 ? Geometry::TRIANGLE : Geometry::SQUARE;
+      }
+      else
+      {
+         switch (num_fine_pts)
+         {
+            case 4: fine_geom[i] = Geometry::TETRAHEDRON; break;
+            case 5: fine_geom[i] = Geometry::PYRAMID; break;
+            case 6: fine_geom[i] = Geometry::PRISM; break;
+            case 8: fine_geom[i] = Geometry::CUBE; break;
+         }
+      }
+
+      const FiniteElement *fine_fe =
+         fec->FiniteElementForGeometry(fine_geom[i]);
+      fine_ndofs[i] = fine_fe->GetDof();
+   }
+
+   IsoparametricTransformation isotr;
+
+   // calculate local interpolation matrices for all refinement types
+   localP.SetSize(fine_ndofs, ldof, nmat);
+   for (int i = 0; i < nmat; i++)
+   {
+      isotr.SetIdentityTransformation(fine_geom[i]);
       isotr.SetPointMat(pmats(i));
-      fe->GetLocalInterpolation(isotr, localP(i));
+      if (geom == fine_geom[i])
+      {
+         fe->GetLocalInterpolation(isotr, localP(i));
+      }
+      else
+      {
+         const FiniteElement *fine_fe =
+            fec->FiniteElementForGeometry(fine_geom[i]);
+         fine_fe->GetTransferMatrix(*fe, isotr, localP(i));
+      }
    }
 }
 
@@ -1798,7 +1839,7 @@ SparseMatrix* FiniteElementSpace::RefinementMatrix(int old_ndofs,
    Mesh::GeometryList elem_geoms(*mesh);
    if (!IsVariableOrder())
    {
-      DenseTensor localP[Geometry::NumGeom];
+      DenseMatrixStack localP[Geometry::NumGeom];
       for (int i = 0; i < elem_geoms.Size(); i++)
       {
          GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
@@ -1831,7 +1872,8 @@ FiniteElementSpace::RefinementOperator::RefinementOperator(
    {
       for (int i = 0; i < elem_geoms.Size(); i++)
       {
-         fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
+         fespace->GetLocalRefinementMatrices(elem_geoms[i],
+                                             localP[elem_geoms[i]]);
       }
    }
 
@@ -1943,18 +1985,20 @@ void FiniteElementSpace::RefinementOperator::Mult(const Vector &x,
    {
       const Embedding &emb = trans_ref.embeddings[k];
       const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
+      const Geometry::Type par_geom = (Geometry::Type)emb.par_geom;
+      const Geometry::Type child_geom = (Geometry::Type)emb.child_geom;
       if (fespace->IsVariableOrder())
       {
          const FiniteElement *fe = fespace->GetFE(k);
          isotr.SetIdentityTransformation(geom);
          const int ldof = fe->GetDof();
          eP.SetSize(ldof, ldof);
-         const DenseTensor &pmats = trans_ref.point_matrices[geom];
+         const DenseMatrixStack &pmats = trans_ref.point_matrices[geom];
          isotr.SetPointMat(pmats(emb.matrix));
          fe->GetLocalInterpolation(isotr, eP);
       }
-      const DenseMatrix &lP = (fespace->IsVariableOrder()) ? eP : localP[geom](
-                                 emb.matrix);
+      const DenseMatrix &lP = (fespace->IsVariableOrder()) ?
+                              eP : localP[par_geom](emb.matrix);
 
       subY.SetSize(lP.Height());
 
@@ -2035,7 +2079,7 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
          isotr.SetIdentityTransformation(geom);
          const int ldof = fe->GetDof();
          eP.SetSize(ldof);
-         const DenseTensor &pmats = trans_ref.point_matrices[geom];
+         const DenseMatrixStack &pmats = trans_ref.point_matrices[geom];
          isotr.SetPointMat(pmats(emb.matrix));
          fe->GetLocalInterpolation(isotr, eP);
       }
@@ -2249,7 +2293,7 @@ FiniteElementSpace::DerefinementOperator::DerefinementOperator(
          f_fes->fec->FiniteElementForGeometry(geom);
       const FiniteElement *coarse_fe =
          c_fes->fec->FiniteElementForGeometry(geom);
-      const DenseTensor &pmats = rtrans.point_matrices[geom];
+      const DenseMatrixStack &pmats = rtrans.point_matrices[geom];
 
       lP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), pmats.SizeK());
       lM.SetSize(fine_fe->GetDof(),   fine_fe->GetDof(), pmats.SizeK());
@@ -2367,7 +2411,7 @@ void FiniteElementSpace::GetLocalDerefinementMatrices(Geometry::Type geom,
 
    const CoarseFineTransformations &dtrans =
       mesh->ncmesh->GetDerefinementTransforms();
-   const DenseTensor &pmats = dtrans.point_matrices[geom];
+   const DenseMatrixStack &pmats = dtrans.point_matrices[geom];
 
    const int nmat = pmats.SizeK();
    const int ldof = fe->GetDof();
@@ -2430,7 +2474,7 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
       if (IsVariableOrder())
       {
          fe = GetFE(emb.parent);
-         const DenseTensor &pmats = dtrans.point_matrices[geom];
+         const DenseMatrixStack &pmats = dtrans.point_matrices[geom];
          const int ldof = fe->GetDof();
 
          IsoparametricTransformation isotr;
@@ -2485,27 +2529,59 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
 
 void FiniteElementSpace::GetLocalRefinementMatrices(
    const FiniteElementSpace &coarse_fes, Geometry::Type geom,
-   DenseTensor &localP) const
+   DenseMatrixStack &localP) const
 {
    // Assumptions: see the declaration of the method.
-
-   const FiniteElement *fine_fe = fec->FiniteElementForGeometry(geom);
    const FiniteElement *coarse_fe =
       coarse_fes.fec->FiniteElementForGeometry(geom);
+   const int dim = coarse_fe->GetDim();
 
    const CoarseFineTransformations &rtrans = mesh->GetRefinementTransforms();
-   const DenseTensor &pmats = rtrans.point_matrices[geom];
+   const DenseMatrixStack &pmats = rtrans.point_matrices[geom];
 
    int nmat = pmats.SizeK();
 
    IsoparametricTransformation isotr;
    isotr.SetIdentityTransformation(geom);
 
+   // Determine the number of DoFs in each fine element
+   Array<int> fine_ndofs(nmat);
+   Array<Geometry::Type> fine_geom(nmat);
+   for (int i = 0; i < nmat; i++)
+   {
+      const int num_fine_pts = pmats(i).Width();
+      if (dim == 1)
+      {
+         fine_geom[i] = Geometry::SEGMENT;
+      }
+      else if (dim == 2)
+      {
+         fine_geom[i] =
+            num_fine_pts == 3 ? Geometry::TRIANGLE : Geometry::SQUARE;
+      }
+      else
+      {
+         switch (num_fine_pts)
+         {
+            case 4: fine_geom[i] = Geometry::TETRAHEDRON; break;
+            case 5: fine_geom[i] = Geometry::PYRAMID; break;
+            case 6: fine_geom[i] = Geometry::PRISM; break;
+            case 8: fine_geom[i] = Geometry::CUBE; break;
+         }
+      }
+
+      const FiniteElement *fine_fe =
+         fec->FiniteElementForGeometry(fine_geom[i]);
+      fine_ndofs[i] = fine_fe->GetDof();
+   }
+
    // Calculate the local interpolation matrices for all refinement types
-   localP.SetSize(fine_fe->GetDof(), coarse_fe->GetDof(), nmat);
+   localP.SetSize(fine_ndofs, coarse_fe->GetDof(), nmat);
    for (int i = 0; i < nmat; i++)
    {
       isotr.SetPointMat(pmats(i));
+      const FiniteElement *fine_fe =
+         fec->FiniteElementForGeometry(fine_geom[i]);
       fine_fe->GetTransferMatrix(*coarse_fe, isotr, localP(i));
    }
 }
@@ -4058,7 +4134,7 @@ void FiniteElementSpace::GetTransferOperator(
       {
          Mesh::GeometryList elem_geoms(*mesh);
 
-         DenseTensor localP[Geometry::NumGeom];
+         DenseMatrixStack localP[Geometry::NumGeom];
          for (int i = 0; i < elem_geoms.Size(); i++)
          {
             GetLocalRefinementMatrices(coarse_fes, elem_geoms[i],

fem/fespace.hpp

@@ -532,7 +532,7 @@ protected:
    class RefinementOperator : public Operator
    {
       const FiniteElementSpace* fespace;
-      DenseTensor localP[Geometry::NumGeom];
+      DenseMatrixStack localP[Geometry::NumGeom];
       Table* old_elem_dof; // Owned.
       Table* old_elem_fos; // Owned.
 
@@ -583,7 +583,7 @@ protected:
    SparseMatrix *RefinementMatrix_main(const int coarse_ndofs,
                                        const Table &coarse_elem_dof,
                                        const Table *coarse_elem_fos,
-                                       const DenseTensor localP[]) const;
+                                       const DenseMatrixStack localP[]) const;
 
    /* This method returns the Refinement matrix (i.e., the embedding)
       from a coarse variable-order fes to a fine fes (after a geometric refinement) */
@@ -591,7 +591,7 @@ protected:
                                                const Table &coarse_elem_dof) const;
 
    void GetLocalRefinementMatrices(Geometry::Type geom,
-                                   DenseTensor &localP) const;
+                                   DenseMatrixStack &localP) const;
    void GetLocalDerefinementMatrices(Geometry::Type geom,
                                      DenseTensor &localR) const;
 
@@ -614,7 +614,7 @@ protected:
        NOT variable-order spaces. */
    void GetLocalRefinementMatrices(const FiniteElementSpace &coarse_fes,
                                    Geometry::Type geom,
-                                   DenseTensor &localP) const;
+                                   DenseMatrixStack &localP) const;
 
    /// Help function for constructors + Load().
    void Constructor(Mesh *mesh, NURBSExtension *ext,

fem/pderefmat_op.cpp

@@ -463,7 +463,7 @@ ParDerefineMatrixOp::ParDerefineMatrixOp(ParFiniteElementSpace &fespace_,
          if (fespace->IsVariableOrder())
          {
             const FiniteElement *fe = fespace->GetFE(emb.parent);
-            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const DenseMatrixStack &pmats = dtrans.point_matrices[geom];
             const int ldof = fe->GetDof();
 
             IsoparametricTransformation isotr;

fem/transfer.cpp

@@ -169,7 +169,7 @@ const Operator &InterpolationGridTransfer::ForwardOperator()
    {
       Mesh::GeometryList elem_geoms(*ran_fes.GetMesh());
 
-      DenseTensor localP[Geometry::NumGeom];
+      DenseMatrixStack localP[Geometry::NumGeom];
       for (int i = 0; i < elem_geoms.Size(); i++)
       {
          ran_fes.GetLocalRefinementMatrices(dom_fes, elem_geoms[i],
@@ -399,7 +399,7 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
          });
 
          emb_tr.SetIdentityTransformation(geom);
-         const DenseTensor &pmats = cf_tr.point_matrices[geom];
+         const DenseMatrixStack &pmats = cf_tr.point_matrices[geom];
 
          // Collect the basis functions
          for (int iref = 0; iref < nref; ++iref)
@@ -546,7 +546,7 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
       ElementTransformation *tr_ho = fes_ho.GetElementTransformation(iho);
 
       emb_tr.SetIdentityTransformation(geom);
-      const DenseTensor &pmats = cf_tr.point_matrices[geom];
+      const DenseMatrixStack &pmats = cf_tr.point_matrices[geom];
 
       DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
 
@@ -1617,7 +1617,7 @@ std::unique_ptr<SparseMatrix>>
       ElementTransformation *tr_ho = fes_ho.GetElementTransformation(iho);
 
       emb_tr.SetIdentityTransformation(geom);
-      const DenseTensor& pmats = cf_tr.point_matrices[geom];
+      const DenseMatrixStack& pmats = cf_tr.point_matrices[geom];
 
       int nedof_ho = fe_ho.GetDof();
       int nedof_lor = fe_lor.GetDof();

linalg/densemat.cpp

@@ -4395,6 +4395,16 @@ DenseTensor &DenseTensor::operator=(const DenseTensor &other)
    return *this;
 }
 
+void DenseMatrixStack::InitOffset()
+{
+   off.SetSize(nk);
+   if (nk > 0) { off[0] = 0; }
+   for (int k = 1; k < nk; k++)
+   {
+      off[k] = off[k - 1] + ni[k - 1] * nj[k - 1];
+   }
+}
+
 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL)
 {
    BatchedLinAlg::LUFactor(Mlu, P);

linalg/densemat.hpp

@@ -24,6 +24,7 @@ class DenseMatrix : public Matrix
 {
    friend class DenseTensor;
    friend class DenseMatrixInverse;
+   friend class DenseMatrixStack;
 
 private:
    Memory<real_t> data;
@@ -1312,6 +1313,324 @@ public:
    ~DenseTensor() { tdata.Delete(); }
 };
 
+/// Three index container (stack of matrices)
+class DenseMatrixStack
+{
+private:
+   mutable DenseMatrix Mk;
+   Memory<real_t> tdata;
+   int tsize;
+   Array<int> off, ni, nj;
+   int nk;
+
+   void InitOffset();
+
+public:
+   DenseMatrixStack()
+      : Mk(NULL, 0, 0)
+   {
+      tsize = 0;
+      nk = 0;
+   }
+
+   DenseMatrixStack(const Array<int>& ik, const Array<int>& jk, int k)
+      : Mk(NULL, 0, 0), ni(ik), nj(jk)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize);
+   }
+
+   DenseMatrixStack(real_t *d, const Array<int>& ik, const Array<int>& jk,
+                    int k)
+      : Mk(NULL, 0, 0), ni(ik), nj(jk)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.Wrap(d, tsize, false);
+   }
+
+   DenseMatrixStack(const Array<int>& ik, const Array<int>& jk, int k,
+                    MemoryType mt)
+      : Mk(NULL, 0, 0), ni(ik), nj(jk)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize, mt);
+   }
+
+   /// Copy constructor: deep copy
+   DenseMatrixStack(const DenseMatrixStack &other)
+      : Mk(NULL, 0, 0), off(other.off), ni(other.ni), nj(other.nj), nk(other.nk)
+   {
+      tsize = other.tsize;
+      if (tsize > 0)
+      {
+         tdata.New(tsize, other.tdata.GetMemoryType());
+         tdata.CopyFrom(other.tdata, tsize);
+      }
+   }
+
+   int SizeI(int k) const { return ni[k]; }
+   int SizeJ(int k) const { return nj[k]; }
+   int SizeK() const { return nk; }
+
+   int TotalSize() const { return tsize; }
+
+   void SetSize(int i, int j, int k, MemoryType mt_ = MemoryType::PRESERVE)
+   {
+      const MemoryType mt = mt_ == MemoryType::PRESERVE ?
+                            tdata.GetMemoryType() : mt_;
+      tdata.Delete();
+      Mk.UseExternalData(NULL, i, j);
+      ni.SetSize(k);
+      nj.SetSize(k);
+      if (k > 0) { ni = i; }
+      if (k > 0) { nj = j; }
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize, mt);
+   }
+
+   void SetSize(int i, const Array<int>& jk, int k,
+                MemoryType mt_ = MemoryType::PRESERVE)
+   {
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      const MemoryType mt = mt_ == MemoryType::PRESERVE ?
+                            tdata.GetMemoryType() : mt_;
+      tdata.Delete();
+      Mk.UseExternalData(NULL, 0, 0);
+      ni.SetSize(k);
+      if (k > 0) { ni = i; }
+      nj = jk;
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize, mt);
+   }
+
+   void SetSize(const Array<int>& ik, int j, int k,
+                MemoryType mt_ = MemoryType::PRESERVE)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      const MemoryType mt = mt_ == MemoryType::PRESERVE ?
+                            tdata.GetMemoryType() : mt_;
+      tdata.Delete();
+      Mk.UseExternalData(NULL, 0, 0);
+      ni = ik;
+      nj.SetSize(k);
+      if (k > 0) { nj = j; }
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize, mt);
+   }
+
+   void SetSize(const Array<int>& ik, const Array<int>& jk, int k,
+                MemoryType mt_ = MemoryType::PRESERVE)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      const MemoryType mt = mt_ == MemoryType::PRESERVE ?
+                            tdata.GetMemoryType() : mt_;
+      tdata.Delete();
+      Mk.UseExternalData(NULL, 0, 0);
+      ni = ik;
+      nj = jk;
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.New(tsize, mt);
+   }
+
+   void UseExternalData(real_t *ext_data, int i, int j, int k)
+   {
+      tdata.Delete();
+      Mk.UseExternalData(NULL, i, j);
+      ni.SetSize(k);
+      nj.SetSize(k);
+      ni = i;
+      nj = j;
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.Wrap(ext_data, tsize, false);
+   }
+
+   void UseExternalData(real_t *ext_data,
+                        const Array<int>& ik, const Array<int>& jk, int k)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      tdata.Delete();
+      Mk.UseExternalData(NULL, 0, 0);
+      ni = ik;
+      nj = jk;
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      tdata.Wrap(ext_data, tsize, false);
+   }
+
+   /// @brief Reset the DenseTensor to use the given external Memory @a mem and
+   /// dimensions @a i, @a j, and @a k.
+   ///
+   /// If @a own_mem is false, the DenseTensor will not own any of the pointers
+   /// of @a mem.
+   ///
+   /// Note that when @a own_mem is true, the @a mem object can be destroyed
+   /// immediately by the caller but `mem.Delete()` should NOT be called since
+   /// the DenseTensor object takes ownership of all pointers owned by @a mem.
+   void NewMemoryAndSize(const Memory<real_t> &mem,
+                         const Array<int>& ik, const Array<int>& jk, int k,
+                         bool own_mem)
+   {
+      MFEM_ASSERT(ik.Size() == k, "incompatible Array size");
+      MFEM_ASSERT(jk.Size() == k, "incompatible Array size");
+      tdata.Delete();
+      Mk.UseExternalData(NULL, 0, 0);
+      ni = ik;
+      nj = jk;
+      nk = k;
+      InitOffset();
+      tsize = (k > 0) ? off[nk - 1] + ni[nk - 1] * nj[nk - 1] : 0;
+      if (own_mem)
+      {
+         tdata = mem;
+      }
+      else
+      {
+         tdata.MakeAlias(mem, 0, tsize);
+      }
+   }
+
+   /// Sets the tensor elements equal to constant c
+   DenseMatrixStack &operator=(real_t c);
+
+   /// Copy assignment operator (performs a deep copy)
+   DenseMatrixStack &operator=(const DenseTensor &other);
+
+   DenseMatrix &operator()(int k)
+   {
+      return operator()(k, Mk);
+   }
+   const DenseMatrix &operator()(int k) const
+   {
+      return operator()(k, Mk);
+   }
+   DenseMatrix &operator()(int k, DenseMatrix& buff)
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      buff.UseExternalData(nullptr, SizeI(k), SizeJ(k));
+      buff.data = Memory<real_t>(GetData(k), SizeI(k)*SizeJ(k), false);
+      return buff;
+   }
+   const DenseMatrix &operator()(int k, DenseMatrix& buff) const
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      buff.UseExternalData(nullptr, SizeI(k), SizeJ(k));
+      buff.data = Memory<real_t>(const_cast<real_t*>(GetData(k)),
+                                 SizeI(k)*SizeJ(k),
+                                 false);
+      return buff;
+   }
+
+   real_t &operator()(int i, int j, int k)
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ(k));
+      MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI(k));
+      return tdata[off[k] + i + SizeI(k) * j];
+   }
+
+   const real_t &operator()(int i, int j, int k) const
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ(k));
+      MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI(k));
+      return tdata[off[k] + i + SizeI(k) * j];
+   }
+
+   real_t *GetData(int k)
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      return tdata + off[k];
+   }
+
+   const real_t *GetData(int k) const
+   {
+      MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
+      return tdata + off[k];
+   }
+
+   real_t *Data() { return tdata; }
+
+   const real_t *Data() const { return tdata; }
+
+   Memory<real_t> &GetMemory() { return tdata; }
+   const Memory<real_t> &GetMemory() const { return tdata; }
+
+   /** Matrix-vector product from unassembled element matrices, assuming both
+       'x' and 'y' use the same elem_dof table. */
+   void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;
+
+   void Clear()
+   {
+      tdata.Delete();
+      tsize = 0;
+      Mk.UseExternalData(NULL, 0, 0);
+      off.LoseData();
+      ni.LoseData();
+      nj.LoseData();
+      nk = 0;
+   }
+
+   std::size_t MemoryUsage() const { return nk*Mk.MemoryUsage(); }
+
+   /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
+   const real_t *Read(bool on_dev = true) const
+   { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
+
+   /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
+   const real_t *HostRead() const
+   { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, false); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
+   real_t *Write(bool on_dev = true)
+   { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
+
+   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
+   real_t *HostWrite()
+   { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, false); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
+   real_t *ReadWrite(bool on_dev = true)
+   { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
+
+   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
+   real_t *HostReadWrite()
+   { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, false); }
+
+   void Swap(DenseMatrixStack &t)
+   {
+      mfem::Swap(tdata, t.tdata);
+      mfem::Swap(nk, t.nk);
+      Mk.Swap(t.Mk);
+   }
+
+   ~DenseMatrixStack() { tdata.Delete(); }
+};
+
 /** @brief Compute the LU factorization of a batch of matrices. Calls
     BatchedLinAlg::LUFactor.
 

makefile

@@ -793,7 +793,6 @@ status info:
 	$(info MFEM_MPI_NP            = $(MFEM_MPI_NP))
 	@true
 
-ASTYLE_BIN = astyle
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
 FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)

mesh/element.hpp

@@ -35,6 +35,12 @@ protected:
    /// Element's type from the Finite Element's perspective
    Geometry::Type base_geom;
 
+   /// This is the default data previously used for all element
+   /// types. It is defined here temporarily but each element type
+   /// will eventually define its own data.
+   // static constexpr int def_ref_type_num_children[8] =
+   // { 0, 2, 2, 4, 2, 4, 4, 8 };
+
 public:
 
    /// Constants for the classes derived from Element.
@@ -110,6 +116,9 @@ public:
 
    /// Destroys element.
    virtual ~Element() { }
+
+   // static unsigned int GetNumChildren(int geom_type, int ref_type)
+   //  { return def_ref_type_num_children[ref_type]; }
 };
 
 }

mesh/mesh.cpp

@@ -11775,7 +11775,7 @@ const CoarseFineTransformations &Mesh::GetRefinementTransforms() const
             CoarseFineTr.embeddings[j].matrix = index;
          }
 
-         DenseTensor &pmats = CoarseFineTr.point_matrices[geom];
+         DenseMatrixStack &pmats = CoarseFineTr.point_matrices[geom];
          pmats.SetSize(Dim, Dim+1, static_cast<int>((mat_no.size())));
 
          // calculate the point matrices used

mesh/mesh.hpp

@@ -1597,6 +1597,24 @@ public:
       { mesh.GetGeometries(dim, *this); }
    };
 
+   class PotentialGeometryList : public Array<Geometry::Type>
+   {
+   protected:
+      Geometry::Type geom_buf[Geometry::NumGeom];
+   public:
+      /// Construct a GeometryList of all element geometries in @a mesh.
+      PotentialGeometryList(const Mesh &mesh)
+         : Array<Geometry::Type>(geom_buf, Geometry::NumGeom)
+      {
+         SetSize(0);
+         const int dim = mesh.Dimension();
+         for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
+         {
+            this->Append(Geometry::Type(g));
+         }
+      }
+   };
+
    /// @name Access connectivity for individual mesh entites
    /// @{
 

mesh/ncmesh.cpp

@@ -124,6 +124,16 @@ void NCMesh::GeomInfo::InitGeom(Geometry::Type geom)
    initialized = true;
 }
 
+void NCMesh::InitGeomInfoByDim(int dim, GeomInfo gi[])
+{
+   MFEM_ASSERT(dim > 0 && dim <= Geometry::MaxDim, "unsupported dimension");
+
+   for (int t = Geometry::DimStart[dim]; t < Geometry::DimStart[dim+1]; t++)
+   {
+      gi[t].InitGeom((Geometry::Type)t);
+   }
+}
+
 NCMesh::NCMesh(const Mesh *mesh)
    : shadow(1024, 2048)
 {
@@ -133,6 +143,8 @@ NCMesh::NCMesh(const Mesh *mesh)
    Iso = true;
    Legacy = false;
 
+   InitGeomInfoByDim(Dim, GI);
+
    // create the NCMesh::Element struct for each Mesh element
    for (int i = 0; i < mesh->GetNE(); i++)
    {
@@ -140,13 +152,6 @@ NCMesh::NCMesh(const Mesh *mesh)
 
       Geometry::Type geom = elem->GetGeometryType();
       CheckSupportedGeom(geom);
-      GI[geom].InitGeom(geom);
-
-      // if we have pyramids we will need tets after refinement
-      if (geom == Geometry::PYRAMID)
-      {
-         GI[Geometry::TETRAHEDRON].InitGeom(Geometry::TETRAHEDRON);
-      }
 
       // create NCMesh::Element for this mfem::Element
       int root_id = AddElement(geom, elem->GetAttribute());
@@ -533,7 +538,9 @@ void Refinement::SetScale(const ScaledType &r)
          s[2] = r.second;
          break;
       default:
-         MFEM_ABORT("Unsupported refinement type.");
+         s[0] = r.second;
+         s[1] = r.second;
+         s[2] = r.second;
    }
 }
 
@@ -549,12 +556,14 @@ Refinement::Refinement(int index)
 {
    for (int i=0; i<3; ++i) { s[i] = 0.0; }
    // Default case is XYZ type with scale 0.5.
+   t = XYZ;
    SetScale(ScaledType(XYZ, 0.5));
 }
 
 Refinement::Refinement(int index, const std::initializer_list<ScaledType> &refs)
    : index(index)
 {
+   t = XYZ;
    for (int i=0; i<3; ++i) { s[i] = 0.0; }
    if (refs.size() == 0)
    {
@@ -573,28 +582,27 @@ Refinement::Refinement(int index, const std::initializer_list<ScaledType> &refs)
 Refinement::Refinement(int index, char type, real_t scale)
    : index(index)
 {
+   t = XYZ;
    for (int i=0; i<3; ++i) { s[i] = 0.0; }
    SetScale(ScaledType(type, scale));
 }
 
 char Refinement::GetType() const
 {
-   char t{0}; // Set the X, Y or Z bit
-   for (int i = 0; i < 3; ++i)
-      if (s[i] > real_t{0})
-         t |= (1 << i);
    return t;
 }
 
 void Refinement::Set(int element, char type, real_t scale)
 {
    index = element;
+   t = type;
    for (int i=0; i<3; ++i) { s[i] = 0.0; }
    SetScale(ScaledType(type, scale));
 }
 
 void Refinement::SetType(char type, real_t scale)
 {
+   t = type;
    for (int i=0; i<3; ++i) { s[i] = 0.0; }
    SetScale(ScaledType(type, scale));
 }
@@ -612,6 +620,10 @@ NCMesh::Element::Element(Geometry::Type geom, int attr)
    // keep the code as simple as possible.
 }
 
+unsigned int NCMesh::Element::GetNumChildren(int iref_type) const
+{ return ref_type_num_children[(int)geom][iref_type]; }
+
+
 int NCMesh::NewHexahedron(int n0, int n1, int n2, int n3,
                           int n4, int n5, int n6, int n7,
                           int attr,
@@ -1883,17 +1895,87 @@ void NCMesh::RefineElement(const Refinement & ref)
    }
    else if (el.Geom() == Geometry::TRIANGLE)
    {
-      ref_type = Refinement::XY; // for consistency
+      if ((int)ref_type == 1)
+      {
+         // split first edge
+         int mid01 = nodes.GetId(no[0], no[1]);
 
-      // isotropic split - the only ref_type available for triangles
-      int mid01 = nodes.GetId(no[0], no[1]);
-      int mid12 = nodes.GetId(no[1], no[2]);
-      int mid20 = nodes.GetId(no[2], no[0]);
+         child[0] = NewTriangle(no[0], mid01, no[2], attr, fa[0], -1, fa[2]);
+         child[1] = NewTriangle(mid01, no[1], no[2], attr, fa[0], fa[1], -1);
+      }
+      else if ((int)ref_type == 2)
+      {
+         // split second edge
+         int mid12 = nodes.GetId(no[1], no[2]);
 
-      child[0] = NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
-      child[1] = NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
-      child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
-      child[3] = NewTriangle(mid12, mid20, mid01, attr, -1, -1, -1);
+         child[0] = NewTriangle(no[0], no[1], mid12, attr, fa[0], fa[1], -1);
+         child[1] = NewTriangle(no[0], mid12, no[2], attr, -1, fa[1], fa[2]);
+      }
+      else if ((int)ref_type == 3)
+      {
+         // split first and second edges
+         int mid01 = nodes.GetId(no[0], no[1]);
+         int mid12 = nodes.GetId(no[1], no[2]);
+
+         child[0] = NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
+         child[1] = NewQuadrilateral(no[0], mid01, mid12, no[2],
+                                     attr, fa[0], -1, fa[1], fa[2]);
+      }
+      else if ((int)ref_type == 4)
+      {
+         // split third edge
+         int mid20 = nodes.GetId(no[2], no[0]);
+
+         child[0] = NewTriangle(no[0], no[1], mid20, attr, fa[0], -1, fa[2]);
+         child[1] = NewTriangle(mid20, no[1], no[2], attr, -1, fa[1], fa[2]);
+      }
+      else if ((int)ref_type == 5)
+      {
+         // split first and third edges
+         int mid01 = nodes.GetId(no[0], no[1]);
+         int mid20 = nodes.GetId(no[2], no[0]);
+
+         child[0] = NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
+         child[1] = NewQuadrilateral(mid01, no[1], no[2], mid20,
+                                     attr, fa[0], fa[1], fa[2], -1);
+      }
+      else if ((int)ref_type == 6)
+      {
+         // split second and third edges
+         int mid12 = nodes.GetId(no[1], no[2]);
+         int mid20 = nodes.GetId(no[2], no[0]);
+
+         child[0] = NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
+         child[1] = NewQuadrilateral(no[0], no[1], mid12, mid20,
+                                     attr, fa[0], fa[1], -1, fa[2]);
+      }
+      else if ((int)ref_type == 7)
+      {
+         // isotropic split into four triangles
+         int mid01 = nodes.GetId(no[0], no[1]);
+         int mid12 = nodes.GetId(no[1], no[2]);
+         int mid20 = nodes.GetId(no[2], no[0]);
+
+         child[0] = NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
+         child[1] = NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
+         child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
+         child[3] = NewTriangle(mid12, mid20, mid01, attr, -1, -1, -1);
+      }
+      else if ((int)ref_type == 8)
+      {
+         // isotropic split into three triangles
+
+         // *** the following fails to obtain a useable node index (obviously).
+         //
+         // If we want to allow refinements which add nodes to the
+         // element without adding edge nodes we will need to find a
+         // way to add these face-centered nodes alongside the new edge nodes.
+         int midel = faces.GetId(no[0], no[1], no[2]);
+
+         child[0] = NewTriangle(no[0], no[1], midel, attr, fa[0], -1, -1);
+         child[1] = NewTriangle(no[1], no[2], midel, attr, fa[1], -1, -1);
+         child[2] = NewTriangle(no[2], no[0], midel, attr, fa[2], -1, -1);
+      }
    }
    else if (el.Geom() == Geometry::SEGMENT)
    {
@@ -4622,6 +4704,64 @@ const NCMesh::PointMatrix& NCMesh::GetGeomIdentity(Geometry::Type geom)
    }
 }
 
+int NCMesh::GetPointMatrixNumPoints(Geometry::Type geom,
+                                    const char* ref_path) const
+{
+   while (*ref_path)
+   {
+      int ref_type = *ref_path++;
+      int child = *ref_path++;
+
+      // TODO: do this with the new child transform tables
+
+      if (geom == Geometry::CUBE)
+      {
+         return 8;
+      }
+      else if (geom == Geometry::PRISM)
+      {
+         return 6;
+      }
+      else if (geom == Geometry::PYRAMID)
+      {
+         if (child <= 5) // Pyramid
+         {
+            return 5;
+         }
+         else // Tet
+         {
+            return 4;
+         }
+      }
+      else if (geom == Geometry::TETRAHEDRON)
+      {
+         return 4;
+      }
+      else if (geom == Geometry::SQUARE)
+      {
+         return 4;
+      }
+      else if (geom == Geometry::TRIANGLE)
+      {
+         if ((ref_type == 3|| ref_type == 5 || ref_type == 6) && child == 1)
+         {
+            return 4;
+         }
+         else
+         {
+            return 3;
+         }
+      }
+      else if (geom == Geometry::SEGMENT)
+      {
+         return 2;
+      }
+   }
+
+   // ref_path was an empty string so return the default for this Geometry type
+   return GetPointMatrixNumPoints(geom, "\1");
+}
+
 void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
                             DenseMatrix& matrix) const
 {
@@ -5060,23 +5200,107 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
       }
       else if (geom == Geometry::TRIANGLE)
       {
-         Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
+         if (ref_type == 1) // split first edge
+         {
+            Point mid01(pm(0), pm(1));
 
-         if (child == 0)
-         {
-            pm = PointMatrix(pm(0), mid01, mid20);
+            if (child == 0)
+            {
+               pm = PointMatrix(pm(0), mid01, pm(2));
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(mid01, pm(1), pm(2));
+            }
          }
-         else if (child == 1)
+         else if (ref_type == 2) // split second edge
          {
-            pm = PointMatrix(mid01, pm(1), mid12);
+            Point mid12(pm(1), pm(2));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(pm(0), pm(1), mid12);
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(pm(0), mid12, pm(2));
+            }
          }
-         else if (child == 2)
+         else if (ref_type == 3) // split first and second edges
          {
-            pm = PointMatrix(mid20, mid12, pm(2));
+            Point mid01(pm(0), pm(1));
+            Point mid12(pm(1), pm(2));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(mid01, pm(1), mid12);
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(pm(0), mid01, mid12, pm(2));
+            }
          }
-         else if (child == 3)
+         else if (ref_type == 4) // split third edge
          {
-            pm = PointMatrix(mid12, mid20, mid01);
+            Point mid20(pm(2), pm(0));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(pm(0), pm(1), mid20);
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(mid20, pm(1), pm(2));
+            }
+         }
+         else if (ref_type == 5) // split first and third edges
+         {
+            Point mid01(pm(0), pm(1));
+            Point mid20(pm(2), pm(0));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(pm(0), mid01, mid20);
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(mid01, pm(1), pm(2), mid20);
+            }
+         }
+         else if (ref_type == 6) // split second and third edges
+         {
+            Point mid12(pm(1), pm(2));
+            Point mid20(pm(2), pm(0));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(mid20, mid12, pm(2));
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(pm(0), pm(1), mid12, mid20);
+            }
+         }
+         else if (ref_type == 7) // iso split
+         {
+            Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
+
+            if (child == 0)
+            {
+               pm = PointMatrix(pm(0), mid01, mid20);
+            }
+            else if (child == 1)
+            {
+               pm = PointMatrix(mid01, pm(1), mid12);
+            }
+            else if (child == 2)
+            {
+               pm = PointMatrix(mid20, mid12, pm(2));
+            }
+            else if (child == 3)
+            {
+               pm = PointMatrix(mid12, mid20, mid01);
+            }
          }
       }
       else if (geom == Geometry::SEGMENT)
@@ -5095,6 +5319,7 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
    }
 
    // write the points to the matrix
+   matrix.SetSize(pm(0).dim, pm.np);
    for (int i = 0; i < pm.np; i++)
    {
       for (int j = 0; j < pm(i).dim; j++)
@@ -5171,6 +5396,7 @@ void NCMesh::TraverseRefinements(int elem, int coarse_index,
                                  std::string &ref_path, RefPathMap &map) const
 {
    const Element &el = elements[elem];
+   const Element &coarse_el = elements[coarse_index];
    if (!el.ref_type)
    {
       int &matrix = map[ref_path];
@@ -5180,6 +5406,8 @@ void NCMesh::TraverseRefinements(int elem, int coarse_index,
       emb.parent = coarse_index;
       emb.matrix = matrix - 1;
       emb.geom = el.Geom();
+      emb.par_geom = coarse_el.Geom();
+      emb.child_geom = el.Geom();
       emb.ghost = IsGhost(el);
    }
    else
@@ -5234,13 +5462,19 @@ const CoarseFineTransformations& NCMesh::GetRefinementTransforms() const
          if (used_geoms & (1 << g))
          {
             Geometry::Type geom = Geometry::Type(g);
-            const PointMatrix &identity = GetGeomIdentity(geom);
 
+            Array<int> num_pts(path_map[g].size());
+
+            RefPathMap::iterator it;
+            int pos = 0;
+            for (it = path_map[g].begin(); it != path_map[g].end(); ++it, ++pos)
+            {
+               num_pts[pos] = GetPointMatrixNumPoints(geom, it->first.c_str());
+            }
             transforms.point_matrices[g]
-            .SetSize(Dim, identity.np, static_cast<int>(path_map[g].size()));
+            .SetSize(Dim, num_pts, static_cast<int>(path_map[g].size()));
 
             // calculate the point matrices
-            RefPathMap::iterator it;
             for (it = path_map[g].begin(); it != path_map[g].end(); ++it)
             {
                GetPointMatrix(geom, it->first.c_str(),
@@ -6551,6 +6785,8 @@ NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
       input >> ident;
    }
 
+   InitGeomInfoByDim(Dim, GI);
+
    // load elements
    MFEM_VERIFY(ident == "elements", "Invalid mesh file: " << ident);
    input >> count;
@@ -6569,7 +6805,6 @@ NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
       if (geom >= 0)
       {
          CheckSupportedGeom(type);
-         GI[geom].InitGeom(type);
 
          input >> ref_type;
          MFEM_VERIFY(ref_type >= 0 && ref_type < 8, "");
@@ -6577,7 +6812,7 @@ NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
 
          if (ref_type) // refined element
          {
-            for (int j = 0; j < ref_type_num_children[ref_type]; j++)
+            for (unsigned int j = 0; j < el.GetNumChildren(ref_type); j++)
             {
                input >> el.child[j];
             }
@@ -6727,7 +6962,7 @@ void NCMesh::LoadCoarseElements(std::istream &input)
       if (Dim == 3 && ref_type != 7) { iso = false; }
 
       // load child IDs and make parent-child links
-      int nch = ref_type_num_children[ref_type];
+      int nch = el.GetNumChildren(ref_type);
       for (int i = 0, id; i < nch; i++)
       {
          input >> id;
@@ -6793,6 +7028,8 @@ void NCMesh::LoadLegacyFormat(std::istream &input, int &curved, int &is_nc)
    MFEM_VERIFY(ident == "dimension", "invalid mesh file");
    input >> Dim;
 
+   InitGeomInfoByDim(Dim, GI);
+
    // load elements
    skip_comment_lines(input, '#');
    input >> ident;
@@ -6805,7 +7042,6 @@ void NCMesh::LoadLegacyFormat(std::istream &input, int &curved, int &is_nc)
 
       Geometry::Type type = Geometry::Type(geom);
       CheckSupportedGeom(type);
-      GI[geom].InitGeom(type);
 
       int eid = AddElement(type, attr);
       MFEM_ASSERT(eid == i, "");

mesh/ncmesh.hpp

@@ -40,8 +40,9 @@ namespace mfem
 struct Refinement
 {
    int index; ///< Mesh element number
-   enum : char { X = 1, Y = 2, Z = 4, XY = 3, XZ = 5, YZ = 6, XYZ = 7 };
+   enum : char { X = 1, Y = 2, Z = 4, XY = 3, XZ = 5, YZ = 6, XYZ = 7, T8 = 8 };
    using ScaledType = std::pair<char, real_t>;
+   char t;
    real_t s[3];  /// Refinement scale in each dimension
    Refinement() = default;
    /// Refinement type XYZ, with scale 0.5.
@@ -77,6 +78,9 @@ struct Embedding
    unsigned geom : 4;
    unsigned matrix : 27;
 
+   unsigned par_geom : 4;
+   unsigned child_geom : 4;
+
    /// For internal use: 0 if regular fine element, 1 if parallel ghost element.
    unsigned ghost : 1;
 
@@ -93,7 +97,7 @@ struct CoarseFineTransformations
 
    /** A "dictionary" of matrices for IsoparametricTransformation. Use
        Embedding::{geom,matrix} to access a fine element point matrix. */
-   DenseTensor point_matrices[Geometry::NumGeom];
+   DenseMatrixStack point_matrices[Geometry::NumGeom];
 
    /** Invert the 'embeddings' array: create a Table with coarse elements as
        rows and fine elements as columns. If 'want_ghosts' is false, parallel
@@ -676,6 +680,8 @@ protected: // implementation
       Geometry::Type Geom() const { return Geometry::Type(geom); }
       bool IsLeaf() const { return !ref_type && (parent != -2); }
       int GetAttribute() const { return attribute; }
+
+      unsigned int GetNumChildren(int iref_type) const;
    };
 
 
@@ -1247,6 +1253,9 @@ protected:
 
    static const PointMatrix& GetGeomIdentity(Geometry::Type geom);
 
+   int GetPointMatrixNumPoints(Geometry::Type geom,
+                               const char* ref_path) const;
+
    void GetPointMatrix(Geometry::Type geom, const char* ref_path,
                        DenseMatrix& matrix) const;
 
@@ -1409,6 +1418,8 @@ protected:
 
    static GeomInfo GI[Geometry::NumGeom];
 
+   static void InitGeomInfoByDim(int dim, GeomInfo gi[]);
+
    /// This is used for a NURBS mesh with this NCMesh as its patch topology.
    VertexToKnotSpan vertex_to_knotspan;
 

mesh/ncmesh_tables.hpp

@@ -15,7 +15,26 @@
 namespace mfem
 {
 
-static constexpr int ref_type_num_children[8] = { 0, 2, 2, 4, 2, 4, 4, 8 };
+static constexpr int def_ref_type_num_children[] = { 0, 2, 2, 4, 2, 4, 4, 8 };
+static constexpr int tri_ref_type_num_children[] = { 0,
+                                                     2, 2, 2,
+                                                     2, 2, 2,
+                                                     4,  3,
+                                                     4, 4, 4,
+                                                     4, 4, 4,
+                                                     3
+                                                   };
+static const int* const ref_type_num_children[Geometry::NUM_GEOMETRIES] =
+{
+   NULL,
+   def_ref_type_num_children,
+   tri_ref_type_num_children,
+   def_ref_type_num_children,
+   def_ref_type_num_children,
+   def_ref_type_num_children,
+   def_ref_type_num_children,
+   def_ref_type_num_children
+};
 
 // derefinement tables
 // The first n numbers in each line are the refined elements that contain
@@ -498,4 +517,4 @@ static const RefTrf** geom_child[7] =
 
 } // namespace mfem
 
-#endif // MFEM_NCMESH_TABLES
+#endif // MFEM_NCMESH_TABLES

tests/unit/fem/test_derefine.cpp

@@ -144,6 +144,10 @@ int PolyCoeff::order_ = -1;
 
 void test_derefine_L2_element(int order, Element::Type el_type, int basis_type)
 {
+   CAPTURE(order);
+   CAPTURE(el_type);
+   CAPTURE(basis_type);
+
    Mesh mesh;
    if (dimension == 1)
    {
@@ -194,6 +198,9 @@ void test_derefine_L2_element(int order, Element::Type el_type, int basis_type)
 
    double mass_coarse = integrate(&x);
 
+   CAPTURE(mass_fine);
+   CAPTURE(mass_coarse);
+
    // conservation check
    REQUIRE( fabs(mass_fine-mass_coarse) < 1.e-12 );