Initial commit, CalcJShape implemented for all RT_Triangle, Mult is modified for DenseTensors and new (not working yet) integrators added

CHANGELOG

@@ -13,13 +13,6 @@ Version 3.4.1 (development)
 - Added support for reading linear and quadratic 2D quadrilateral and triangular
   Cubit meshes.
 
-- The tetrahedral mesh refinement algorithm in serial and in parallel now
-  follows precisely the paper:
-     D. Arnold, A. Mukherjee, and L. Pouly, "Locally Adapted Tetrahedral Meshes
-     Using Bisection", SIAM J. Sci. Comput., 22(2), 431–448.
-  This guarantees that the shape regularity of the elements will be preserved
-  under refinement.
-
 
 Version 3.4, released on May 29, 2018
 =====================================

config/test.mk

@@ -38,7 +38,7 @@ export TIME='%es %MkB %x'; \
 set -- $$($(1) $(SHELL) -c "$(2)" 2>&1); while [ "$$#" -gt 3 ]; do shift; done
 endef
 define TIMECMD.NOTGNU
-set -- $$($(1) -l $(SHELL) -c "{ $(2); } > /dev/null 2>&1" 2>&1; echo $$?); \
+set -- $$($(1) -l $(SHELL) -c "$(2)" 2>&1; echo $$?); \
 set -- "$$1"s "$$(($$7/1024))"kB "$${60}"
 endef
 define TIMECMD.BASH

data/inline-tri.mesh

@@ -1,7 +1,7 @@
 MFEM INLINE mesh v1.0
 
 type = tri
-nx = 4
-ny = 4
+nx = 1
+ny = 1
 sx = 1.0
 sy = 1.0

examples/ex13p.cpp

@@ -110,7 +110,6 @@ int main(int argc, char *argv[])
    {
       pmesh->UniformRefinement();
    }
-   pmesh->ReorientTetMesh();
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Nedelec finite elements of the specified order.

examples/ex4.cpp

@@ -94,8 +94,8 @@ int main(int argc, char *argv[])
    //    largest number that gives a final mesh with no more than 25,000
    //    elements.
    {
-      int ref_levels =
-         (int)floor(log(25000./mesh->GetNE())/log(2.)/dim);
+      int ref_levels = 2;
+         //(int)floor(log(25000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();
@@ -136,17 +136,33 @@ int main(int argc, char *argv[])
    //    when eliminating the non-homogeneous boundary condition to modify the
    //    r.h.s. vector b.
    GridFunction x(fespace);
+   GridFunction xp(fespace);
    VectorFunctionCoefficient F(sdim, F_exact);
-   x.ProjectCoefficient(F);
+   xp.ProjectCoefficient(F);
+   VectorGridFunctionCoefficient Xp(&xp);
+//   xp.Print();
+   
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << xp
+               << "window_title 'Projected'\n" << flush;
+   }
 
    // 8. Set up the bilinear form corresponding to the H(div) diffusion operator
    //    grad alpha div + beta I, by adding the div-div and the mass domain
    //    integrators.
-   Coefficient *alpha = new ConstantCoefficient(1.0);
-   Coefficient *beta  = new ConstantCoefficient(1.0);
+   Coefficient *alpha = new ConstantCoefficient(0.0);
+   Coefficient *beta  = new ConstantCoefficient(0.0);
    BilinearForm *a = new BilinearForm(fespace);
-   a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
-   a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
+//   a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
+//   a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
+   a->AddDomainIntegrator(new VectorFEDiffusionIntegrator);
+   a->AddBdrFaceIntegrator(new VectorFEDGDiffusionIntegrator(-1., 10000.));
+   a->AddInteriorFaceIntegrator(new VectorFEDGDiffusionIntegrator(-1., 10000.));
 
    // 9. Assemble the bilinear form and the corresponding linear system,
    //    applying any necessary transformations such as: eliminating boundary
@@ -166,31 +182,43 @@ int main(int argc, char *argv[])
                              ess_tdof_list);
    }
    a->Assemble();
+   a->Finalize();
 
-   SparseMatrix A;
-   Vector B, X;
-   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+   SparseMatrix A = a->SpMat();
+//   b->Print();
+//   printf("lala\n");
+   Vector B(*b), X(*b);
+//   B.Print();
+//   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+   
+   Vector A_times_x = xp;
+   A.Mult(xp, A_times_x);
+   A_times_x.Add(-1.0, *b);
+   cout << "norm of Ax-b: " << A_times_x.Norml2() << endl << flush;
 
-   cout << "Size of linear system: " << A.Height() << endl;
+//   cout << "Size of linear system: " << A.Height() << endl;
 
 #ifndef MFEM_USE_SUITESPARSE
    // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
    //     solve the system A X = B with PCG.
    GSSmoother M(A);
-   PCG(A, M, B, X, 1, 10000, 1e-20, 0.0);
+   int maxit = 1000;
+   PCG(A, M, B, X, 1, maxit, 1e-20, 0.0);
 #else
    // 10. If compiled with SuiteSparse support, use UMFPACK to solve the system.
    UMFPackSolver umf_solver;
    umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
    umf_solver.SetOperator(A);
-   umf_solver.Mult(B, X);
+   umf_solver.Mult(B, x);
 #endif
 
    // 11. Recover the solution as a finite element grid function.
    a->RecoverFEMSolution(X, *b, x);
 
    // 12. Compute and print the L^2 norm of the error.
-   cout << "\n|| F_h - F ||_{L^2} = " << x.ComputeL2Error(F) << '\n' << endl;
+   cout << "\n|| F_h - F ||_{L^2} = " << x.ComputeL2Error(F) << endl;
+   cout << "\n|| F_p - F ||_{L^2} = " << xp.ComputeL2Error(F) << endl;
+   cout << "\n|| F_h - F_p ||_{L^2} = " << x.ComputeL2Error(Xp) << endl;
 
    // 13. Save the refined mesh and the solution. This output can be viewed
    //     later using GLVis: "glvis -m refined.mesh -g sol.gf".
@@ -210,7 +238,8 @@ int main(int argc, char *argv[])
       int  visport   = 19916;
       socketstream sol_sock(vishost, visport);
       sol_sock.precision(8);
-      sol_sock << "solution\n" << *mesh << x << flush;
+      sol_sock << "solution\n" << *mesh << x 
+               << "window_title 'Approximate'\n" << flush;
    }
 
    // 15. Free the used memory.
@@ -237,8 +266,8 @@ void F_exact(const Vector &p, Vector &F)
    double y = p(1);
    // double z = (dim == 3) ? p(2) : 0.0;
 
-   F(0) = cos(kappa*x)*sin(kappa*y);
-   F(1) = cos(kappa*y)*sin(kappa*x);
+   F(0) = x*y*(x+y-1.);//-2. * x * (-1. + x + 2. * y);//cos(kappa*x)*sin(kappa*y);
+   F(1) = 0.;//-2. * (-1. + y) * (-1. + x + 2. * y);//cos(kappa*y)*sin(kappa*x);
    if (dim == 3)
    {
       F(2) = 0.0;
@@ -256,8 +285,8 @@ void f_exact(const Vector &p, Vector &f)
 
    double temp = 1 + 2*kappa*kappa;
 
-   f(0) = temp*cos(kappa*x)*sin(kappa*y);
-   f(1) = temp*cos(kappa*y)*sin(kappa*x);
+   f(0) = - 2.*x - 2.*y;//1.;//temp*cos(kappa*x)*sin(kappa*y);
+   f(1) = 0.;//temp*cos(kappa*y)*sin(kappa*x);
    if (dim == 3)
    {
       f(2) = 0;

fem/bilininteg.cpp

@@ -2055,6 +2055,345 @@ void DivDivIntegrator::AssembleElementMatrix(
    }
 }
 
+void VectorFEDiffusionIntegrator::AssembleElementMatrix(
+   const FiniteElement &el,
+   ElementTransformation &Trans,
+   DenseMatrix &elmat)
+{
+   int dim = el.GetDim();
+   int dof = el.GetDof();
+   double c;
+
+   jshape.SetSize(dof,dim,dim);
+   elmat.SetSize(dof);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+       int order = 2 * el.GetOrder() - 2;
+       ir = &IntRules.Get(el.GetGeomType(), order);
+   }
+
+   elmat = 0.0;
+
+   for (int i = 0; i < ir -> GetNPoints(); i++)
+   {
+       const IntegrationPoint &ip = ir->IntPoint(i);
+
+       el.CalcJShape(ip, jshape);
+
+       Trans.SetIntPoint(&ip);
+       c = ip.weight / Trans.Weight(); // ip.weight * Trans.Weight() ?
+
+       if (Q)
+       {
+          c *= Q -> Eval(Trans, ip);
+       }
+
+       for (int row = 0; row < dof; row++)
+       {
+          for (int col = 0; col < dof; col++)
+          {
+              elmat(row,col) += c*(jshape(row,0,0)*jshape(col,0,0) +
+                                   jshape(row,0,1)*jshape(col,0,1) +
+                                   jshape(row,1,0)*jshape(col,1,0) +
+                                   jshape(row,1,1)*jshape(col,1,1) );
+          }
+       }
+   }
+}
+
+
+void VectorFEDGDiffusionIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1,
+   const FiniteElement &el2,
+   FaceElementTransformations &Trans,
+   DenseMatrix &elmat)
+{
+   int dim,ndof1, ndof2, ndofs;
+   bool kappa_is_nonzero = (kappa !=0.);
+   double w, wq = 0.0;
+
+   dim = el1.GetDim();
+   ndof1 = el1.GetDof();
+
+   if (dim != 2)
+   {
+       mfem_error("This integrator is only for 2 dimensional RT Triangle elements");
+   }
+
+   nor.SetSize(dim);
+   ni.SetSize(dim);
+   adjJ.SetSize(dim);
+
+   shape1.SetSize(ndof1, dim);
+   shape1on.SetSize(ndof1, dim, dim);
+   jshape1.SetSize(ndof1, dim, dim);
+
+   if (Trans.Elem2No >=0)
+   {
+       ndof2 = el2.GetDof();
+       shape2.SetSize(ndof2, dim);
+       shape2on.SetSize(ndof2, dim, dim);
+       jshape2.SetSize(ndof2, dim, dim);
+   }
+   else
+   {
+      ndof2 = 0;
+   }
+
+   ndofs = ndof1 + ndof2;
+   elmat.SetSize(ndofs);
+   elmat = 0.0;
+
+   if (kappa_is_nonzero)
+   {
+      jmat.SetSize(ndofs);
+      jmat = 0.;
+   }
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order;
+      if (ndof2)
+      {
+         order = 2*max(el1.GetOrder(), el2.GetOrder());
+      }
+      else
+      {
+         order = 2*el1.GetOrder();
+      }
+      ir = &IntRules.Get(Trans.FaceGeom, order);
+   }
+
+   for (int p = 0; p < ir -> GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+      IntegrationPoint eip1, eip2;
+
+      Trans.Loc1.Transform(ip, eip1);
+      Trans.Face->SetIntPoint(&ip);
+
+      CalcOrtho(Trans.Face->Jacobian(), nor);
+
+      el1.CalcVShape(eip1, shape1);
+      el1.CalcJShape(eip1, jshape1);
+      Trans.Elem1->SetIntPoint(&eip1);
+      w = ip.weight/Trans.Elem1->Weight();
+
+      if (ndof2)
+      { 
+         w /= 2.;
+      }
+
+      if (Q)
+      { 
+         w *= Q->Eval(*Trans.Elem1, eip1);
+      }
+      ni.Set(w, nor);
+      nh = ni; // Just to initialize
+
+      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(ni, nh);
+      
+      if (kappa_is_nonzero)
+      {
+         wq = ni*nor;
+      }
+      
+      for (int row = 0; row < ndof1; row++)
+      { 
+         shape1on(row,0,0) = shape1(row,0)*nh(0); shape1on(row,0,1) = shape1(row,0)*nh(1);
+         shape1on(row,1,0) = shape1(row,1)*nh(0); shape1on(row,1,1) = shape1(row,1)*nh(1);
+      }
+
+      for (int i = 0; i < ndof1; i++)
+      { 
+         for (int j = 0; j < ndof1; j++)
+         {
+            elmat(i,j) = ( shape1on(i,0,0)*jshape1(j,0,0) +
+                           shape1on(i,0,1)*jshape1(j,0,1) +
+                           shape1on(i,1,0)*jshape1(j,1,0) +
+                           shape1on(i,1,1)*jshape1(j,1,1) );
+         }
+      }
+
+      if (ndof2)
+      { 
+         Trans.Loc2.Transform(ip, eip2);
+         el2.CalcVShape(eip2, shape2);
+         el2.CalcJShape(eip2, jshape2);
+         Trans.Elem2->SetIntPoint(&eip2);
+         w = ip.weight/Trans.Elem2->Weight()/2.; 
+
+         if (Q)
+         {
+            w *= Q->Eval(*Trans.Elem2, eip2);
+         }
+         ni.Set(w, nor);
+
+         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
+         adjJ.Mult(ni, nh);
+         
+         if (kappa_is_nonzero)
+         {
+            wq += ni*nor;
+         }
+         
+         for (int row = 0; row < ndof1; row++)
+         {
+            shape2on(row,0,0) = shape2(row,0)*nh(0); shape2on(row,0,1) = shape2(row,0)*nh(1);
+            shape2on(row,1,0) = shape2(row,1)*nh(0); shape2on(row,1,1) = shape2(row,1)*nh(1);
+         }
+         // To adjust for the normal direction
+         // Add this one (+=) ?
+         for (int i = 0; i < ndof1; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(i, ndof1 + j) += ( shape1on(i,0,0)*jshape2(j,0,0) +
+                                        shape1on(i,0,1)*jshape2(j,0,1) +
+                                        shape1on(i,1,0)*jshape2(j,1,0) +
+                                        shape1on(i,1,1)*jshape2(j,1,1) ); 
+            }
+         // Subtract these two (-=) ?
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof1; j++)
+            {
+               elmat(ndof1 + i, j) -= ( shape2on(i,0,0)*jshape1(j,0,0) +
+                                        shape2on(i,0,1)*jshape1(j,0,1) +
+                                        shape2on(i,1,0)*jshape1(j,1,0) +
+                                        shape2on(i,1,1)*jshape1(j,1,1) ); 
+            }
+
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(ndof1 + i, ndof1 + j) -= ( shape2on(i,0,0)*jshape2(j,0,0) +
+                                                shape2on(i,0,1)*jshape2(j,0,1) +
+                                                shape2on(i,1,0)*jshape2(j,1,0) +
+                                                shape2on(i,1,1)*jshape2(j,1,1) ); 
+            }
+      }
+      
+      if (kappa_is_nonzero)
+      {
+ //////////////////////////////////////////     
+ //////////////////////////////     
+////////////////////////////////////////////////      
+      w = ip.weight/Trans.Elem1->Weight();
+
+      if (ndof2)
+      { 
+         w /= 2.;
+      }
+      
+      nh = nor; // Just to initialize
+
+      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(nor, nh);
+      
+      for (int row = 0; row < ndof1; row++)
+      { 
+         shape1on(row,0,0) = shape1(row,0)*nh(0); shape1on(row,0,1) = shape1(row,0)*nh(1);
+         shape1on(row,1,0) = shape1(row,1)*nh(0); shape1on(row,1,1) = shape1(row,1)*nh(1);
+      }
+///////////////////////////////////////////      
+ ///////////////////////////////////////////     
+//////////////////////////////////////////      
+      
+      
+         // only assemble the lower triangular part of jmat
+         wq = w*kappa;
+         for (int i = 0; i < ndof1; i++)
+         {
+            for (int j = 0; j <= i; j++)
+            {
+               jmat(i,j) += wq * ( shape1on(i,0,0) * shape1on(j,0,0) +
+                                   shape1on(i,0,1) * shape1on(j,0,1) +
+                                   shape1on(i,1,0) * shape1on(j,1,0) +
+                                   shape1on(i,1,1) * shape1on(j,1,1) );
+//                 jmat(i,j) += wq * ( shape1(i,0) * shape1(j,0) +
+//                                     shape1(i,1) * shape1(j,1) );
+            }
+         }
+         if (ndof2)
+         {
+         
+////////////////////////////         
+////////////////////////////
+////////////////////////////
+         
+         w = ip.weight/Trans.Elem2->Weight()/2.; 
+
+         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
+         adjJ.Mult(nor, nh);
+         
+         for (int row = 0; row < ndof1; row++)
+         {
+            shape2on(row,0,0) = shape2(row,0)*nh(0); shape2on(row,0,1) = shape2(row,0)*nh(1);
+            shape2on(row,1,0) = shape2(row,1)*nh(0); shape2on(row,1,1) = shape2(row,1)*nh(1);
+         }
+///////////////////////////
+///////////////////////////
+///////////////////////////
+         
+         
+            for (int i = 0; i < ndof2; i++)
+            {
+               for (int j = 0; j < ndof1; j++)
+               {
+                  jmat(ndof1+i, j) -= wq * ( shape2on(i,0,0) * shape1on(j,0,0) +
+                                             shape2on(i,0,1) * shape1on(j,0,1) +
+                                             shape2on(i,1,0) * shape1on(j,1,0) +
+                                             shape2on(i,1,1) * shape1on(j,1,1) );
+//                 jmat(ndof1+i,j) -= wq * ( shape2(i,0) * shape1(j,0) +
+//                                           shape2(i,1) * shape1(j,1) );
+               }
+               for (int j = 0; j <= i; j++)
+               {
+                  jmat(ndof1+i, ndof1+j) += wq * ( shape2on(i,0,0) * shape2on(j,0,0) +
+                                                   shape2on(i,0,1) * shape2on(j,0,1) +
+                                                   shape2on(i,1,0) * shape2on(j,1,0) +
+                                                   shape2on(i,1,1) * shape2on(j,1,1) );
+//                 jmat(ndof1+i,ndof1+j) += wq * ( shape2(i,0) * shape2(j,0) +
+//                                                 shape2(i,1) * shape2(j,1) );
+               }
+            }
+         }
+      } 
+      
+      if (kappa_is_nonzero)
+      {
+         for (int i = 0; i < ndofs; i++)
+         {
+            for (int j = 0; j < i; j++)
+            {
+               double aij = elmat(i,j), aji = elmat(j,i), mij = jmat(i,j);
+               elmat(i,j) = sigma*aji - aij + mij;
+               elmat(j,i) = sigma*aij - aji + mij;
+            }
+            elmat(i,i) = (sigma - 1.)*elmat(i,i) + jmat(i,i);
+         }
+      }
+      else
+      {
+         for (int i = 0; i < ndofs; i++)
+         {
+            for (int j = 0; j < i; j++)
+            {
+               double aij = elmat(i,j), aji = elmat(j,i);
+               elmat(i,j) = sigma*aji - aij;
+               elmat(j,i) = sigma*aij - aji;
+            }
+            elmat(i,i) *= (sigma - 1.);
+         }
+      }
+//      elmat.Threshold(1.e-16);
+//      elmat.Print();
+   }
+
+}
 
 void VectorDiffusionIntegrator::AssembleElementMatrix(
    const FiniteElement &el,

fem/bilininteg.hpp

@@ -1978,6 +1978,50 @@ public:
                                       DenseMatrix &elmat);
 };
 
+/// \int_{K} Q \grad u : \grad v for RT elements
+class VectorFEDiffusionIntegrator: public BilinearFormIntegrator
+{
+private:
+   Coefficient *Q;
+
+#ifndef MFEM_THREAD_SAFE
+   DenseTensor jshape;
+#endif
+
+public:
+   VectorFEDiffusionIntegrator() { Q = NULL; }
+   VectorFEDiffusionIntegrator(Coefficient &q) : Q(&q) { }
+
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+};
+
+/// DG Diffusion integrator for RT elements
+class VectorFEDGDiffusionIntegrator: public BilinearFormIntegrator
+{
+private:
+   Coefficient *Q;
+   double sigma, kappa;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector nor, ni, nh;
+   DenseMatrix shape1, shape2, jmat, adjJ;
+   DenseTensor jshape1, jshape2, shape1on, shape2on;
+#endif
+
+public:
+   VectorFEDGDiffusionIntegrator(const double s, const double k)
+	   : Q(NULL), sigma(s), kappa(k) { }
+   VectorFEDGDiffusionIntegrator(Coefficient &q, const double s, const double k)
+	   : Q(&q), sigma(s), kappa(k) { }
+
+   virtual void AssembleFaceMatrix(const FiniteElement &el1,
+		                   const FiniteElement &el2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &elmat);
+};
+
 /** Integrator for
       (Q grad u, grad v) = sum_i (Q grad u_i, grad v_i) e_i e_i^T
     for FE spaces defined by 'dim' copies of a scalar FE space. Where e_i

fem/fe.cpp

@@ -58,6 +58,13 @@ void FiniteElement::CalcDivShape (
                "   is not implemented for this class!");
 }
 
+void FiniteElement::CalcJShape (
+   const IntegrationPoint &ip, DenseTensor &jshape) const
+{
+   mfem_error ("FiniteElement::CalcJShape (ip, ...)\n"
+   	       "   is not implemented for this class!");
+}
+
 void FiniteElement::CalcPhysDivShape(
    ElementTransformation &Trans, Vector &div_shape) const
 {
@@ -2665,6 +2672,21 @@ void RT0TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
    divshape(2) = 2.;
 }
 
+void RT0TriangleFiniteElement::CalcJShape(const IntegrationPoint &ip,
+                                          DenseTensor &jshape) const
+{
+   double x = ip.x, y = ip.y;
+
+   jshape(0,0,0) = 1.; jshape(0,0,1) = 0.;
+   jshape(0,1,0) = 0.; jshape(0,1,1) = 1.;
+
+   jshape(1,0,0) = 1.; jshape(1,0,1) = 0.;
+   jshape(1,1,0) = 0.; jshape(1,1,1) = 1.;
+
+   jshape(2,0,0) = 1.; jshape(2,0,1) = 0.;
+   jshape(2,1,0) = 0.; jshape(2,1,1) = 1.;
+}
+
 const double RT0TriangleFiniteElement::nk[3][2] =
 { {0, -1}, {1, 1}, {-1, 0} };
 
@@ -2884,22 +2906,22 @@ void RT1TriangleFiniteElement::CalcVShape(const IntegrationPoint &ip,
 {
    double x = ip.x, y = ip.y;
 
-   shape(0,0) = -2 * x * (-1 + x + 2 * y);
-   shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
-   shape(1,0) =  2 * x * (x - y);
-   shape(1,1) =  2 * (x - y) * (-1 + y);
-   shape(2,0) =  2 * x * (-1 + 2 * x + y);
-   shape(2,1) =  2 * y * (-1 + 2 * x + y);
-   shape(3,0) =  2 * x * (-1 + x + 2 * y);
-   shape(3,1) =  2 * y * (-1 + x + 2 * y);
-   shape(4,0) = -2 * (-1 + x) * (x - y);
-   shape(4,1) =  2 * y * (-x + y);
-   shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
-   shape(5,1) = -2 * y * (-1 + 2 * x + y);
-   shape(6,0) = -3 * x * (-2 + 2 * x + y);
-   shape(6,1) = -3 * y * (-1 + 2 * x + y);
-   shape(7,0) = -3 * x * (-1 + x + 2 * y);
-   shape(7,1) = -3 * y * (-2 + x + 2 * y);
+   shape(0,0) = -2. * x * (-1. + x + 2. * y);
+   shape(0,1) = -2. * (-1. + y) * (-1. + x + 2. * y);
+   shape(1,0) =  2. * x * (x - y);
+   shape(1,1) =  2. * (x - y) * (-1. + y);
+   shape(2,0) =  2. * x * (-1. + 2. * x + y);
+   shape(2,1) =  2. * y * (-1. + 2. * x + y);
+   shape(3,0) =  2. * x * (-1. + x + 2. * y);
+   shape(3,1) =  2. * y * (-1. + x + 2. * y);
+   shape(4,0) = -2. * (-1. + x) * (x - y);
+   shape(4,1) =  2. * y * (-x + y);
+   shape(5,0) = -2. * (-1. + x) * (-1. + 2. * x + y);
+   shape(5,1) = -2. * y * (-1. + 2. * x + y);
+   shape(6,0) = -3. * x * (-2. + 2. * x + y);
+   shape(6,1) = -3. * y * (-1. + 2. * x + y);
+   shape(7,0) = -3. * x * (-1. + x + 2. * y);
+   shape(7,1) = -3. * y * (-2. + x + 2. * y);
 }
 
 void RT1TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
@@ -2907,14 +2929,44 @@ void RT1TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
 {
    double x = ip.x, y = ip.y;
 
-   divshape(0) = -2 * (-4 + 3 * x + 6 * y);
-   divshape(1) =  2 + 6 * x - 6 * y;
-   divshape(2) = -4 + 12 * x + 6 * y;
-   divshape(3) = -4 + 6 * x + 12 * y;
-   divshape(4) =  2 - 6 * x + 6 * y;
-   divshape(5) = -2 * (-4 + 6 * x + 3 * y);
-   divshape(6) = -9 * (-1 + 2 * x + y);
-   divshape(7) = -9 * (-1 + x + 2 * y);
+   divshape(0) = -2. * (-4. + 3. * x + 6. * y);
+   divshape(1) =  2. + 6. * x - 6. * y;
+   divshape(2) = -4. + 12. * x + 6. * y;
+   divshape(3) = -4. + 6. * x + 12. * y;
+   divshape(4) =  2. - 6. * x + 6. * y;
+   divshape(5) = -2. * (-4. + 6. * x + 3. * y);
+   divshape(6) = -9. * (-1. + 2. * x + y);
+   divshape(7) = -9. * (-1. + x + 2. * y);
+}
+
+void RT1TriangleFiniteElement::CalcJShape(const IntegrationPoint &ip,
+                                          DenseTensor &jshape) const
+{
+   double x = ip.x, y = ip.y;
+
+   jshape(0,0,0) = 2. - 4.*y - 4.*x; jshape(0,0,1) = -4.*x;
+   jshape(0,1,0) = 2. - 2.*y; jshape(0,1,1) = 6. - 8.*y - 2.*x;
+
+   jshape(1,0,0) = 4.*x - 2.*y; jshape(1,0,1) = -2.*x;
+   jshape(1,1,0) = 2.*y - 2.; jshape(1,1,1) = 2.*x - 4.*y + 2.;
+
+   jshape(2,0,0) = 8.*x + 2.*y - 2.; jshape(2,0,1) = 2.*x;
+   jshape(2,1,0) = 4.*y; jshape(2,1,1) = 4.*x + 4.*y - 2.;
+
+   jshape(3,0,0) = 4.*x + 4.*y - 2.; jshape(3,0,1) = 4.*x;
+   jshape(3,1,0) = 2.*y; jshape(3,1,1) = 2.*x + 8.*y - 2.;
+
+   jshape(4,0,0) = 2.*y - 4.*x + 2.; jshape(4,0,1) = 2.*x - 2.;
+   jshape(4,1,0) = -2.*y; jshape(4,1,1) = 4.*y - 2.*x;
+
+   jshape(5,0,0) = 6. - 2.*y - 8.*x; jshape(5,0,1) = 2. - 2.*x;
+   jshape(5,1,0) = -4.*y; jshape(5,1,1) = 2. - 4.*y - 4.*x;
+
+   jshape(6,0,0) = 6. - 3.*y - 12.*x; jshape(6,0,1) = -3.*x;
+   jshape(6,1,0) = -6.*y; jshape(6,1,1) = 3. - 6.*y - 6.*y;
+
+   jshape(7,0,0) = 3. - 6.*y - 6.*x; jshape(7,0,1) = -6.*x;
+   jshape(7,1,0) = -3.*y; jshape(7,1,1) = 6. - 12.*y - 3.*x;
 }
 
 const double RT1TriangleFiniteElement::nk[8][2] =
@@ -3320,6 +3372,41 @@ void RT2TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
    }
 }
 
+void RT2TriangleFiniteElement::CalcJShape(const IntegrationPoint &ip,
+                                          DenseTensor &jshape) const
+{
+   double x = ip.x, y = ip.y;
+
+   double Bxx[15] = {0., 0., 1., 0., 0., 0., 2.*x, 0., y, 0., 0., 0., 3.*x*x,
+                     2.*x*y, y*y
+                    };
+           
+   double Bxy[15] = {0., 0., 0., 0., 1., 0., 0., 0., x, 0., 2.*y, 0., 0., x*x, 
+                     2.*x*y
+                    };
+                   
+   double Byx[15] = {0., 0., 0., 1., 0., 0., 0., 2.*x, 0., y, 0., 0., 2.*x*y,
+                     y*y, 0.
+                    };
+                   
+   double Byy[15] = {0., 0., 0., 0., 0., 1., 0., 0., 0., x, 0., 2.*y, x*x,
+                     2.*x*y, 3.*y*y
+                    };
+
+   for (int i = 0; i < 15; i++)
+   {
+      double cxx = 0.0, cxy = 0.0, cyx = 0.0, cyy = 0.0;
+      for (int j = 0; j < 15; j++)
+      {
+         cxx += M[i][j] * Bxx[j]; cxy += M[i][j] * Bxy[j];
+         cyx += M[i][j] * Byx[j]; cyy += M[i][j] * Byy[j];
+      }
+      jshape(i,0,0) = cxx; jshape(i,0,1) = cxy;
+      jshape(i,1,0) = cyx; jshape(i,1,1) = cyy;
+   }
+
+}
+
 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
 
 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
@@ -9826,6 +9913,51 @@ void RT_TriangleElement::CalcDivShape(const IntegrationPoint &ip,
    Ti.Mult(divu, divshape);
 }
 
+void RT_TriangleElement::CalcJShape(const IntegrationPoint &ip,
+                                          DenseTensor &jshape) const
+{
+   const int p = Order - 1;
+
+   Vector shape_x(p + 1),  shape_y(p + 1),  shape_l(p + 1);
+   Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
+   DenseTensor ju(Dof, Dim, Dim);
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)
+      {
+         int k = p - i - j;
+         //double s = shape_x(i)*shape_y(j)*shape_l(p-i-j);
+         ju(o,0,0) = shape_y(j)*(dshape_x(i)*shape_l(k)-shape_x(i)*dshape_l(k));
+         ju(o,0,1) = shape_x(i)*(dshape_y(j)*shape_l(k)-shape_y(j)*dshape_l(k));
+         ju(o,1,0) = 0.;
+         ju(o,1,1) = 0.;
+         o++;
+         ju(o,0,0) = 0.;
+         ju(o,0,1) = 0.;
+         ju(o,1,0) = shape_y(j)*(dshape_x(i)*shape_l(k)-shape_x(i)*dshape_l(k));
+         ju(o,1,1) = shape_x(i)*(dshape_y(j)*shape_l(k)-shape_y(j)*dshape_l(k));
+         o++;
+      }
+   for (int i = 0; i <= p; i++)
+   {
+      int j = p - i;
+      double s = shape_x(i)*shape_y(j);
+      double sx = dshape_x(i)*shape_y(j);
+      double sy = shape_x(i)*dshape_y(j);
+      ju(o,0,0) = s + (ip.x - c)*sx;
+      ju(o,0,1) = (ip.x - c)*sy;
+      ju(o,1,0) = (ip.y - c)*sx;
+      ju(o,1,1) = s + (ip.y - c)*sy;
+      o++;
+   }
+
+   Ti.Mult(ju, jshape);
+}
 
 const double RT_TetrahedronElement::nk[12] =
 { 1,1,1,  -1,0,0,  0,-1,0,  0,0,-1 };

fem/fe.hpp

@@ -271,6 +271,9 @@ public:
 
    // virtual functions for finite elements on vector spaces
 
+   virtual void CalcJShape(const IntegrationPoint &ip,
+                           DenseTensor &jshape) const;
+
    /** @brief Evaluate the values of all shape functions of a *vector* finite
        element in reference space at the given point @a ip. */
    /** Each row of the result DenseMatrix @a shape contains the components of
@@ -1066,6 +1069,9 @@ public:
 
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
+                             
+   virtual void CalcJShape(const IntegrationPoint &ip,
+                           DenseTensor &jshape) const;
 
    virtual void GetLocalInterpolation (ElementTransformation &Trans,
                                        DenseMatrix &I) const;
@@ -1120,6 +1126,9 @@ public:
 
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
+                             
+   virtual void CalcJShape(const IntegrationPoint &ip,
+                           DenseTensor &jshape) const;
 
    virtual void GetLocalInterpolation (ElementTransformation &Trans,
                                        DenseMatrix &I) const;
@@ -1173,6 +1182,9 @@ public:
 
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
+                             
+   virtual void CalcJShape(const IntegrationPoint &ip,
+                           DenseTensor &jshape) const;
 };
 
 class RT2QuadFiniteElement : public VectorFiniteElement
@@ -2204,6 +2216,7 @@ class RT_TriangleElement : public VectorFiniteElement
    mutable Vector shape_x, shape_y, shape_l;
    mutable Vector dshape_x, dshape_y, dshape_l;
    mutable DenseMatrix u;
+   mutable DenseTensor ju;
    mutable Vector divu;
 #endif
    Array<int> dof2nk;
@@ -2218,6 +2231,10 @@ public:
    { CalcVShape_RT(Trans, shape); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
+                             
+   virtual void CalcJShape(const IntegrationPoint &ip,
+                           DenseTensor &jshape) const;
+                           
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
                                       DenseMatrix &I) const
    { LocalInterpolation_RT(*this, nk, dof2nk, Trans, I); }

fem/fe_coll.cpp

@@ -1468,9 +1468,6 @@ const
 
 H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
 {
-   MFEM_VERIFY(p >= 1, "H1_FECollection requires order >= 1.");
-   MFEM_VERIFY(dim >= 0 && dim <= 3, "H1_FECollection requires 0 <= dim <= 3.");
-
    const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1;
 
    int pt_type = BasisType::GetQuadrature1D(btype);
@@ -1711,8 +1708,6 @@ H1_Trace_FECollection::H1_Trace_FECollection(const int p, const int dim,
 L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
                                  const int map_type)
 {
-   MFEM_VERIFY(p >= 0, "L2_FECollection requires order >= 0.");
-
    b_type = BasisType::Check(btype);
    const char *prefix = NULL;
    switch (map_type)
@@ -1887,8 +1882,6 @@ RT_FECollection::RT_FECollection(const int p, const int dim,
                                  const int cb_type, const int ob_type)
    : ob_type(ob_type)
 {
-   MFEM_VERIFY(p >= 0, "RT_FECollection requires order >= 0.");
-
    int cp_type = BasisType::GetQuadrature1D(cb_type);
    int op_type = BasisType::GetQuadrature1D(ob_type);
 
@@ -2164,9 +2157,6 @@ DG_Interface_FECollection::DG_Interface_FECollection(const int p, const int dim,
 ND_FECollection::ND_FECollection(const int p, const int dim,
                                  const int cb_type, const int ob_type)
 {
-   MFEM_VERIFY(p >= 1, "ND_FECollection requires order >= 1.");
-   MFEM_VERIFY(dim >= 1 && dim <= 3, "ND_FECollection requires 1 <= dim <= 3.");
-
    const int pm1 = p - 1, pm2 = p - 2;
 
    if (cb_type == BasisType::GaussLobatto &&

fem/intrules.hpp

@@ -31,7 +31,6 @@ public:
 
    void Set(const double *p, const int dim)
    {
-      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
       x = p[0];
       if (dim > 1)
       {
@@ -45,7 +44,6 @@ public:
 
    void Get(double *p, const int dim) const
    {
-      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
       p[0] = x;
       if (dim > 1)
       {

general/sort_pairs.hpp

@@ -58,11 +58,6 @@ public:
    A one;
    B two;
    C three;
-
-   Triple() { }
-
-   Triple(const A &one, const B &two, const C &three)
-      : one(one), two(two), three(three) { }
 };
 
 /// @brief Lexicographic comparison operator for class Triple.

linalg/densemat.cpp

@@ -4273,6 +4273,12 @@ void DenseMatrixInverse::Mult(const DenseMatrix &B, DenseMatrix &X) const
    lu.Solve(width, X.Width(), X.Data());
 }
 
+void DenseMatrixInverse::Mult(const DenseTensor &B, DenseTensor &X) const
+{
+   X = B;
+   lu.Solve(width, X.SizeJ()*X.SizeK(), X.Data());
+}
+
 void DenseMatrixInverse::TestInversion()
 {
    DenseMatrix C(width);
@@ -4487,4 +4493,19 @@ DenseTensor &DenseTensor::operator=(double c)
    return *this;
 }
 
+DenseTensor &DenseTensor::operator=(const DenseTensor &t)
+{
+   SetSize(t.SizeI(), t.SizeJ(), t.SizeK());
+
+   const int hwd = t.SizeI()*t.SizeJ()*t.SizeK();
+   tdata = new double[hwd];
+   std::memcpy(tdata, t.tdata, sizeof(double) * hwd);
+//   for (int i = 0; i < hwd; i++)
+//   {
+//      tdata[i] = t.tdata[i];
+//   }
+
+   return *this;
+}
+
 }

linalg/densemat.hpp

@@ -543,6 +543,107 @@ public:
                        const double *X2, double *Y1) const;
 };
 
+class Table;
+
+/// Rank 3 tensor (array of matrices)
+class DenseTensor
+{
+private:
+   DenseMatrix Mk;
+   double *tdata;
+   int nk;
+   bool own_data;
+
+public:
+   DenseTensor()
+   {
+      nk = 0;
+      tdata = NULL;
+      own_data = true;
+   }
+
+   DenseTensor(int i, int j, int k)
+      : Mk(NULL, i, j)
+   {
+      nk = k;
+      tdata = new double[i*j*k];
+      own_data = true;
+   }
+
+   /// Copy constructor: deep copy
+   DenseTensor(const DenseTensor& other)
+      : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk), own_data(true)
+   {
+      const int size = Mk.Height()*Mk.Width()*nk;
+      if (size > 0)
+      {
+         tdata = new double[size];
+         std::memcpy(tdata, other.tdata, sizeof(double) * size);
+      }
+      else
+      {
+         tdata = NULL;
+      }
+   }
+
+   int SizeI() const { return Mk.Height(); }
+   int SizeJ() const { return Mk.Width(); }
+   int SizeK() const { return nk; }
+
+   void SetSize(int i, int j, int k)
+   {
+      if (own_data) { delete [] tdata; }
+      Mk.UseExternalData(NULL, i, j);
+      nk = k;
+      tdata = new double[i*j*k];
+      own_data = true;
+   }
+
+   void UseExternalData(double *ext_data, int i, int j, int k)
+   {
+      if (own_data) { delete [] tdata; }
+      Mk.UseExternalData(NULL, i, j);
+      nk = k;
+      tdata = ext_data;
+      own_data = false;
+   }
+
+   /// Sets the tensor elements equal to constant c
+   DenseTensor &operator=(double c);
+   
+   /// Sets the tensor size and elements equal to another tensor
+   DenseTensor &operator=(const DenseTensor &t);
+
+   DenseMatrix &operator()(int k) { Mk.data = GetData(k); return Mk; }
+   const DenseMatrix &operator()(int k) const
+   { return const_cast<DenseTensor&>(*this)(k); }
+
+   double &operator()(int i, int j, int k)
+   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }
+   const double &operator()(int i, int j, int k) const
+   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }
+
+   double *GetData(int k) { return tdata+k*Mk.Height()*Mk.Width(); }
+
+   double *Data() { 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;
+
+   /** Tensor contraction on the last two indices */
+   void DoubleDot(const DenseTensor &T, DenseMatrix &A);
+
+   void Clear()
+   { UseExternalData(NULL, 0, 0, 0); }
+
+   long MemoryUsage() const { return nk*Mk.MemoryUsage(); }
+
+   ~DenseTensor()
+   {
+      if (own_data) { delete [] tdata; }
+   }
+};
 
 /** Data type for inverse of square dense matrix.
     Stores LU factors */
@@ -580,6 +681,9 @@ public:
    /// Multiply the inverse matrix by another matrix: X = A^{-1} B.
    void Mult(const DenseMatrix &B, DenseMatrix &X) const;
 
+   /// Multiply the inverse matrix by a tensor: X = A^{-1} B.
+   void Mult(const DenseTensor &B, DenseTensor &X) const;
+
    /// Compute and return the inverse matrix in Ainv.
    void GetInverseMatrix(DenseMatrix &Ainv) const
    {
@@ -651,103 +755,6 @@ public:
    ~DenseMatrixSVD();
 };
 
-class Table;
-
-/// Rank 3 tensor (array of matrices)
-class DenseTensor
-{
-private:
-   DenseMatrix Mk;
-   double *tdata;
-   int nk;
-   bool own_data;
-
-public:
-   DenseTensor()
-   {
-      nk = 0;
-      tdata = NULL;
-      own_data = true;
-   }
-
-   DenseTensor(int i, int j, int k)
-      : Mk(NULL, i, j)
-   {
-      nk = k;
-      tdata = new double[i*j*k];
-      own_data = true;
-   }
-
-   /// Copy constructor: deep copy
-   DenseTensor(const DenseTensor& other)
-      : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk), own_data(true)
-   {
-      const int size = Mk.Height()*Mk.Width()*nk;
-      if (size > 0)
-      {
-         tdata = new double[size];
-         std::memcpy(tdata, other.tdata, sizeof(double) * size);
-      }
-      else
-      {
-         tdata = NULL;
-      }
-   }
-
-   int SizeI() const { return Mk.Height(); }
-   int SizeJ() const { return Mk.Width(); }
-   int SizeK() const { return nk; }
-
-   void SetSize(int i, int j, int k)
-   {
-      if (own_data) { delete [] tdata; }
-      Mk.UseExternalData(NULL, i, j);
-      nk = k;
-      tdata = new double[i*j*k];
-      own_data = true;
-   }
-
-   void UseExternalData(double *ext_data, int i, int j, int k)
-   {
-      if (own_data) { delete [] tdata; }
-      Mk.UseExternalData(NULL, i, j);
-      nk = k;
-      tdata = ext_data;
-      own_data = false;
-   }
-
-   /// Sets the tensor elements equal to constant c
-   DenseTensor &operator=(double c);
-
-   DenseMatrix &operator()(int k) { Mk.data = GetData(k); return Mk; }
-   const DenseMatrix &operator()(int k) const
-   { return const_cast<DenseTensor&>(*this)(k); }
-
-   double &operator()(int i, int j, int k)
-   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }
-   const double &operator()(int i, int j, int k) const
-   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }
-
-   double *GetData(int k) { return tdata+k*Mk.Height()*Mk.Width(); }
-
-   double *Data() { 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()
-   { UseExternalData(NULL, 0, 0, 0); }
-
-   long MemoryUsage() const { return nk*Mk.MemoryUsage(); }
-
-   ~DenseTensor()
-   {
-      if (own_data) { delete [] tdata; }
-   }
-};
-
-
 // Inline methods
 
 inline double &DenseMatrix::operator()(int i, int j)

linalg/sparsemat.hpp

@@ -406,6 +406,9 @@ public:
    /// Print various sparse matrix staticstics.
    void PrintInfo(std::ostream &out) const;
 
+   /// Walks the sparse matrix
+   int Walk(int &i, int &j, double &a);
+
    /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
    double IsSymmetric() const;
 

mesh/element.hpp

@@ -17,7 +17,6 @@
 #include "../general/table.hpp"
 #include "../linalg/densemat.hpp"
 #include "../fem/geom.hpp"
-#include "../general/hash.hpp"
 
 namespace mfem
 {
@@ -81,7 +80,7 @@ public:
    virtual void MarkEdge(const DSTable &v_to_v, const int *length) {}
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const { return 0; }
+   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const { return 0; }
 
    /// Set current coarse-fine transformation number.
    virtual void ResetTransform(int tr) {}

mesh/mesh.cpp

@@ -2389,7 +2389,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
       faces[i] = (face) ? face->Duplicate(this) : NULL;
    }
    mesh.faces_info.Copy(faces_info);
-   mesh.nc_faces_info.Copy(nc_faces_info);
 
    // Do NOT copy the element-to-element Table, el_to_el
    el_to_el = NULL;
@@ -2419,16 +2418,9 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    }
 
    // Deep copy the NCMesh.
-#ifdef MFEM_USE_MPI
-   if (dynamic_cast<const ParMesh*>(&mesh))
-   {
-      ncmesh = NULL; // skip; will be done in ParMesh copy ctor
-   }
-   else
-#endif
-   {
-      ncmesh = mesh.ncmesh ? new NCMesh(*mesh.ncmesh) : NULL;
-   }
+   // TODO: ParNCMesh; ParMesh has a separate 'pncmesh' pointer, and 'ncmesh'
+   //       is initialized from it. Need ParNCMesh copy constructor.
+   ncmesh = mesh.ncmesh ? new NCMesh(*mesh.ncmesh) : NULL;
 
    // Duplicate the Nodes, including the FiniteElementCollection and the
    // FiniteElementSpace
@@ -5913,47 +5905,56 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
    }
    else if (Dim == 3) // ---------------------------------------------------
    {
-      // 1. Hash table of vertex to vertex connections corresponding to refined
-      //    edges.
-      HashTable<Hashed2> v_to_v;
+      // 1. Get table of vertex to vertex connections.
+      DSTable v_to_v(NumOfVertices);
+      GetVertexToVertexTable(v_to_v);
 
-      // 2. Do the red refinement.
+      // 2. Get edge to element connections in arrays edge1 and edge2
+      nedges = v_to_v.NumberOfEntries();
+      int *middle = new int[nedges];
+
+      for (i = 0; i < nedges; i++)
+      {
+         middle[i] = -1;
+      }
+
+      // 3. Do the red refinement.
       int ii;
       switch (type)
       {
          case 1:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
          case 2:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
 
-               Bisection(NumOfElements - 1, v_to_v);
-               Bisection(marked_el[i], v_to_v);
+               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
          case 3:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
 
                ii = NumOfElements - 1;
-               Bisection(ii, v_to_v);
-               Bisection(NumOfElements - 1, v_to_v);
-               Bisection(ii, v_to_v);
+               Bisection(ii, v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
+               Bisection(ii, v_to_v, NULL, NULL, middle);
 
-               Bisection(marked_el[i], v_to_v);
-               Bisection(NumOfElements-1, v_to_v);
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
       }
 
-      // 3. Do the green refinement (to get conforming mesh).
+      // 4. Do the green refinement (to get conforming mesh).
       int need_refinement;
       // int need_refinement, onoe, max_gen = 0;
       do
@@ -5967,10 +5968,10 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
             // ((Tetrahedron *)elements[i])->
             // ParseRefinementFlag(redges, type, flag);
             // if (flag > max_gen)  max_gen = flag;
-            if (elements[i]->NeedRefinement(v_to_v))
+            if (elements[i]->NeedRefinement(v_to_v, middle))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v);
+               Bisection(i, v_to_v, NULL, NULL, middle);
             }
          }
       }
@@ -5978,23 +5979,38 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       // mfem::out << "Maximum generation: " << max_gen << endl;
 
-      // 4. Update the boundary elements.
+      // 5. Update the boundary elements.
       do
       {
          need_refinement = 0;
          for (i = 0; i < NumOfBdrElements; i++)
-            if (boundary[i]->NeedRefinement(v_to_v))
+            if (boundary[i]->NeedRefinement(v_to_v, middle))
             {
                need_refinement = 1;
-               BdrBisection(i, v_to_v);
+               Bisection(i, v_to_v, middle);
             }
       }
       while (need_refinement == 1);
 
-      NumOfVertices = vertices.Size();
+      // 6. Un-mark the Pf elements.
+      int refinement_edges[2], type, flag;
+      for (i = 0; i < NumOfElements; i++)
+      {
+         Tetrahedron* el = (Tetrahedron*) elements[i];
+         el->ParseRefinementFlag(refinement_edges, type, flag);
+
+         if (type == Tetrahedron::TYPE_PF)
+         {
+            el->CreateRefinementFlag(refinement_edges, Tetrahedron::TYPE_PU,
+                                     flag);
+         }
+      }
+
       NumOfBdrElements = boundary.Size();
 
-      // 5. Update element-to-edge and element-to-face relations.
+      // 7. Free the allocated memory.
+      delete [] middle;
+
       DeleteLazyTables();
       if (el_to_edge != NULL)
       {
@@ -6572,21 +6588,7 @@ void Mesh::Bisection(int i, const DSTable &v_to_v,
       }
       NumOfElements++;
    }
-   else
-   {
-      MFEM_ABORT("Bisection for now works only for triangles.");
-   }
-}
-
-void Mesh::Bisection(int i, HashTable<Hashed2> &v_to_v)
-{
-   int *vert;
-   int v[2][4], v_new, bisect, t;
-   Element *el = elements[i];
-   Vertex V;
-
-   t = el->GetType();
-   if (t == Element::TETRAHEDRON)
+   else if (t == Element::TETRAHEDRON)
    {
       int j, type, new_type, old_redges[2], new_redges[2][2], flag;
       Tetrahedron *tet = (Tetrahedron *) el;
@@ -6597,19 +6599,32 @@ void Mesh::Bisection(int i, HashTable<Hashed2> &v_to_v)
       vert = tet->GetVertices();
 
       // 1. Get the index for the new vertex in v_new.
-      bisect = v_to_v.FindId(vert[0], vert[1]);
+      bisect = v_to_v(vert[0], vert[1]);
       if (bisect == -1)
       {
-         v_new = NumOfVertices + v_to_v.GetId(vert[0],vert[1]);
+         tet->ParseRefinementFlag(old_redges, type, flag);
+         mfem::err << "Error in Bisection(...) of tetrahedron!" << endl
+                   << "   redge[0] = " << old_redges[0]
+                   << "   redge[1] = " << old_redges[1]
+                   << "   type = " << type
+                   << "   flag = " << flag << endl;
+         mfem_error();
+      }
+
+      if (middle[bisect] == -1)
+      {
+         v_new = NumOfVertices++;
          for (j = 0; j < 3; j++)
          {
             V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
          }
          vertices.Append(V);
+
+         middle[bisect] = v_new;
       }
       else
       {
-         v_new = NumOfVertices + bisect;
+         v_new = middle[bisect];
       }
 
       // 2. Set the node indices for the new elements in v[2][4] so that
@@ -6693,11 +6708,11 @@ void Mesh::Bisection(int i, HashTable<Hashed2> &v_to_v)
    }
    else
    {
-      MFEM_ABORT("Bisection with HashTable for now works only for tetrahedra.");
+      MFEM_ABORT("Bisection for now works only for triangles & tetrahedra.");
    }
 }
 
-void Mesh::BdrBisection(int i, const HashTable<Hashed2> &v_to_v)
+void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
 {
    int *vert;
    int v[2][3], v_new, bisect, t;
@@ -6711,9 +6726,9 @@ void Mesh::BdrBisection(int i, const HashTable<Hashed2> &v_to_v)
       vert = tri->GetVertices();
 
       // 1. Get the index for the new vertex in v_new.
-      bisect = v_to_v.FindId(vert[0], vert[1]);
+      bisect = v_to_v(vert[0], vert[1]);
       MFEM_ASSERT(bisect >= 0, "");
-      v_new = NumOfVertices + bisect;
+      v_new = middle[bisect];
       MFEM_ASSERT(v_new != -1, "");
 
       // 2. Set the node indices for the new elements in v[0] and v[1] so that
@@ -6729,8 +6744,7 @@ void Mesh::BdrBisection(int i, const HashTable<Hashed2> &v_to_v)
    }
    else
    {
-      MFEM_ABORT("Bisection of boundary elements with HashTable works only for"
-                 " triangles!");
+      MFEM_ABORT("Bisection of boundary elements works only for triangles!");
    }
 }
 

mesh/mesh.hpp

@@ -250,14 +250,11 @@ protected:
                         int *edge1, int *edge2, int *middle)
    { Bisection(i, v_to_v, edge1, edge2, middle); }
 
-   /// Bisect a triangle: element with index @a i is bisected.
+   /** Bisection. Element with index i is bisected. */
    void Bisection(int i, const DSTable &, int *, int *, int *);
 
-   /// Bisect a tetrahedron: element with index @a i is bisected.
-   void Bisection(int i, HashTable<Hashed2> &);
-
-   /// Bisect a boundary triangle: boundary element with index @a i is bisected.
-   void BdrBisection(int i, const HashTable<Hashed2> &);
+   /** Bisection. Boundary element with index i is bisected. */
+   void Bisection(int i, const DSTable &, int *);
 
    /** Uniform Refinement. Element with index i is refined uniformly. */
    void UniformRefinement(int i, const DSTable &, int *, int *, int *);

mesh/pmesh.cpp

@@ -62,17 +62,9 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
    // If pmesh has a ParNURBSExtension, it was copied by the Mesh copy ctor, so
    // there is no need to do anything here.
 
-   // Copy ParNCMesh, if present
-   if (pmesh.pncmesh)
-   {
-      pncmesh = new ParNCMesh(*pmesh.pncmesh);
-      pncmesh->OnMeshUpdated(this);
-   }
-   else
-   {
-      pncmesh = NULL;
-   }
-   ncmesh = pncmesh;
+   MFEM_VERIFY(pmesh.pncmesh == NULL,
+               "copy of parallel non-conforming meshes is not implemented");
+   pncmesh = NULL;
 
    // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
    // and the FiniteElementSpace (as a ParFiniteElementSpace)
@@ -1222,76 +1214,48 @@ int ParMesh::GetEdgeSplittings(Element *edge, const DSTable &v_to_v,
    }
 }
 
-void ParMesh::GetFaceSplittings(Element *face, const HashTable<Hashed2> &v_to_v,
-                                Array<unsigned> &codes)
+// For a triangular face with (correctly ordered) vertices v[0], v[1], v[2]
+// return a number with the following meaning:
+// 0 - the face was not refined
+// 1 - the face was refined once by splitting v[0],v[1]
+// 2 - the face was refined twice by splitting v[0],v[1] and then v[1],v[2]
+// 3 - the face was refined twice by splitting v[0],v[1] and then v[0],v[2]
+// 4 - the face was refined three times (as in 2+3)
+int ParMesh::GetFaceSplittings(Element *face, const DSTable &v_to_v,
+                               int *middle)
 {
-   const int *v = face->GetVertices();
-   typedef Triple<int,int,int> face_t;
-   Array<face_t> face_stack;
+   int m, right = 0;
+   int number_of_splittings = 0;
+   int *v = face->GetVertices();
 
-   unsigned code = 0;
-   face_stack.Append(face_t(v[0], v[1], v[2]));
-   for (unsigned bit = 0; face_stack.Size() > 0; bit++)
+   if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
    {
-      if (bit == 8*sizeof(unsigned))
+      number_of_splittings++;
+      if ((m = v_to_v(v[1], v[2])) != -1 && middle[m] != -1)
       {
-         codes.Append(code);
-         code = bit = 0;
+         right = 1;
+         number_of_splittings++;
+      }
+      if ((m = v_to_v(v[2], v[0])) != -1 && middle[m] != -1)
+      {
+         number_of_splittings++;
       }
 
-      const face_t &f = face_stack.Last();
-      int mid = v_to_v.FindId(f.one, f.two);
-      if (mid == -1)
+      switch (number_of_splittings)
       {
-         // leave a 0 at bit 'bit'
-         face_stack.DeleteLast();
-      }
-      else
-      {
-         code += (1 << bit); // set bit 'bit' to 1
-         mid += NumOfVertices;
-         face_stack.Append(face_t(f.three, f.one, mid));
-         face_t &r = face_stack[face_stack.Size()-2];
-         r = face_t(r.two, r.three, mid);
+         case 2:
+            if (right == 0)
+            {
+               number_of_splittings++;
+            }
+            break;
+         case 3:
+            number_of_splittings++;
+            break;
       }
    }
-   codes.Append(code);
-}
 
-bool ParMesh::DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
-                                   const Array<unsigned> &codes, int &pos)
-{
-   typedef Triple<int,int,int> face_t;
-   Array<face_t> face_stack;
-
-   bool need_refinement = 0;
-   face_stack.Append(face_t(v[0], v[1], v[2]));
-   for (unsigned bit = 0, code = codes[pos++]; face_stack.Size() > 0; bit++)
-   {
-      if (bit == 8*sizeof(unsigned))
-      {
-         code = codes[pos++];
-         bit = 0;
-      }
-
-      if ((code & (1 << bit)) == 0) { face_stack.DeleteLast(); continue; }
-
-      const face_t &f = face_stack.Last();
-      int mid = v_to_v.FindId(f.one, f.two);
-      if (mid == -1)
-      {
-         mid = v_to_v.GetId(f.one, f.two);
-         int ind[2] = { f.one, f.two };
-         vertices.Append(Vertex());
-         AverageVertices(ind, 2, vertices.Size()-1);
-         need_refinement = 1;
-      }
-      mid += NumOfVertices;
-      face_stack.Append(face_t(f.three, f.one, mid));
-      face_t &r = face_stack[face_stack.Size()-2];
-      r = face_t(r.two, r.three, mid);
-   }
-   return need_refinement;
+   return number_of_splittings;
 }
 
 void ParMesh::GenerateOffsets(int N, HYPRE_Int loc_sizes[],
@@ -2149,6 +2113,8 @@ void ParMesh::ReorientTetMesh()
 
 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 {
+   int i, j;
+
    if (pncmesh)
    {
       MFEM_ABORT("Local and nonconforming refinements cannot be mixed.");
@@ -2167,63 +2133,76 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
          uniform_refinement = 1;
       }
 
-      // 1. Hash table of vertex to vertex connections corresponding to refined
-      //    edges.
-      HashTable<Hashed2> v_to_v;
+      // 1. Get table of vertex to vertex connections.
+      DSTable v_to_v(NumOfVertices);
+      GetVertexToVertexTable(v_to_v);
 
-      // 2. Do the red refinement.
+      // 2. Create a marker array for all edges (vertex to vertex connections).
+      Array<int> middle(v_to_v.NumberOfEntries());
+      middle = -1;
+
+      // 3. Do the red refinement.
       switch (type)
       {
          case 1:
-            for (int i = 0; i < marked_el.Size(); i++)
+            for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
          case 2:
-            for (int i = 0; i < marked_el.Size(); i++)
+            for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
 
-               Bisection(NumOfElements - 1, v_to_v);
-               Bisection(marked_el[i], v_to_v);
+               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
          case 3:
-            for (int i = 0; i < marked_el.Size(); i++)
+            for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
 
-               int j = NumOfElements - 1;
-               Bisection(j, v_to_v);
-               Bisection(NumOfElements - 1, v_to_v);
-               Bisection(j, v_to_v);
+               j = NumOfElements - 1;
+               Bisection(j, v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
+               Bisection(j, v_to_v, NULL, NULL, middle);
 
-               Bisection(marked_el[i], v_to_v);
-               Bisection(NumOfElements-1, v_to_v);
-               Bisection(marked_el[i], v_to_v);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
             }
             break;
       }
 
-      // 3. Do the green refinement (to get conforming mesh).
+      // 4. Do the green refinement (to get conforming mesh).
       int need_refinement;
-      int max_faces_in_group = 0;
-      // face_splittings identify how the shared faces have been split
-      Array<unsigned> *face_splittings = new Array<unsigned>[GetNGroups()-1];
-      for (int i = 0; i < GetNGroups()-1; i++)
+      int refined_edge[5][3] =
       {
-         const int faces_in_group = GroupNFaces(i+1);
-         face_splittings[i].Reserve(faces_in_group);
+         {0, 0, 0},
+         {1, 0, 0},
+         {1, 1, 0},
+         {1, 0, 1},
+         {1, 1, 1}
+      };
+      int faces_in_group, max_faces_in_group = 0;
+      // face_splittings identify how the shared faces have been split
+      int **face_splittings = new int*[GetNGroups()-1];
+      for (i = 0; i < GetNGroups()-1; i++)
+      {
+         faces_in_group = GroupNFaces(i+1);
+         face_splittings[i] = new int[faces_in_group];
          if (faces_in_group > max_faces_in_group)
          {
             max_faces_in_group = faces_in_group;
          }
       }
-      int neighbor;
-      Array<unsigned> iBuf(max_faces_in_group);
+      int neighbor, *iBuf = new int[max_faces_in_group];
 
-      MPI_Request *requests = new MPI_Request[GetNGroups()-1];
+      Array<int> group_faces;
+
+      MPI_Request request;
       MPI_Status  status;
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
@@ -2232,12 +2211,12 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       do
       {
          need_refinement = 0;
-         for (int i = 0; i < NumOfElements; i++)
+         for (i = 0; i < NumOfElements; i++)
          {
-            if (elements[i]->NeedRefinement(v_to_v))
+            if (elements[i]->NeedRefinement(v_to_v, middle))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v);
+               Bisection(i, v_to_v, NULL, NULL, middle);
             }
          }
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
@@ -2260,85 +2239,96 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             const int tag = 293;
 
             // (a) send the type of interface splitting
-            int req_count = 0;
-            for (int i = 0; i < GetNGroups()-1; i++)
+            for (i = 0; i < GetNGroups()-1; i++)
             {
-               const int *group_faces = group_sface.GetRow(i);
-               const int faces_in_group = group_sface.RowSize(i);
+               group_sface.GetRow(i, group_faces);
+               faces_in_group = group_faces.Size();
                // it is enough to communicate through the faces
                if (faces_in_group == 0) { continue; }
 
-               face_splittings[i].SetSize(0);
-               for (int j = 0; j < faces_in_group; j++)
+               for (j = 0; j < faces_in_group; j++)
                {
-                  GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
-                                    face_splittings[i]);
+                  face_splittings[i][j] =
+                     GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
+                                       middle);
                }
                const int *nbs = gtopo.GetGroup(i+1);
                neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
-               MPI_Isend(face_splittings[i], face_splittings[i].Size(),
-                         MPI_UNSIGNED, neighbor, tag, MyComm,
-                         &requests[req_count++]);
+               MPI_Isend(face_splittings[i], faces_in_group, MPI_INT,
+                         neighbor, tag, MyComm, &request);
             }
 
             // (b) receive the type of interface splitting
-            for (int i = 0; i < GetNGroups()-1; i++)
+            for (i = 0; i < GetNGroups()-1; i++)
             {
-               const int *group_faces = group_sface.GetRow(i);
-               const int faces_in_group = group_sface.RowSize(i);
+               group_sface.GetRow(i, group_faces);
+               faces_in_group = group_faces.Size();
                if (faces_in_group == 0) { continue; }
 
                const int *nbs = gtopo.GetGroup(i+1);
                neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
-               MPI_Probe(neighbor, tag, MyComm, &status);
-               int count;
-               MPI_Get_count(&status, MPI_UNSIGNED, &count);
-               iBuf.SetSize(count);
-               MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag, MyComm,
-                        MPI_STATUS_IGNORE);
+               MPI_Recv(iBuf, faces_in_group, MPI_INT, neighbor,
+                        tag, MyComm, &status);
 
-               for (int j = 0, pos = 0; j < faces_in_group; j++)
+               for (j = 0; j < faces_in_group; j++)
                {
-                  const int *v = shared_faces[group_faces[j]]->GetVertices();
-                  need_refinement |= DecodeFaceSplittings(v_to_v, v, iBuf, pos);
+                  if (iBuf[j] == face_splittings[i][j]) { continue; }
+
+                  int *v = shared_faces[group_faces[j]]->GetVertices();
+                  for (int k = 0; k < 3; k++)
+                  {
+                     if (refined_edge[iBuf[j]][k] != 1 ||
+                         refined_edge[face_splittings[i][j]][k] != 0)
+                     { continue; }
+
+                     int ind[2] = { v[k], v[(k+1)%3] };
+                     int ii = v_to_v(ind[0], ind[1]);
+                     if (middle[ii] == -1)
+                     {
+                        need_refinement = 1;
+                        middle[ii] = NumOfVertices++;
+                        vertices.Append(Vertex());
+                        AverageVertices(ind, 2, vertices.Size()-1);
+                     }
+                  }
                }
             }
 
-            int nr = need_refinement;
-            MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
-
-            MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
+            i = need_refinement;
+            MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
          }
       }
       while (need_refinement == 1);
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
+      i = ref_loops_all;
+      MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
+      if (MyRank == 0)
       {
-         int i = ref_loops_all;
-         MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
-         if (MyRank == 0)
-         {
-            mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
-                      << ref_loops_all << ", ref_loops_par = " << ref_loops_par
-                      << '\n' << endl;
-         }
+         mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
+                   << ref_loops_all << ", ref_loops_par = " << ref_loops_par
+                   << '\n' << endl;
       }
 #endif
 
-      delete [] requests;
-      iBuf.DeleteAll();
+      delete [] iBuf;
+      for (i = 0; i < GetNGroups()-1; i++)
+      {
+         delete [] face_splittings[i];
+      }
       delete [] face_splittings;
 
-      // 4. Update the boundary elements.
+
+      // 5. Update the boundary elements.
       do
       {
          need_refinement = 0;
-         for (int i = 0; i < NumOfBdrElements; i++)
+         for (i = 0; i < NumOfBdrElements; i++)
          {
-            if (boundary[i]->NeedRefinement(v_to_v))
+            if (boundary[i]->NeedRefinement(v_to_v, middle))
             {
                need_refinement = 1;
-               BdrBisection(i, v_to_v);
+               Bisection(i, v_to_v, middle);
             }
          }
       }
@@ -2352,16 +2342,31 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       DeleteLazyTables();
 
-      // 5. Update the groups after refinement.
+      // 5a. Update the groups after refinement.
       if (el_to_face != NULL)
       {
-         RefineGroups(v_to_v);
+         RefineGroups(v_to_v, middle);
          // GetElementToFaceTable(); // Called by RefineGroups
          GenerateFaces();
       }
-      NumOfVertices = vertices.Size();
 
-      // 6. Update element-to-edge relations.
+      // 6. Un-mark the Pf elements.
+      int refinement_edges[2], type, flag;
+      for (i = 0; i < NumOfElements; i++)
+      {
+         Tetrahedron* el = (Tetrahedron*) elements[i];
+         el->ParseRefinementFlag(refinement_edges, type, flag);
+
+         if (type == Tetrahedron::TYPE_PF)
+         {
+            el->CreateRefinementFlag(refinement_edges, Tetrahedron::TYPE_PU,
+                                     flag);
+         }
+      }
+
+      // 7. Free the allocated memory.
+      middle.DeleteAll();
+
       if (el_to_edge != NULL)
       {
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
@@ -2388,15 +2393,15 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       int *edge2  = new int[nedges];
       int *middle = new int[nedges];
 
-      for (int i = 0; i < nedges; i++)
+      for (i = 0; i < nedges; i++)
       {
          edge1[i] = edge2[i] = middle[i] = -1;
       }
 
-      for (int i = 0; i < NumOfElements; i++)
+      for (i = 0; i < NumOfElements; i++)
       {
          int *v = elements[i]->GetVertices();
-         for (int j = 0; j < 3; j++)
+         for (j = 0; j < 3; j++)
          {
             int ind = v_to_v(v[j], v[(j+1)%3]);
             (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
@@ -2404,7 +2409,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       }
 
       // 3. Do the red refinement.
-      for (int i = 0; i < marked_el.Size(); i++)
+      for (i = 0; i < marked_el.Size(); i++)
       {
          RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
       }
@@ -2414,7 +2419,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       int edges_in_group, max_edges_in_group = 0;
       // edge_splittings identify how the shared edges have been split
       int **edge_splittings = new int*[GetNGroups()-1];
-      for (int i = 0; i < GetNGroups()-1; i++)
+      for (i = 0; i < GetNGroups()-1; i++)
       {
          edges_in_group = GroupNEdges(i+1);
          edge_splittings[i] = new int[edges_in_group];
@@ -2438,14 +2443,12 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       do
       {
          need_refinement = 0;
-         for (int i = 0; i < nedges; i++)
-         {
+         for (i = 0; i < nedges; i++)
             if (middle[i] != -1 && edge1[i] != -1)
             {
                need_refinement = 1;
                GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
             }
-         }
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
          ref_loops_all++;
 #endif
@@ -2465,14 +2468,14 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             // MPI_Barrier(MyComm);
 
             // (a) send the type of interface splitting
-            for (int i = 0; i < GetNGroups()-1; i++)
+            for (i = 0; i < GetNGroups()-1; i++)
             {
                group_sedge.GetRow(i, group_edges);
                edges_in_group = group_edges.Size();
                // it is enough to communicate through the edges
                if (edges_in_group != 0)
                {
-                  for (int j = 0; j < edges_in_group; j++)
+                  for (j = 0; j < edges_in_group; j++)
                   {
                      edge_splittings[i][j] =
                         GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
@@ -2493,7 +2496,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             }
 
             // (b) receive the type of interface splitting
-            for (int i = 0; i < GetNGroups()-1; i++)
+            for (i = 0; i < GetNGroups()-1; i++)
             {
                group_sedge.GetRow(i, group_edges);
                edges_in_group = group_edges.Size();
@@ -2511,18 +2514,15 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                   MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
                            MPI_ANY_TAG, MyComm, &status);
 
-                  for (int j = 0; j < edges_in_group; j++)
-                  {
+                  for (j = 0; j < edges_in_group; j++)
                      if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
                      {
                         int *v = shared_edges[group_edges[j]]->GetVertices();
                         int ii = v_to_v(v[0], v[1]);
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
                         if (middle[ii] != -1)
-                        {
                            mfem_error("ParMesh::LocalRefinement (triangles) : "
                                       "Oops!");
-                        }
 #endif
                         need_refinement = 1;
                         middle[ii] = NumOfVertices++;
@@ -2532,30 +2532,27 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                         }
                         vertices.Append(V);
                      }
-                  }
                }
             }
 
-            int nr = need_refinement;
-            MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
+            i = need_refinement;
+            MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
          }
       }
       while (need_refinement == 1);
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
+      i = ref_loops_all;
+      MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
+      if (MyRank == 0)
       {
-         int i = ref_loops_all;
-         MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
-         if (MyRank == 0)
-         {
-            mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
-                      << ref_loops_all << ", ref_loops_par = " << ref_loops_par
-                      << '\n' << endl;
-         }
+         mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
+                   << ref_loops_all << ", ref_loops_par = " << ref_loops_par
+                   << '\n' << endl;
       }
 #endif
 
-      for (int i = 0; i < GetNGroups()-1; i++)
+      for (i = 0; i < GetNGroups()-1; i++)
       {
          delete [] edge_splittings[i];
       }
@@ -2566,7 +2563,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 5. Update the boundary elements.
       int v1[2], v2[2], bisect, temp;
       temp = NumOfBdrElements;
-      for (int i = 0; i < temp; i++)
+      for (i = 0; i < temp; i++)
       {
          int *v = boundary[i]->GetVertices();
          bisect = v_to_v(v[0], v[1]);
@@ -2582,10 +2579,8 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
             }
             else
-            {
                mfem_error("Only bisection of segment is implemented for bdr"
                           " elem.");
-            }
          }
       }
       NumOfBdrElements = boundary.Size();
@@ -2616,9 +2611,9 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       elements.SetSize(NumOfElements);
       CoarseFineTr.embeddings.SetSize(NumOfElements);
 
-      for (int j = 0; j < marked_el.Size(); j++)
+      for (j = 0; j < marked_el.Size(); j++)
       {
-         int i = marked_el[j];
+         i = marked_el[j];
          Segment *c_seg = (Segment *)elements[i];
          int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
          int new_v = cnv + j, new_e = cne + j;
@@ -2967,226 +2962,6 @@ void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
    }
 }
 
-void ParMesh::RefineGroups(const HashTable<Hashed2> &v_to_v)
-{
-   int i, attr, ind, *v;
-
-   int group;
-   Array<int> group_verts, group_edges, group_faces;
-
-   // To update the groups after a refinement, we observe that:
-   // - every (new and old) vertex, edge and face belongs to exactly one group
-   // - the refinement does not create new groups
-   // - a new vertex appears only as the middle of a refined edge
-   // - a face can be refined multiple times producing new edges and faces
-
-   Array<Segment *> sedge_stack;
-   Array<Triangle *> sface_stack;
-
-   Array<int> I_group_svert, J_group_svert;
-   Array<int> I_group_sedge, J_group_sedge;
-   Array<int> I_group_sface, J_group_sface;
-
-   I_group_svert.SetSize(GetNGroups()+1);
-   I_group_sedge.SetSize(GetNGroups()+1);
-   if (Dim == 3)
-   {
-      I_group_sface.SetSize(GetNGroups()+1);
-   }
-
-   I_group_svert[0] = I_group_svert[1] = 0;
-   I_group_sedge[0] = I_group_sedge[1] = 0;
-   if (Dim == 3)
-   {
-      I_group_sface[0] = I_group_sface[1] = 0;
-   }
-
-   for (group = 0; group < GetNGroups()-1; group++)
-   {
-      // Get the group shared objects
-      group_svert.GetRow(group, group_verts);
-      group_sedge.GetRow(group, group_edges);
-      group_sface.GetRow(group, group_faces);
-
-      // Check which edges have been refined
-      for (i = 0; i < group_sedge.RowSize(group); i++)
-      {
-         v = shared_edges[group_edges[i]]->GetVertices();
-         ind = v_to_v.FindId(v[0], v[1]);
-         if (ind == -1) { continue; }
-
-         // This shared edge is refined: walk the whole refinement tree
-         attr = shared_edges[group_edges[i]]->GetAttribute();
-         do
-         {
-            ind += NumOfVertices;
-            // Add new shared vertex
-            group_verts.Append(svert_lvert.Append(ind)-1);
-            // Put the right sub-edge on top of the stack
-            sedge_stack.Append(new Segment(ind, v[1], attr));
-            // The left sub-edge replaces the original edge
-            v[1] = ind;
-            ind = v_to_v.FindId(v[0], ind);
-         }
-         while (ind != -1);
-         // Process all edges in the edge stack
-         do
-         {
-            Segment *se = sedge_stack.Last();
-            v = se->GetVertices();
-            ind = v_to_v.FindId(v[0], v[1]);
-            if (ind == -1)
-            {
-               // The edge 'se' is not refined
-               sedge_stack.DeleteLast();
-               // Add new shared edge
-               shared_edges.Append(se);
-               group_edges.Append(sedge_ledge.Append(-1)-1);
-            }
-            else
-            {
-               // The edge 'se' is refined
-               ind += NumOfVertices;
-               // Add new shared vertex
-               group_verts.Append(svert_lvert.Append(ind)-1);
-               // Put the left sub-edge on top of the stack
-               sedge_stack.Append(new Segment(v[0], ind, attr));
-               // The right sub-edge replaces the original edge
-               v[0] = ind;
-            }
-         }
-         while (sedge_stack.Size() > 0);
-      }
-
-      // Check which faces have been refined
-      for (i = 0; i < group_sface.RowSize(group); i++)
-      {
-         v = shared_faces[group_faces[i]]->GetVertices();
-         ind = v_to_v.FindId(v[0], v[1]);
-         if (ind == -1) { continue; }
-
-         // This shared face is refined: walk the whole refinement tree
-         attr = shared_faces[group_faces[i]]->GetAttribute();
-         const int edge_attr = 1;
-         do
-         {
-            ind += NumOfVertices;
-            // Add the refinement edge to the edge stack
-            sedge_stack.Append(new Segment(v[2], ind, edge_attr));
-            // Put the right sub-triangle on top of the face stack
-            sface_stack.Append(new Triangle(v[1], v[2], ind, attr));
-            // The left sub-triangle replaces the original one
-            v[1] = v[0]; v[0] = v[2]; v[2] = ind;
-            ind = v_to_v.FindId(v[0], v[1]);
-         }
-         while (ind != -1);
-         // Process all faces (triangles) in the face stack
-         do
-         {
-            Triangle *st = sface_stack.Last();
-            v = st->GetVertices();
-            ind = v_to_v.FindId(v[0], v[1]);
-            if (ind == -1)
-            {
-               // The triangle 'st' is not refined
-               sface_stack.DeleteLast();
-               // Add new shared face
-               shared_faces.Append(st);
-               group_faces.Append(sface_lface.Append(-1)-1);
-            }
-            else
-            {
-               // The triangle 'st' is refined
-               ind += NumOfVertices;
-               // Add the refinement edge to the edge stack
-               sedge_stack.Append(new Segment(v[2], ind, edge_attr));
-               // Put the left sub-triangle on top of the face stack
-               sface_stack.Append(new Triangle(v[2], v[0], ind, attr));
-               // The right sub-triangle replaces the original one
-               v[0] = v[1]; v[1] = v[2]; v[2] = ind;
-            }
-         }
-         while (sface_stack.Size() > 0);
-         // Process all edges in the edge stack (same code as above)
-         do
-         {
-            Segment *se = sedge_stack.Last();
-            v = se->GetVertices();
-            ind = v_to_v.FindId(v[0], v[1]);
-            if (ind == -1)
-            {
-               // The edge 'se' is not refined
-               sedge_stack.DeleteLast();
-               // Add new shared edge
-               shared_edges.Append(se);
-               group_edges.Append(sedge_ledge.Append(-1)-1);
-            }
-            else
-            {
-               // The edge 'se' is refined
-               ind += NumOfVertices;
-               // Add new shared vertex
-               group_verts.Append(svert_lvert.Append(ind)-1);
-               // Put the left sub-edge on top of the stack
-               sedge_stack.Append(new Segment(v[0], ind, attr));
-               // The right sub-edge replaces the original edge
-               v[0] = ind;
-            }
-         }
-         while (sedge_stack.Size() > 0);
-      }
-
-      I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
-      I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
-      if (Dim == 3)
-      {
-         I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
-      }
-
-      J_group_svert.Append(group_verts);
-      J_group_sedge.Append(group_edges);
-      if (Dim == 3)
-      {
-         J_group_sface.Append(group_faces);
-      }
-   }
-
-   // Fix the local numbers of shared edges and faces: sedge_ledge, sface_lface
-   {
-      NumOfVertices = vertices.Size();
-      DSTable new_v_to_v(NumOfVertices);
-      GetVertexToVertexTable(new_v_to_v);
-      for (i = 0; i < shared_edges.Size(); i++)
-      {
-         v = shared_edges[i]->GetVertices();
-         sedge_ledge[i] = new_v_to_v(v[0], v[1]);
-      }
-   }
-   if (Dim == 3)
-   {
-      STable3D *faces_tbl = GetElementToFaceTable(1);
-      for (i = 0; i < shared_faces.Size(); i++)
-      {
-         v = shared_faces[i]->GetVertices();
-         sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
-      }
-      delete faces_tbl;
-   }
-
-   group_svert.SetIJ(I_group_svert, J_group_svert);
-   group_sedge.SetIJ(I_group_sedge, J_group_sedge);
-   if (Dim == 3)
-   {
-      group_sface.SetIJ(I_group_sface, J_group_sface);
-   }
-   I_group_svert.LoseData(); J_group_svert.LoseData();
-   I_group_sedge.LoseData(); J_group_sedge.LoseData();
-   if (Dim == 3)
-   {
-      I_group_sface.LoseData(); J_group_sface.LoseData();
-   }
-}
-
 void ParMesh::QuadUniformRefinement()
 {
    DeleteFaceNbrData();

mesh/pmesh.hpp

@@ -63,12 +63,8 @@ protected:
 
    /// Return a number(0-1) identifying how the given edge has been split
    int GetEdgeSplittings(Element *edge, const DSTable &v_to_v, int *middle);
-   /// Append codes identifying how the given face has been split to @a codes
-   void GetFaceSplittings(Element *face, const HashTable<Hashed2> &v_to_v,
-                          Array<unsigned> &codes);
-
-   bool DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
-                             const Array<unsigned> &codes, int &pos);
+   /// Return a number(0-4) identifying how the given face has been split
+   int GetFaceSplittings(Element *face, const DSTable &v_to_v, int *middle);
 
    void GetFaceNbrElementTransformation(
       int i, IsoparametricTransformation *ElTr);
@@ -189,10 +185,8 @@ public:
    /// Utility function: sum integers from all processors (Allreduce).
    virtual long ReduceInt(int value) const;
 
-   /// Update the groups after triangle refinement
+   /// Update the groups after tet refinement
    void RefineGroups(const DSTable &v_to_v, int *middle);
-   /// Update the groups after tetrahedron refinement
-   void RefineGroups(const HashTable<Hashed2> &v_to_v);
 
    /// Load balance the mesh. NC meshes only.
    void Rebalance();

mesh/pncmesh.cpp

@@ -46,16 +46,6 @@ ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh)
    // branches that only contain someone else's leaves (see Prune())
 }
 
-ParNCMesh::ParNCMesh(const ParNCMesh &other)
-// copy primary data only
-   : NCMesh(other)
-   , MyComm(other.MyComm)
-   , NRanks(other.NRanks)
-   , MyRank(other.MyRank)
-{
-   Update(); // mark all secondary stuff for recalculation
-}
-
 ParNCMesh::~ParNCMesh()
 {
    ClearAuxPM();

mesh/pncmesh.hpp

@@ -63,8 +63,6 @@ class ParNCMesh : public NCMesh
 public:
    ParNCMesh(MPI_Comm comm, const NCMesh& ncmesh);
 
-   ParNCMesh(const ParNCMesh &other);
-
    virtual ~ParNCMesh();
 
    /** An override of NCMesh::Refine, which is called eventually, after making

mesh/tetrahedron.cpp

@@ -153,14 +153,22 @@ void Tetrahedron::GetMarkedFace(const int face, int *fv)
    }
 }
 
-int Tetrahedron::NeedRefinement(HashTable<Hashed2> &v_to_v) const
+int Tetrahedron::NeedRefinement(DSTable &v_to_v, int *middle) const
 {
-   if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[0], indices[3]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[1], indices[3]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[2], indices[3]) != -1) { return 1; }
+   int m;
+
+   if ((m = v_to_v(indices[0], indices[1])) != -1)
+      if (middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[1], indices[2])) != -1)
+      if (middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[2], indices[0])) != -1)
+      if (middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[0], indices[3])) != -1)
+      if (middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[1], indices[3])) != -1)
+      if (middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[2], indices[3])) != -1)
+      if (middle[m] != -1) { return 1; }
    return 0;
 }
 

mesh/tetrahedron.hpp

@@ -64,7 +64,7 @@ public:
    void SetRefinementFlag(int rf) { refinement_flag = rf; }
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const;
+   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const;
 
    /// Set the vertices according to the given input.
    virtual void SetVertices(const int *ind);

mesh/triangle.cpp

@@ -34,11 +34,13 @@ Triangle::Triangle(int ind1, int ind2, int ind3, int attr)
    transform = 0;
 }
 
-int Triangle::NeedRefinement(HashTable<Hashed2> &v_to_v) const
+int Triangle::NeedRefinement(DSTable &v_to_v, int *middle) const
 {
-   if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
-   if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
+   int m;
+
+   if ((m = v_to_v(indices[0], indices[1])) != -1 && middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[1], indices[2])) != -1 && middle[m] != -1) { return 1; }
+   if ((m = v_to_v(indices[2], indices[0])) != -1 && middle[m] != -1) { return 1; }
    return 0;
 }
 

mesh/triangle.hpp

@@ -42,7 +42,7 @@ public:
    virtual int GetType() const { return Element::TRIANGLE; }
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const;
+   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const;
 
    /// Set the vertices according to the given input.
    virtual void SetVertices(const int *ind);

miniapps/electromagnetics/makefile

@@ -51,22 +51,6 @@ $(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
 %.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
 
-# Rules to copy the *.mesh files - needed for running the sample runs when
-# building out-of-source:
-ifneq ($(SRC),)
-JOULE_MESH_FILES = cylinder-hex.mesh cylinder-tet.mesh
-ifeq ($(MFEM_USE_NETCDF),YES)
-   JOULE_MESH_FILES += cylinder-hex-q2.gen cylinder-tet-p2.gen
-endif
-TESLA_MESH_FILES = square-angled-pipe.mesh
-VOLTA_MESH_FILES = llnl.mesh
-$(JOULE_MESH_FILES) $(TESLA_MESH_FILES) $(VOLTA_MESH_FILES): %: $(SRC)%
-	ln -sf $(<) .
-joule: | $(JOULE_MESH_FILES)
-tesla: | $(TESLA_MESH_FILES)
-volta: | $(VOLTA_MESH_FILES)
-endif
-
 MFEM_TESTS = MINIAPPS
 include $(MFEM_TEST_MK)
 
@@ -83,7 +67,7 @@ maxwell-test-par: maxwell
 	-abcs '-1' -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5')
 joule-test-par: joule
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-m cylinder-hex.mesh -p rod -tf 3)
+	-p rod -tf 3 -m $(SRC)cylinder-hex.mesh)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/meshing/makefile

@@ -42,15 +42,6 @@ endif
 
 all: $(MINIAPPS)
 
-# Rules to copy the *.mesh files - needed for running the sample runs when
-# building out-of-source:
-ifneq ($(SRC),)
-MESH_FILES = blade.mesh icf.mesh
-$(MESH_FILES): %: $(SRC)%
-	ln -sf $(<) .
-mesh-optimizer pmesh-optimizer: | $(MESH_FILES)
-endif
-
 MFEM_TESTS = MINIAPPS
 include $(MFEM_TEST_MK)
 
@@ -61,9 +52,10 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 %-test-seq: %
 	@$(call mfem-test-file,$<,, Meshing miniapp,$(<).mesh)
 mesh-optimizer-test-seq: mesh-optimizer
-	@$(call mfem-test,$<,, Meshing miniapp)
+	@$(call mfem-test,$<,, Meshing miniapp,-m $(SRC)icf.mesh)
 pmesh-optimizer-test-par: pmesh-optimizer
-	@$(call mfem-test,$<, $(RUN_MPI), Parallel meshing miniapp)
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel meshing miniapp,\
+	-m $(SRC)icf.mesh)
 
 # Testing: Specific execution options
 mesh-explorer-test-seq: