Merge pull request #1 from camierjs/minimal-surf-dev-rewritten

Minimal surf dev - Costa

.gitignore

@@ -161,6 +161,8 @@ miniapps/meshing/shaper
 miniapps/meshing/extruder
 miniapps/meshing/mesh-optimizer
 miniapps/meshing/pmesh-optimizer
+miniapps/meshing/mesh-minimal-surface
+miniapps/meshing/pmesh-minimal-surface
 
 miniapps/meshing/mobius-strip.mesh
 miniapps/meshing/klein-bottle.mesh

fem/bilininteg.cpp

@@ -2135,17 +2135,17 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
    ElementTransformation &Trans,
    DenseMatrix &elmat)
 {
-   int dim = el.GetDim();
-   int dof = el.GetDof();
-
-   double norm;
+   const int dim = el.GetDim();
+   const int dof = el.GetDof();
+   const int sdim = Trans.GetSpaceDim();
+   const bool square = (dim == sdim);
+   double w;
 
-   elmat.SetSize (dim * dof);
+   elmat.SetSize(sdim * dof);
 
-   Jinv.  SetSize (dim);
-   dshape.SetSize (dof, dim);
-   gshape.SetSize (dof, dim);
-   pelmat.SetSize (dof);
+   dshape.SetSize(dof, dim);
+   dshapedxt.SetSize(dof, sdim);
+   pelmat.SetSize(dof);
 
    const IntegrationRule *ir = IntRule;
    if (ir == NULL)
@@ -2162,36 +2162,29 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
       }
    }
 
-   elmat = 0.0;
+   pelmat = 0.0;
 
    for (int i = 0; i < ir -> GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
-
       el.CalcDShape (ip, dshape);
-
       Trans.SetIntPoint (&ip);
-      norm = ip.weight * Trans.Weight();
-      CalcInverse (Trans.Jacobian(), Jinv);
-
-      Mult (dshape, Jinv, gshape);
-
-      MultAAt (gshape, pelmat);
-
-      if (Q)
-      {
-         norm *= Q -> Eval (Trans, ip);
-      }
-
-      pelmat *= norm;
-
-      for (int d = 0; d < dim; d++)
+      w = Trans.Weight();
+      w = ip.weight / (square ? w : w*w*w);
+      // AdjugateJacobian = / adj(J),         if J is square
+      //                    \ adj(J^t.J).J^t, otherwise
+      Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
+      if (Q) {  w *= Q -> Eval (Trans, ip); }
+      AddMult_a_AAt(w, dshapedxt, pelmat);
+   }
+   for (int d = 0; d < sdim; d++)
+   {
+      for (int k = 0; k < dof; k++)
       {
-         for (int k = 0; k < dof; k++)
-            for (int l = 0; l < dof; l++)
-            {
-               elmat (dof*d+k, dof*d+l) += pelmat (k, l);
-            }
+         for (int l = 0; l < dof; l++)
+         {
+            elmat(dof*d+k, dof*d+l) = pelmat(k, l);
+         }
       }
    }
 }

fem/bilininteg.hpp

@@ -2391,14 +2391,12 @@ class VectorDiffusionIntegrator : public BilinearFormIntegrator
    // PA extension
    const DofToQuad *maps;         ///< Not owned
    const GeometricFactors *geom;  ///< Not owned
-   int dim, ne, dofs1D, quad1D;
+   int dim, sdim, ne, dofs1D, quad1D;
    Vector pa_data;
 
 private:
-   DenseMatrix Jinv;
-   DenseMatrix dshape;
-   DenseMatrix gshape;
-   DenseMatrix pelmat;
+   DenseMatrix dshape, dshapedxt, pelmat;
+   DenseMatrix Jinv, gshape;
 
 public:
    VectorDiffusionIntegrator() { Q = NULL; }

fem/bilininteg_diffusion.cpp

@@ -81,12 +81,20 @@ static void OccaPADiffusionSetup3D(const int D1D,
 #endif // MFEM_USE_OCCA
 
 // PA Diffusion Assemble 2D kernel
+template<const int T_SDIM>
 static void PADiffusionSetup2D(const int Q1D,
                                const int NE,
                                const Array<double> &w,
                                const Vector &j,
                                const Vector &c,
-                               Vector &d)
+                               Vector &d);
+template<>
+void PADiffusionSetup2D<2>(const int Q1D,
+                           const int NE,
+                           const Array<double> &w,
+                           const Vector &j,
+                           const Vector &c,
+                           Vector &d)
 {
    const int NQ = Q1D*Q1D;
    const bool const_c = c.Size() == 1;
@@ -112,6 +120,48 @@ static void PADiffusionSetup2D(const int Q1D,
    });
 }
 
+// PA Diffusion Assemble 2D kernel with 3D node coords
+template<>
+void PADiffusionSetup2D<3>(const int Q1D,
+                           const int NE,
+                           const Array<double> &w,
+                           const Vector &j,
+                           const Vector &c,
+                           Vector &d)
+{
+   constexpr int DIM = 2;
+   constexpr int SDIM = 3;
+   const int NQ = Q1D*Q1D;
+   const bool const_c = c.Size() == 1;
+
+   auto W = w.Read();
+   auto J = Reshape(j.Read(), NQ, SDIM, DIM, NE);
+   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
+   auto D = Reshape(d.Write(), NQ, 3, NE);
+   MFEM_FORALL(e, NE,
+   {
+      for (int q = 0; q < NQ; ++q)
+      {
+         const double wq = W[q];
+         const double J11 = J(q,0,0,e);
+         const double J21 = J(q,1,0,e);
+         const double J31 = J(q,2,0,e);
+         const double J12 = J(q,0,1,e);
+         const double J22 = J(q,1,1,e);
+         const double J32 = J(q,2,1,e);
+         const double E = J11*J11 + J21*J21 + J31*J31;
+         const double G = J12*J12 + J22*J22 + J32*J32;
+         const double F = J11*J12 + J21*J22 + J31*J32;
+         const double iw = 1.0 / sqrt(E*G - F*F);
+         const double coeff = const_c ? C(0,0) : C(q,e);
+         const double alpha = wq * coeff * iw;
+         D(q,0,e) =  alpha * G; // 1,1
+         D(q,1,e) = -alpha * F; // 1,2
+         D(q,2,e) =  alpha * E; // 2,2
+      }
+   });
+}
+
 // PA Diffusion Assemble 3D kernel
 static void PADiffusionSetup3D(const int Q1D,
                                const int NE,
@@ -166,6 +216,7 @@ static void PADiffusionSetup3D(const int Q1D,
 }
 
 static void PADiffusionSetup(const int dim,
+                             const int sdim,
                              const int D1D,
                              const int Q1D,
                              const int NE,
@@ -184,7 +235,8 @@ static void PADiffusionSetup(const int dim,
          return;
       }
 #endif // MFEM_USE_OCCA
-      PADiffusionSetup2D(Q1D, NE, W, J, C, D);
+      if (sdim == 2) { PADiffusionSetup2D<2>(Q1D, NE, W, J, C, D); }
+      if (sdim == 3) { PADiffusionSetup2D<3>(Q1D, NE, W, J, C, D); }
    }
    if (dim == 3)
    {
@@ -224,6 +276,7 @@ void DiffusionIntegrator::SetupPA(const FiniteElementSpace &fes,
    dim = mesh->Dimension();
    ne = fes.GetNE();
    geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
+   const int sdim = mesh->SpaceDimension();
    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
@@ -252,8 +305,8 @@ void DiffusionIntegrator::SetupPA(const FiniteElementSpace &fes,
          }
       }
    }
-   PADiffusionSetup(dim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J, coeff,
-                    pa_data);
+   PADiffusionSetup(dim, sdim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J,
+                    coeff, pa_data);
 }
 
 void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)