Merge branch 'tmop-surf-param' into mixed-bdr-face-integ

eliminated fixed dofs from the system
moved things inside the library
simplified the grad and hessian (wip) computations
fixed the line search

fem/bilinearform.cpp

@@ -1362,6 +1362,9 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
    boundary_integs = mbf->boundary_integs;
    boundary_integs_marker = mbf->boundary_integs_marker;
 
+   boundary_face_integs        = mbf->boundary_face_integs;
+   boundary_face_integs_marker = mbf->boundary_face_integs_marker;
+
    trace_face_integs = mbf->trace_face_integs;
 
    boundary_trace_face_integs = mbf->boundary_trace_face_integs;
@@ -1518,7 +1521,8 @@ void MixedBilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
 void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
 {
    boundary_face_integs.Append(bfi);
-   boundary_face_integs_marker.Append(NULL); // NULL marker means apply everywhere
+   // Active on all boundary faces.
+   boundary_face_integs_marker.Append(NULL);
 }
 
 void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
@@ -1710,7 +1714,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
       bdr_attr_marker = 0;
       for (int k = 0; k < boundary_face_integs.Size(); k++)
       {
-         if (boundary_face_integs_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == nullptr)
          {
             bdr_attr_marker = 1;
             break;
@@ -1748,9 +1752,12 @@ void MixedBilinearForm::Assemble(int skip_zeros)
                if (boundary_face_integs_marker[k] &&
                    (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                           *test_fe2,
+               // boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
+               //                                             *test_fe2,
+               //                                             *ftr, elemmat);
+               boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1,
                                                            *ftr, elemmat);
+
                TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
                mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
             }

fem/bilinearform.hpp

@@ -901,30 +901,31 @@ public:
    }
 
    /// Adds a domain integrator. Assumes ownership of @a bfi.
+   /// Assembled by BilinearFormIntegrator::AssembleElementMatrix2().
    void AddDomainIntegrator(BilinearFormIntegrator *bfi);
-
-   /// Adds a domain integrator. Assumes ownership of @a bfi.
+   /// Adds a domain integrator restricted to certain elements.
    void AddDomainIntegrator(BilinearFormIntegrator *bfi,
                             Array<int> &elem_marker);
 
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
+   /// Assembled by BilinearFormIntegrator::AssembleElementMatrix2().
    void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);
-
-   /// Adds a boundary integrator. Assumes ownership of @a bfi.
-   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
+   /// Adds a boundary integrator restricted to certain boundary elements.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi,
                               Array<int> &bdr_marker);
 
    /// Adds an interior face integrator. Assumes ownership of @a bfi.
    void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds a boundary face integrator. Assumes ownership of @a bfi.
+   /// Assembled by BilinearFormIntegrator::AssembleFaceMatrix().
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi);
-
-   /// Adds a boundary face integrator. Assumes ownership of @a bfi.
+   /// Adds a boundary face integrator restricted to certain boundary faces.
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
                              Array<int> &bdr_marker);
 
    /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
+       Assembled by BilinearFormIntegrator::AssembleFaceMatrix().
 
        This type of integrator assembles terms over all faces of the mesh using
        the face FE from the trial space and the two adjacent volume FEs from
@@ -932,9 +933,9 @@ public:
    void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   /// Assembled by BilinearFormIntegrator::AssembleFaceMatrix().
    void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi);
-
-   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   /// Adds a boundary trace face integrator over certain boundary faces.
    void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi,
                                   Array<int> &bdr_marker);
 
@@ -962,6 +963,7 @@ public:
    /** @brief Access all boundary markers added with AddBdrFaceIntegrator().
        If no marker was specified when the integrator was added, the
        corresponding pointer (to Array<int>) will be NULL. */
+
    Array<Array<int>*> *GetBFBFI_Marker() { return &boundary_face_integs_marker; }
 
    /// Access all integrators added with AddTraceFaceIntegrator().

fem/bilininteg.cpp

@@ -837,6 +837,43 @@ void MixedScalarVectorIntegrator::AssembleElementMatrix2(
    }
 }
 
+void BoundaryMixedForceIntegrator::
+AssembleFaceMatrix(const FiniteElement &trial_fe, const FiniteElement &test_fe,
+                   FaceElementTransformations &Tr, DenseMatrix &elmat)
+{
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      int order = trial_fe.GetOrder() + test_fe.GetOrder();
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+
+   const int nqp_face  = IntRule->GetNPoints();
+   const int vdim      = Q.GetVDim();
+   const int dof_trial = trial_fe.GetDof();
+   const int dof_test  = test_fe.GetDof();
+
+   elmat.SetSize(dof_test * vdim, dof_trial);
+   elmat = 0.0;
+   DenseMatrix loc_force(dof_test, vdim);
+   Vector shape_trial(dof_trial), shape_test(dof_test),
+          Vloc_force(loc_force.Data(), dof_test * vdim);
+   Vector qcoeff(vdim);
+
+   for (int q = 0; q < nqp_face; q++)
+   {
+      const IntegrationPoint &ip_f = IntRule->IntPoint(q);
+      Tr.SetAllIntPoints(&ip_f);
+      const IntegrationPoint &ip_e = Tr.GetElement1IntPoint();
+
+      test_fe.CalcShape(ip_e, shape_test);
+      trial_fe.CalcShape(ip_e, shape_trial);
+      Q.Eval(qcoeff, Tr, ip_f);
+
+      MultVWt(shape_test, qcoeff, loc_force);
+      AddMultVWt(Vloc_force, shape_trial, elmat);
+   }
+}
 
 void GradientIntegrator::AssembleElementMatrix2(
    const FiniteElement &trial_fe, const FiniteElement &test_fe,
@@ -1506,6 +1543,46 @@ void BoundaryMassIntegrator::AssembleFaceMatrix(
    }
 }
 
+void BoundaryVectorMassIntegrator::
+AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2,
+                   FaceElementTransformations &Tr, DenseMatrix &elmat)
+{
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      int order = 2 * el1.GetOrder();
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+
+   const int nqp_face = IntRule->GetNPoints();
+   const int dof = el1.GetDof();
+   elmat.SetSize(dof * vdim);
+   elmat = 0.0;
+
+   mcoeff.SetSize(vdim);
+   shape.SetSize(dof);
+   partelmat.SetSize(dof);
+   for (int q = 0; q < nqp_face; q++)
+   {
+      const IntegrationPoint &ip_f = IntRule->IntPoint(q);
+      // Set the integration point in the face and the neighboring elements
+      Tr.SetAllIntPoints(&ip_f);
+
+      MQ->Eval(mcoeff, Tr, ip_f);
+      el1.CalcShape(Tr.GetElement1IntPoint(), shape);
+      MultVVt(shape, partelmat);
+
+      for (int i = 0; i < vdim; i++)
+      {
+         for (int j = 0; j < vdim; j++)
+         {
+            elmat.AddMatrix(mcoeff(i,j), partelmat, dof*i, dof*j);
+         }
+      }
+   }
+}
+
+
 void ConvectionIntegrator::AssembleElementMatrix(
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
 {

fem/bilininteg.hpp

@@ -56,7 +56,7 @@ public:
 
    /// Method defining partial assembly on NURBS patches.
    /** The result of the partial assembly is stored internally so that it can be
-       used later in the method AddMultNURBSPA(). */
+       used later in the method AdbdMultNURBSPA(). */
    virtual void AssembleNURBSPA(const FiniteElementSpace &fes);
 
    virtual void AssemblePABoundary(const FiniteElementSpace &fes);
@@ -161,7 +161,12 @@ public:
    /** Compute the local matrix representation of a bilinear form
        $a(u,v)$ defined on different trial (given by $u$) and test
        (given by $v$) spaces. The rows in the local matrix correspond
-       to the test dofs and the columns -- to the trial dofs. */
+       to the test dofs and the columns -- to the trial dofs.
+
+       - Used with MixedBilinearForm::AddDomainIntegrator() to assemble over
+         volumetric elements.
+       - Used with MixedBilinearForm::AddBoundaryIntegrator() to assemble over
+         boundary elements. */
    virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                        const FiniteElement &test_fe,
                                        ElementTransformation &Trans,
@@ -174,11 +179,25 @@ public:
                                     const FiniteElementSpace &fes,
                                     SparseMatrix*& smat);
 
+   /** @brief  Integration on faces, interior or boundary, when one uses
+               information from the neighbor volumetric elements. Depending on
+               the context, the arguments are used in different ways. The output
+               matrix @a elmat is always based on the volumetric DOFs.
+
+       - Used with BilinearForm::AddInteriorFaceIntegrator(), where @a el1 and
+         @a el2 are the FiniteElements for both sides of the internal face.
+       - Used with BilinearForm::AddBdrFaceIntegrator(), where @a el1 is for the
+         volumetric neighbor of the boundary face, while @a el2 is dummy and
+         should not be used.
+       - Used with MixedBilinearForm::AddBdrFaceIntegrator(), where @a el1 and
+         @a el2 are the trial and test FiniteElements, respectively, for the
+         volumetric neighbor element of the boundary face. */
    virtual void AssembleFaceMatrix(const FiniteElement &el1,
                                    const FiniteElement &el2,
                                    FaceElementTransformations &Trans,
                                    DenseMatrix &elmat);
 
+
    virtual void AssembleFaceMatrix(const FiniteElement &trial_fe1,
                                    const FiniteElement &test_fe1,
                                    const FiniteElement &trial_fe2,
@@ -2113,6 +2132,30 @@ protected:
    }
 };
 
+/** Integrates (Q*u, v) over boundary faces, where u=(u1,...,un) is a vector
+    trial function (each component is a scalar FE), Q is a VectorCoefficient,
+    Q*u is the Hadamard product, and v is a scalar test function.
+
+    - SetIntRule() is expected to take a face-based quadrature rule.
+    - The VectorCoefficient Q is evaluated through a FaceElementTransformation,
+      allowing to utilize both face-based information (e.g. face normals) and
+      volumetric element information (e.g. volumetric deformation). */
+class BoundaryMixedForceIntegrator : public BilinearFormIntegrator
+{
+protected:
+   VectorCoefficient &Q;
+
+public:
+   BoundaryMixedForceIntegrator(VectorCoefficient &vc) : Q(vc) { }
+
+   /// Expected use is with MixedBilinearForm::AddBdrFaceIntegrator(), where
+   /// @a el1 and @a el2 are for the (mixed) volumetric neighbor of the face.
+  void AssembleFaceMatrix(const FiniteElement &trial_fe,
+                          const FiniteElement &test_fe,
+                          FaceElementTransformations &Tr,
+                          DenseMatrix &elmat) override;
+};
+
 /** Class for integrating the bilinear form $a(u,v) := (Q \nabla u, v)$ where $Q$ is a
     scalar coefficient, $u$ is in ($H^1$), and $v$ is a vector with components
     $v_i$ in ($H^1$) or ($L^2$).
@@ -2464,9 +2507,17 @@ protected:
    }
 };
 
-/** Mass integrator $(u, v)$ restricted to the boundary of a domain */
+/** @brief Mass integrator $(u, v)$ restricted to the boundary of a domain.
+
+    - SetIntRule() is expected to take a face-based quadrature rule.
+    - The Coefficient Q is evaluated through a FaceElementTransformation,
+      allowing to utilize both face-based information (e.g. face normals) and
+      volumetric element information (e.g. volumetric deformation). */
 class BoundaryMassIntegrator : public MassIntegrator
 {
+protected:
+   int nf;
+
 public:
    BoundaryMassIntegrator(Coefficient &q) : MassIntegrator(q) { }
 
@@ -2596,15 +2647,16 @@ public:
     by scalar FE through standard transformation. */
 class VectorMassIntegrator: public BilinearFormIntegrator
 {
+protected:
    int vdim = -1, Q_order = 0;
    Vector shape, te_shape, vec;
    DenseMatrix partelmat;
    DenseMatrix mcoeff;
 
-protected:
    Coefficient *Q = nullptr;
    VectorCoefficient *VQ = nullptr;
    MatrixCoefficient *MQ = nullptr;
+
    // PA extension
    const DofToQuad *maps;         ///< Not owned
    const GeometricFactors *geom;  ///< Not owned
@@ -2661,6 +2713,36 @@ public:
                          (int, int, int));
 };
 
+/** @brief Integrates (Q u, v) over the boundary of the domain, where
+    u=(u1,...,un) and v=(v1,...,vn); ui and vi are defined by scalar FE through
+    standard transformation.
+
+    - SetIntRule() is expected to take a face-based quadrature rule.
+    - The Coefficient Q is evaluated through a FaceElementTransformation,
+      allowing to utilize both face-based information (e.g. face normals) and
+      volumetric element information (e.g. volumetric deformation). */
+class BoundaryVectorMassIntegrator : public VectorMassIntegrator
+{
+protected:
+   int nf;
+
+public:
+   /// The given MatrixCoefficient fully couples the vector components, i.e.,
+   /// the local (dof x vdim) matrices have no zero blocks.
+   BoundaryVectorMassIntegrator(MatrixCoefficient &mc)
+      : VectorMassIntegrator(mc) { }
+
+   /// Expected use is with BilinearForm::AddBdrFaceIntegrator(), where @a el1
+   /// is for the volumetric neighbor of the boundary face, @a el2 is not used.
+   void AssembleFaceMatrix(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Tr,
+                           DenseMatrix &elmat) override;
+
+   void AssemblePABoundaryFaces(const FiniteElementSpace &fes) override;
+   void AddMultPA(const Vector &x, Vector &y) const override;
+};
+
 
 /** Class for integrating $(\nabla \cdot u, p)$ where $u$ is a vector field given by
     VectorFiniteElement through Piola transformation (for Raviart-Thomas elements); $p$ is

fem/integ/bilininteg_bdrvecmass_pa.cpp

@@ -0,0 +1,196 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../../general/forall.hpp"
+#include "../bilininteg.hpp"
+
+namespace mfem
+{
+
+template<const int T_D1D = 0, const int T_Q1D = 0>
+static void PABdrVectorMassApply2D(const int NE,
+                                   const Array<double> &B_,
+                                   const Array<double> &Bt_,
+                                   const Vector &d_,
+                                   const Vector &x_,
+                                   Vector &y_,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int VDIM = 2;
+   MFEM_VERIFY(T_D1D ? T_D1D : d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(T_Q1D ? T_Q1D : q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+   auto B  = Reshape(B_.Read(), Q1D, D1D);
+   auto Bt = Reshape(Bt_.Read(), D1D, Q1D);
+   auto D  = Reshape(d_.Read(), VDIM, VDIM, Q1D, Q1D, NE);
+   auto x  = Reshape(x_.Read(), D1D, D1D, VDIM, NE);
+   auto y  = Reshape(y_.ReadWrite(), D1D, D1D, VDIM, NE);
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      // the following variables are evaluated at compile time
+      constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+      constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
+
+      double sol_xy[max_Q1D][max_Q1D][VDIM];
+      for (int c = 0; c < VDIM; ++c)
+      {
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int cc = 0; cc < VDIM; ++cc)
+               {
+                  sol_xy[qy][qx][cc] = 0.0;
+               }
+            }
+         }
+
+         // dof -> quad.
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            double sol_x[max_Q1D][VDIM];
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               for (int cc = 0; cc < VDIM; ++cc)
+               {
+                  sol_x[qy][cc] = 0.0;
+               }
+            }
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               for (int cc = 0; cc < VDIM; cc++)
+               {
+                  const double s = x(dx,dy,cc,e);
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     sol_x[qx][cc] += B(qx,dx)* s;
+                  }
+               }
+            }
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               const double d2q = B(qy,dy);
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  for (int cc = 0; cc < VDIM; cc++)
+                  {
+                     sol_xy[qy][qx][cc] += d2q * sol_x[qx][cc];
+                  }
+               }
+            }
+         }
+
+         // quad data.
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int cc = 0; cc < VDIM; cc++)
+               {
+                  sol_xy[qy][qx][cc] *= D(c,cc,qx,qy,e);
+               }
+            }
+         }
+
+         // quad -> dof.
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            double sol_x[max_D1D][VDIM];
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               for (int cc = 0; cc < VDIM; cc++)
+               {
+                  sol_x[dx][cc] = 0.0;
+               }
+            }
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int cc = 0; cc < VDIM; cc++)
+               {
+                  const double s = sol_xy[qy][qx][cc];
+                  for (int dx = 0; dx < D1D; ++dx)
+                  {
+                     sol_x[dx][cc] += Bt(dx,qx) * s;
+                  }
+               }
+            }
+            for (int dy = 0; dy < D1D; ++dy)
+            {
+               const double q2d = Bt(dy,qy);
+               for (int dx = 0; dx < D1D; ++dx)
+               {
+                  for (int cc = 0; cc < VDIM; cc++)
+                  {
+                     y(dx,dy,c,e) += q2d * sol_x[dx][cc];
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
+template<const int T_D1D = 0, const int T_Q1D = 0>
+static void PABdrVectorMassApply3D(const int NE,
+                                   const Array<double> &B_,
+                                   const Array<double> &Bt_,
+                                   const Vector &q_,
+                                   const Vector &x_,
+                                   Vector &y_,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
+{
+
+}
+
+void BoundaryVectorMassIntegrator::
+     AssemblePABoundaryFaces(const FiniteElementSpace &fes)
+{
+   nf = fes.GetNFbyType(FaceType::Boundary);
+   if (nf == 0) { return; }
+
+   Mesh *mesh = fes.GetMesh();
+   const FiniteElement &el =
+      *fes.GetTraceElement(0, fes.GetMesh()->GetFaceGeometry(0));
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      int order = 2 * fes.FEColl()->GetOrder();
+      ir = &IntRules.Get(mesh->GetFaceGeometry(0), order);
+   }
+
+   dim = mesh->Dimension();
+
+   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
+}
+
+void BoundaryVectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
+{
+   if (dim == 2)
+   {
+      MFEM_VERIFY(vdim == 2, "Not implemented for genereal vdim");
+      return PABdrVectorMassApply2D(ne, maps->B, maps->Bt, pa_data,
+                                    x, y, dofs1D, quad1D);
+   }
+   if (dim == 3)
+   {
+      MFEM_VERIFY(vdim == 3, "Not implemented for genereal vdim");
+      return PABdrVectorMassApply3D(ne, maps->B, maps->Bt, pa_data,
+                                    x, y, dofs1D, quad1D);
+   }
+}
+
+} // namespace mfem

fem/tmop.cpp

@@ -19,7 +19,7 @@
 namespace mfem
 {
 
-/* AD related definitions below ========================================*/
+/* AD related definitions below */
 
 /// MFEM native AD-type for first derivatives
 using AD1Type = future::dual<real_t, real_t>;
@@ -485,7 +485,216 @@ void ADHessian(std::function<AD2Type(std::vector<AD2Type>&,
    }
    return;
 }
-/* end AD related definitions ========================================*/
+/* end AD related definitions */
+
+// I1 = |M|^2 / det(M).
+double TMOP_QualityMetric::Dim2Invariant1(const DenseMatrix &M)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   return M.FNorm2() / M.Det();
+}
+
+// I2 = det(M).
+double TMOP_QualityMetric::Dim2Invariant2(const DenseMatrix &M)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   return M.Det();
+}
+
+// dI1_dM = [ 2 det(M) M - |M|^2 adj(M)^T ] / det(M)^2.
+void TMOP_QualityMetric::Dim2Invariant1_dM(const DenseMatrix &M,
+                                           DenseMatrix &dM)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   const double fnorm2 = M.FNorm2(), det = M.Det();
+
+   Dim2Invariant2_dM(M, dM);
+   dM *= - fnorm2/(det*det);
+   dM.Add(2.0/det, M);
+}
+
+// dI2_dM = d(det(M))_dM = adj(M)^T.
+void TMOP_QualityMetric::Dim2Invariant2_dM(const DenseMatrix &M,
+                                           DenseMatrix &dM)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   dM(0, 0) =  M(1, 1); dM(0, 1) = -M(1, 0);
+   dM(1, 0) = -M(0, 1); dM(1, 1) =  M(0, 0);
+}
+
+// (dI1_dM)_d(Mij) = d[(2 det(M) M - |M|^2 adj(M)^T) / det(M)^2]_d[Mij].
+void TMOP_QualityMetric::Dim2Invariant1_dMdM(const DenseMatrix &M, int i, int j,
+                                             DenseMatrix &dMdM)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   // Compute d(det(M))_d(Mij), d(|M|^2)_d(Mij).
+   DenseMatrix dI(2);
+   Dim2Invariant2_dM(M, dI);
+   const double ddet   = dI(i,j);
+   const double dfnorm2 = 2.0 * M(i,j);
+
+   const double det    = M.Det();
+   const double det2   = det * det;
+   const double fnorm2 = M.FNorm2();
+
+   DenseMatrix dM(2); dM = 0.0; dM(i, j) = 1.0;
+   DenseMatrix ddI(2);
+   Dim2Invariant2_dMdM(M, i, j, ddI);
+   for (int r = 0; r < 2; r++)
+   {
+      for (int c = 0; c < 2; c++)
+      {
+         dMdM(r,c) =
+             (det2 *
+                  (2.0 * ddet * M(r,c) + 2.0 * det * dM(r,c)
+                   - dfnorm2 * dI(r,c) - fnorm2 * ddI(r,c))
+              - 2.0 * det * ddet *
+                    (2.0 * det * M(r,c) - fnorm2 * dI(r,c)) ) / (det2 * det2);
+      }
+   }
+}
+
+// (dI2_dM)_d(Mij) = ...
+void TMOP_QualityMetric::Dim2Invariant2_dMdM(const DenseMatrix &M, int i, int j,
+                                             DenseMatrix &dMdM)
+{
+   MFEM_ASSERT(M.Height() == 2 && M.Width() == 2, "Incorrect dimensions!");
+
+   dMdM = 0.0;
+   dMdM(1-i,1-j) = (i == j) ? 1.0 : -1.0;
+}
+
+// I1 = |M|^2/ det(M)^(2/3).
+double TMOP_QualityMetric::Dim3Invariant1(const DenseMatrix &M)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   const double fnorm = M.FNorm(), det = M.Det();
+   return fnorm * fnorm / pow(det, 2.0/3.0);
+}
+
+
+// I2 = |adj(M)|^2 / det(M)^(4/3).
+double TMOP_QualityMetric::Dim3Invariant2(const DenseMatrix &M)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   DenseMatrix Madj(3);
+   CalcAdjugate(M, Madj);
+
+   const double fnorm = Madj.FNorm(), det = M.Det();
+   return fnorm * fnorm / pow(det, 4.0/3.0);
+}
+
+
+// I3 = det(M).
+double TMOP_QualityMetric::Dim3Invariant3(const DenseMatrix &M)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   return M.Det();
+}
+
+// dI1_dM = [ 2 det(M) M - 2/3 |M|^2 det(M)^(-1/3) adj(M)^T ] / det(M)^4/3.
+void TMOP_QualityMetric::Dim3Invariant1_dM(const DenseMatrix &M,
+                                           DenseMatrix &dM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   DenseMatrix Madj(3);
+   CalcAdjugate(M, Madj);
+   const double fnorm = M.FNorm(), det = M.Det();
+
+   Dim3Invariant3_dM(M, dM);
+   dM *= -(2./3.)* fnorm * fnorm * pow(det, -1./3.);
+   dM.Add(2.0 * pow(det, 2./3.), M);
+   dM *= 1.0 / pow(det, 4./3.);
+}
+
+// dI2_dM = [ -4/3 |adj(M)|^2  det(M)^(1/3) adj(M)^T ] / det(M)^(8/3).
+void TMOP_QualityMetric::Dim3Invariant2_dM(const DenseMatrix &M,
+                                           DenseMatrix &dM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   DenseMatrix Madj(3);
+   // dM will have Madj^t because it is the third invariant's derivative.
+   CalcAdjugate(M, Madj);
+   const double fnorm = Madj.FNorm(), det = M.Det();
+
+   Dim3Invariant3_dM(M, dM);
+   dM *= -(4./3.)* fnorm * fnorm * pow(det, 1./3.);
+   dM *= 1.0 / (pow(det, 8./3.));
+}
+
+// dI3_dM = d(det(M))_dM = adj(M)^T.
+void TMOP_QualityMetric::Dim3Invariant3_dM(const DenseMatrix &M,
+                                           DenseMatrix &dM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   CalcAdjugateTranspose(M, dM);
+}
+
+void TMOP_QualityMetric::Dim3Invariant1_dMdM(const DenseMatrix &M, int i, int j,
+                                             DenseMatrix &dMdM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   DenseMatrix dI(3);
+   Dim3Invariant3_dM(M, dI);
+   const double fnorm  = M.FNorm(), det = M.Det();
+
+   DenseMatrix dM(3); dM = 0.0; dM(i, j) = 1.0;
+   for (int r = 0; r < 3; r++)
+   {
+      for (int c = 0; c < 3; c++)
+      {
+         dMdM(r,c) = (2.0 * det * det * dM(r,c)
+                       + dI(i,j) * (10./9.) * fnorm * fnorm * dI(r,c)
+                       - (4./3.) * dI(i,j) * det * M(r,c)
+                       - (4./3.) * det * M(i,j) * dI(r,c))
+                      / pow(det, 8./3.);
+      }
+   }
+}
+
+void TMOP_QualityMetric::Dim3Invariant2_dMdM(const DenseMatrix &M, int i, int j,
+                                             DenseMatrix &dMdM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   DenseMatrix dI(3);
+   Dim3Invariant3_dM(M, dI);
+   DenseMatrix Madj(3);
+   CalcAdjugate(M, Madj);
+   const double det   = M.Det();
+   const double fnorm = Madj.FNorm();
+
+   DenseMatrix dM(3); dM = 0.0; dM(i, j) = 1.0;
+   for (int r = 0; r < 3; r++)
+   {
+      for (int c = 0; c < 3; c++)
+      {
+         dMdM(r,c) = (28./9.) * fnorm * fnorm *
+                      det * det * dI(i,j) * dI(r,c) / pow(det, 16./3.);
+      }
+   }
+}
+
+// (dI3_dM)_d(Mij) = 0.
+void TMOP_QualityMetric::Dim3Invariant3_dMdM(const DenseMatrix &M, int i, int j,
+                                             DenseMatrix &dMdM)
+{
+   MFEM_ASSERT(M.Height() == 3 && M.Width() == 3, "Incorrect dimensions!");
+
+   dMdM(i, j) = 0.0;
+}
 
 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
 
@@ -975,12 +1184,24 @@ real_t TMOP_Metric_002::EvalWMatrixForm(const DenseMatrix &Jpt) const
 
 real_t TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
 {
+   if (use_old_invariants_code)
+   {
+      return 0.5 * Dim2Invariant1(Jpt) - 1.0;
+   }
+
    ie.SetJacobian(Jpt.GetData());
    return 0.5 * ie.Get_I1b() - 1.0;
 }
 
 void TMOP_Metric_002::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
 {
+   if (use_old_invariants_code)
+   {
+      Dim2Invariant1_dM(Jpt, P);
+      P *= 0.5;
+      return;
+   }
+
    ie.SetJacobian(Jpt.GetData());
    P.Set(0.5, ie.Get_dI1b());
 }
@@ -990,11 +1211,83 @@ void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
                                 const real_t weight,
                                 DenseMatrix &A) const
 {
+   if (use_old_invariants_code)
+   {
+      const int dof = DS.Height(), dim = DS.Width();
+      DenseMatrix dI1_dMdM(dim);
+
+      // The first two go over the rows and cols of dP_dJ where P = dW_dJ.
+      for (int r = 0; r < dim; r++)
+      {
+         for (int c = 0; c < dim; c++)
+         {
+            Dim2Invariant1_dMdM(Jpt, r, c, dI1_dMdM);
+
+            // Compute each entry of d(Prc)_dJ.
+            for (int rr = 0; rr < dim; rr++)
+            {
+               for (int cc = 0; cc < dim; cc++)
+               {
+                  const double entry_rr_cc = 0.5 * dI1_dMdM(rr,cc);
+
+                  for (int i = 0; i < dof; i++)
+                  {
+                     for (int j = 0; j < dof; j++)
+                     {
+                        A(i+r*dof, j+rr*dof) +=
+                            weight * DS(i, c) * DS(j, cc) * entry_rr_cc;
+                     }
+                  }
+               }
+            }
+         }
+      }
+      return;
+   }
+
    ie.SetJacobian(Jpt.GetData());
    ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
    ie.Assemble_ddI1b(0.5*weight, A.GetData());
 }
 
+void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
+                                const DenseMatrix &DS,
+                                const DenseMatrix &dx_dt,
+                                const real_t weight,
+                                DenseMatrix &A) const
+{
+   const int dof = DS.Height(), dim = DS.Width();
+   DenseMatrix dI1_dMdM(dim);
+
+   // The first two go over the rows and cols of dP_dJ where P = dW_dJ.
+   for (int r = 0; r < dim; r++)
+   {
+      for (int c = 0; c < dim; c++)
+      {
+         Dim2Invariant1_dMdM(Jpt, r, c, dI1_dMdM);
+
+         // Compute each entry of d(Prc)_dJ.
+         for (int rr = 0; rr < dim; rr++)
+         {
+            for (int cc = 0; cc < dim; cc++)
+            {
+               const double entry_rr_cc = 0.5 * dI1_dMdM(rr,cc);
+
+               for (int i = 0; i < dof; i++)
+               {
+                  for (int j = 0; j < dof; j++)
+                  {
+                     A(i+r*dof, j+rr*dof) +=
+                         weight * DS(i, c)  * dx_dt(i, r)
+                                * DS(j, cc) * dx_dt(j, rr) * entry_rr_cc;
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 real_t TMOP_Metric_004::EvalW(const DenseMatrix &Jpt) const
 {
    ie.SetJacobian(Jpt.GetData());
@@ -4086,6 +4379,8 @@ void TMOP_Integrator::GetSurfaceFittingErrors(const Vector &d_loc,
 {
    MFEM_VERIFY(periodic == false,
                "Fitting is not supported for periodic meshes.");
+   MFEM_VERIFY(tan_analytic_surf == nullptr,
+               "Fitting is not supported for parametrized optimization.");
 
    Vector pos(d_loc.Size());
    if (x_0) { add(*x_0, d_loc, pos); }
@@ -4219,9 +4514,21 @@ real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
    const IntegrationRule &ir = EnergyIntegrationRule(el);
 
+   // Convert parametric coordinates to physical to compute Jpt.
+   // Needed as elfun has the surface parameters for some entries.
+   Vector convertedX(elfun);
+   if (tan_analytic_surf)
+   {
+      Array<int> vdofs;
+      fes_mesh->GetElementVDofs(T.ElementNo, vdofs);
+      tan_analytic_surf->ConvertParamToPhys(vdofs, elfun, convertedX);
+      // Use converted coordinates for PMatI
+      PMatI.UseExternalData(convertedX.GetData(), dof, dim);
+   }
+
    energy = 0.0;
    DenseTensor Jtr(dim, dim, ir.GetNPoints());
-   targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
+   targetC->ComputeElementTargets(el_id, el, ir, convertedX, Jtr);
 
    // Limited case.
    Vector shape, p, p0, d_vals;
@@ -4501,6 +4808,171 @@ real_t TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
    return energy;
 }
 
+// Helper functions used for computing derivatives in TMOP_Integrator.
+namespace
+{
+
+// Stores derivatives of param -> phys map.
+struct TangentialNodeDerivatives
+{
+   DenseMatrix jacobian;
+   DenseTensor hessian;
+
+   TangentialNodeDerivatives() = default;
+
+   explicit TangentialNodeDerivatives(int dim)
+      : jacobian(dim, dim), hessian(dim, dim, dim)
+   {
+      jacobian = 0.0;
+      hessian = 0.0;
+   }
+};
+
+// Extracts the parameter-space coordinates stored for one element node.
+void GetNodeParam(const Vector &elfun, int dof, int dim, int node,
+                  double *param)
+{
+   for (int d = 0; d < dim; d++)
+   {
+      param[d] = elfun(node + d * dof);
+   }
+}
+
+// Builds per-node surface Jacobian and Hessian data for tangential mapping.
+void BuildTangentialNodeDerivatives(const AnalyticCompositeSurface &surfaces,
+                                    const Array<int> &vdofs,
+                                    const Vector &elfun,
+                                    int dof,
+                                    int dim,
+                                    std::vector<TangentialNodeDerivatives> &data)
+{
+   data.assign(dof, TangentialNodeDerivatives(dim));
+
+   for (int i = 0; i < dof; i++)
+   {
+      TangentialNodeDerivatives &node = data[i];
+      const AnalyticSurface *surf = surfaces.GetSurface(vdofs[i]);
+      if (!surf)
+      {
+         for (int d = 0; d < dim; d++)
+         {
+            node.jacobian(d, d) = 1.0;
+         }
+         continue;
+      }
+
+      double param[3] = { 0.0, 0.0, 0.0 };
+      GetNodeParam(elfun, dof, dim, i, param);
+
+      DenseMatrix deriv_1;
+      DenseTensor deriv_2;
+      surf->Deriv_1(param, deriv_1);
+      surf->Deriv_2(param, deriv_2);
+
+      for (int r = 0; r < dim; r++)
+      {
+         for (int c = 0; c < surf->NumParams(); c++)
+         {
+            node.jacobian(r, c) = deriv_1(r, c);
+         }
+      }
+
+      for (int r = 0; r < dim; r++)
+      {
+         for (int c1 = 0; c1 < surf->NumParams(); c1++)
+         {
+            for (int c2 = 0; c2 < surf->NumParams(); c2++)
+            {
+               node.hessian(r, c1, c2) = deriv_2(r, c1, c2);
+            }
+         }
+      }
+   }
+}
+
+// Maps an element vector from physical coordinates to surface parameters.
+void TransformTangentialElementVector(
+   const std::vector<TangentialNodeDerivatives> &data,
+   int dof,
+   int dim,
+   const Vector &phys,
+   Vector &param)
+{
+   param.SetSize(dof * dim);
+   param = 0.0;
+
+   for (int i = 0; i < dof; i++)
+   {
+      for (int a = 0; a < dim; a++)
+      {
+         real_t val = 0.0;
+         for (int c = 0; c < dim; c++)
+         {
+            val += data[i].jacobian(c, a) * phys(i + c * dof);
+         }
+         param(i + a * dof) = val;
+      }
+   }
+}
+
+// Applies the tangential chain rule to convert element Hessian data.
+void TransformTangentialElementGrad(
+   const std::vector<TangentialNodeDerivatives> &data,
+   int dof,
+   int dim,
+   const DenseMatrix &phys_hess,
+   const Vector &phys_grad,
+   DenseMatrix &param_hess)
+{
+   param_hess.SetSize(dof * dim);
+   param_hess = 0.0;
+
+   for (int i = 0; i < dof; i++)
+   {
+      for (int a = 0; a < dim; a++)
+      {
+         for (int j = 0; j < dof; j++)
+         {
+            for (int b = 0; b < dim; b++)
+            {
+               real_t val = 0.0;
+               for (int c = 0; c < dim; c++)
+               {
+                  const real_t jac_ic = data[i].jacobian(c, a);
+                  if (jac_ic == 0.0) { continue; }
+                  for (int r = 0; r < dim; r++)
+                  {
+                     const real_t jac_jr = data[j].jacobian(r, b);
+                     if (jac_jr == 0.0) { continue; }
+                     val += jac_ic * phys_hess(i + c * dof, j + r * dof) *
+                            jac_jr;
+                  }
+               }
+               param_hess(i + a * dof, j + b * dof) += val;
+            }
+         }
+      }
+   }
+
+   for (int i = 0; i < dof; i++)
+   {
+      for (int a = 0; a < dim; a++)
+      {
+         for (int b = 0; b < dim; b++)
+         {
+            real_t val = 0.0;
+            for (int c = 0; c < dim; c++)
+            {
+               val += phys_grad(i + c * dof) * data[i].hessian(c, a, b);
+            }
+            param_hess(i + a * dof, i + b * dof) += val;
+         }
+      }
+   }
+}
+
+}
+
 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
                                             ElementTransformation &T,
                                             const Vector &d_el, Vector &elvect)
@@ -4536,6 +5008,7 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
                                                  Vector &elvect)
 {
    const int dof = el.GetDof(), dim = el.GetDim();
+
    const int el_id = T.ElementNo;
 
    // Form the Vector of node positions, depending on what's the input.
@@ -4560,17 +5033,31 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
    Jpt.SetSize(dim);
    P.SetSize(dim);
    PMatI.UseExternalData(elfun.GetData(), dof, dim);
-   elvect.SetSize(dof*dim);
-   PMatO.UseExternalData(elvect.GetData(), dof, dim);
+   Vector elvect_phys(dof * dim);
+   PMatO.UseExternalData(elvect_phys.GetData(), dof, dim);
 
    const IntegrationRule &ir = ActionIntegrationRule(el);
    const int nqp = ir.GetNPoints();
 
-   elvect = 0.0;
+   // Convert parametric coordinates to physical to compute Jpt.
+   // Needed as elfun has the surface parameters for some entries.
+   Vector convertedX(elfun);
+   MFEM_VERIFY(tan_analytic_surf, "only for analytic surfaces");
+   Array<int> vdofs;
+   fes_mesh->GetElementVDofs(T.ElementNo, vdofs);
+   tan_analytic_surf->ConvertParamToPhys(vdofs, elfun, convertedX);
+   // Use converted coordinates for PMatI
+   PMatI.UseExternalData(convertedX.GetData(), dof, dim);
+
+   std::vector<TangentialNodeDerivatives> tan_node_data;
+   BuildTangentialNodeDerivatives(*tan_analytic_surf, vdofs, elfun, dof, dim,
+                                  tan_node_data);
+
+   elvect_phys = 0.0;
    Vector weights(nqp);
    DenseTensor Jtr(dim, dim, nqp);
    DenseTensor dJtr(dim, dim, dim*nqp);
-   targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, ir, convertedX, Jtr);
 
    // Limited case.
    DenseMatrix pos0;
@@ -4614,7 +5101,7 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
       Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
       if (exact_action)
       {
-         targetC->ComputeElementTargetsGradient(ir, elfun, *Tpr, dJtr);
+         targetC->ComputeElementTargetsGradient(ir, convertedX, *Tpr, dJtr);
       }
    }
 
@@ -4630,6 +5117,7 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
       real_t weight_m = weights(q) * metric_normal;
 
       el.CalcDShape(ip, DSh);
+      // change comes here for DSh?
       Mult(DSh, Jrt, DS);
       MultAtB(PMatI, DS, Jpt);
 
@@ -4637,7 +5125,21 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
 
       if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
       P *= weight_m;
-      AddMultABt(DS, P, PMatO); // w_q det(W) dmu/dx : dA/dx Winv
+
+      // Assemble the physical-coordinate residual first. The constrained
+      // parameter residual is obtained by a local chain-rule transform below.
+      for (int i = 0; i < dof; i++)
+      {
+         for (int c = 0; c < dim; c++)
+         {
+            real_t d = 0.0;
+            for (int k = 0; k < dim; k++)
+            {
+               d += DS(i, k) * P(c, k);
+            }
+            PMatO(i, c) += d;
+         }
+      }
 
       if (exact_action)
       {
@@ -4703,6 +5205,8 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
    if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, *Tpr, ir, weights, PMatO); }
    if (surf_fit_gf || surf_fit_pos) { AssembleElemVecSurfFit(el, *Tpr, PMatO); }
 
+   TransformTangentialElementVector(tan_node_data, dof, dim, elvect_phys, elvect);
+
    delete Tpr;
 }
 
@@ -4739,10 +5243,28 @@ void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
    const IntegrationRule &ir = GradientIntegrationRule(el);
    const int nqp = ir.GetNPoints();
 
-   elmat = 0.0;
+   // Convert parametric coordinates to physical to compute Jpt.
+   // Needed as elfun has the surface parameters for some entries.
+   Vector convertedX(elfun);
+   MFEM_VERIFY(tan_analytic_surf, "only for analytic surfaces");
+   Array<int> vdofs;
+   fes_mesh->GetElementVDofs(T.ElementNo, vdofs);
+   tan_analytic_surf->ConvertParamToPhys(vdofs, elfun, convertedX);
+   // Use converted coordinates for PMatI
+   PMatI.UseExternalData(convertedX.GetData(), dof, dim);
+
+   std::vector<TangentialNodeDerivatives> tan_node_data;
+   BuildTangentialNodeDerivatives(*tan_analytic_surf, vdofs, elfun, dof, dim,
+                                  tan_node_data);
+
+   DenseMatrix elmat_phys(dof * dim);
+   elmat_phys = 0.0;
+   Vector elvect_phys(dof * dim);
+   DenseMatrix PMatGrad(elvect_phys.GetData(), dof, dim);
+   elvect_phys = 0.0;
    Vector weights(nqp);
    DenseTensor Jtr(dim, dim, nqp);
-   targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
+   targetC->ComputeElementTargets(el_id, el, ir, convertedX, Jtr);
 
    // Limited case.
    DenseMatrix pos0, hess;
@@ -4800,7 +5322,22 @@ void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
 
       if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
 
-      metric->AssembleH(Jpt, DS, weight_m, elmat);
+      metric->AssembleH(Jpt, DS, weight_m, elmat_phys);
+
+      metric->EvalP(Jpt, P);
+      P *= weight_m;
+      for (int i = 0; i < dof; i++)
+      {
+         for (int c = 0; c < dim; c++)
+         {
+            real_t d = 0.0;
+            for (int k = 0; k < dim; k++)
+            {
+               d += DS(i, k) * P(c, k);
+            }
+            PMatGrad(i, c) += d;
+         }
+      }
 
       // TODO: derivatives of adaptivity-based targets.
 
@@ -4810,6 +5347,10 @@ void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
          PMatI.MultTranspose(shape, p);
          pos0.MultTranspose(shape, p0);
          weight_m = weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
+         Vector grad_lim(dim);
+         lim_func->Eval_d1(p, p0, d_vals(q), grad_lim);
+         grad_lim *= weight_m;
+         AddMultVWt(shape, grad_lim, PMatGrad);
          lim_func->Eval_d2(p, p0, d_vals(q), hess);
          for (int i = 0; i < dof; i++)
          {
@@ -4821,7 +5362,7 @@ void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
                {
                   for (int d2 = 0; d2 < dim; d2++)
                   {
-                     elmat(d1*dof + i, d2*dof + j) += w * hess(d1, d2);
+                     elmat_phys(d1 * dof + i, d2 * dof + j) += w * hess(d1, d2);
                   }
                }
             }
@@ -4829,8 +5370,17 @@ void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
       }
    }
 
-   if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, *Tpr, ir, weights, elmat); }
-   if (surf_fit_gf || surf_fit_pos) { AssembleElemGradSurfFit(el, *Tpr, elmat);}
+   if (adapt_lim_gf)
+   {
+      AssembleElemGradAdaptLim(el, *Tpr, ir, weights, elmat_phys);
+   }
+   if (surf_fit_gf || surf_fit_pos)
+   {
+      AssembleElemGradSurfFit(el, *Tpr, elmat_phys);
+   }
+
+   TransformTangentialElementGrad(tan_node_data, dof, dim, elmat_phys,
+                                  elvect_phys, elmat);
 
    delete Tpr;
 }

fem/tmop.hpp

@@ -14,6 +14,7 @@
 
 #include "../linalg/invariants.hpp"
 #include "nonlininteg.hpp"
+#include "tmop/AnalyticalSurface.hpp"
 #include "../linalg/dual.hpp"
 
 namespace mfem
@@ -34,6 +35,47 @@ protected:
        for TMOP_QualityMetric%s, because it is not used. */
    void SetTransformation(ElementTransformation &) { }
 
+   /// First invariant of the given 2x2 matrix @a M.
+   static double Dim2Invariant1(const DenseMatrix &M);
+   /// Second invariant of the given 2x2 matrix @a M.
+   static double Dim2Invariant2(const DenseMatrix &M);
+
+   /// 1st derivative of the first invariant for the given 2x2 matrix @a M.
+   static void Dim2Invariant1_dM(const DenseMatrix &M, DenseMatrix &dM);
+   /// 1st derivative of the second invariant for the given 2x2 matrix @a M.
+   static void Dim2Invariant2_dM(const DenseMatrix &M, DenseMatrix &dM);
+
+   /// 2nd derivative of the first invariant for the given 2x2 matrix @a M.
+   static void Dim2Invariant1_dMdM(const DenseMatrix &M, int i, int j,
+                                   DenseMatrix &dMdM);
+   /// 2nd derivative of the second invariant for the given 2x2 matrix @a M.
+   static void Dim2Invariant2_dMdM(const DenseMatrix &M, int i, int j,
+                                   DenseMatrix &dMdM);
+
+   /// First invariant of the given 3x3 matrix @a M.
+   static double Dim3Invariant1(const DenseMatrix &M);
+   /// Second invariant of the given 3x3 matrix @a M.
+   static double Dim3Invariant2(const DenseMatrix &M);
+   /// Third invariant of the given 3x3 matrix @a M.
+   static double Dim3Invariant3(const DenseMatrix &M);
+
+   /// 1st derivative of the first invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant1_dM(const DenseMatrix &M, DenseMatrix &dM);
+   /// 1st derivative of the second invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant2_dM(const DenseMatrix &M, DenseMatrix &dM);
+   /// 1st derivative of the third invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant3_dM(const DenseMatrix &M, DenseMatrix &dM);
+
+   /// 2nd derivative of the first invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant1_dMdM(const DenseMatrix &M, int i, int j,
+                                   DenseMatrix &dMdM);
+   /// 2nd derivative of the second invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant2_dMdM(const DenseMatrix &M, int i, int j,
+                                   DenseMatrix &dMdM);
+   /// 2nd derivative of the third invariant for the given 3x3 matrix @a M.
+   static void Dim3Invariant3_dMdM(const DenseMatrix &M, int i, int j,
+                                   DenseMatrix &dMdM);
+
    /** @brief See AssembleH(). This is a default implementation for the case
        when the 2nd derivatives of the metric are pre-computed and stored into
        @a H. This function is used in combination with AD-based computations. */
@@ -44,6 +86,8 @@ public:
    TMOP_QualityMetric() : Jtr(NULL) { }
    virtual ~TMOP_QualityMetric() { }
 
+   bool use_old_invariants_code = false;
+
    /** @brief Specify the reference-element -> target-element Jacobian matrix
        for the point of interest.
 
@@ -104,6 +148,11 @@ public:
    virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                           const real_t weight, DenseMatrix &A) const = 0;
 
+   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const DenseMatrix &dx_dt,
+                          const real_t weight, DenseMatrix &A) const
+   { MFEM_ABORT("not implemented"); }
+
    /** @brief Return the metric ID. */
    virtual int Id() const { return 0; }
 };
@@ -388,6 +437,10 @@ public:
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const DenseMatrix &dx_dt,
+                          const real_t weight, DenseMatrix &A) const override;
+
    int Id() const override { return 2; }
 };
 
@@ -2062,6 +2115,10 @@ protected:
    Array<int> surf_fit_dof_count;            // Number of dofs per node.
    Array<int> surf_fit_marker_dof_index;     // Indices of nodes to fit.
 
+   // Tangential optimization.
+   AnalyticCompositeSurface *tan_analytic_surf = nullptr;
+   FiniteElementSpace *fes_mesh = nullptr;
+
    DiscreteAdaptTC *discr_tc;
 
    // Parameters for FD-based Gradient & Hessian calculation.
@@ -2159,11 +2216,11 @@ protected:
    void ComputeNormalizationEnergies(const GridFunction &x,
                                      real_t &metric_energy, real_t &lim_energy);
 
-   void AssembleElementVectorExact(const FiniteElement &el,
+   virtual void AssembleElementVectorExact(const FiniteElement &el,
                                    ElementTransformation &T,
                                    const Vector &d_el, Vector &elvect);
 
-   void AssembleElementGradExact(const FiniteElement &el,
+   virtual void AssembleElementGradExact(const FiniteElement &el,
                                  ElementTransformation &T,
                                  const Vector &d_el, DenseMatrix &elmat);
 
@@ -2495,6 +2552,16 @@ public:
       return surf_fit_gf != NULL || surf_fit_pos != NULL;
    }
 
+   // Tangential movement.
+   void EnableTangentialMovement(AnalyticCompositeSurface &surf,
+                                 FiniteElementSpace &mesh_fes)
+   {
+      tan_analytic_surf = &surf;
+      fes_mesh = &mesh_fes;
+   }
+   const AnalyticCompositeSurface *GetAnalyticSurface() const
+   { return tan_analytic_surf; }
+
    /// Update the original/reference nodes used for limiting.
    void SetLimitingNodes(const GridFunction &n0) { lim_nodes0 = &n0; }
 

fem/tmop/AnalyticalSurface.cpp

@@ -0,0 +1,323 @@
+// Copyright (c) 2017, Lawrence LivermoreA National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#include "AnalyticalSurface.hpp"
+
+namespace mfem
+{
+
+AnalyticSurface::AnalyticSurface(const Array<bool> &marker)
+    : dof_marker(marker), dists(0)
+{
+   //geometry->ComputeDistances(coord, pmesh, pfes_mesh);
+}
+
+AnalyticCompositeSurface::AnalyticCompositeSurface
+    (const Array<AnalyticSurface *> &surf)
+    : AnalyticSurface(), surfaces(surf), dof_to_surface(0)
+{
+   UpdateDofToSurface();
+}
+
+// For each DOF, the Array dof_to_surface has its corresponding surface index.
+// -1: the DOF is not associated with any surface.
+// -2: the DOF is associated with more than one surface.
+void AnalyticCompositeSurface::UpdateDofToSurface()
+{
+   if (surfaces.Size() == 0)
+   {
+      dof_to_surface.DeleteAll();
+      return;
+   }
+
+   dof_to_surface.SetSize(surfaces[0]->dof_marker.Size());
+   dof_to_surface = -1;
+   for (int s = 0; s < surfaces.Size(); s++)
+   {
+      for (int i = 0; i < dof_to_surface.Size(); i++)
+      {
+         if (surfaces[s]->dof_marker[i] == true)
+         {
+            if (dof_to_surface[i] == -1)    { dof_to_surface[i] =  s; }
+            else if (dof_to_surface[i] >= 0) { dof_to_surface[i] = -2; }
+         }
+      }
+   }
+}
+
+const AnalyticSurface *AnalyticCompositeSurface::GetSurface(int dof_id) const
+{
+   if (surfaces.Size() == 0)       { return nullptr; }
+   if (dof_to_surface[dof_id] < 0) { return nullptr; }
+
+   return surfaces[dof_to_surface[dof_id]];
+}
+
+
+void AnalyticCompositeSurface::ConvertPhysCoordToParam(const Vector &coord_x,
+                                                       Vector &coord_t)
+{
+   coord_t = coord_x;
+   for (int s = 0; s < surfaces.Size(); s++)
+   {
+      surfaces[s]->ConvertPhysCoordToParam(coord_x, coord_t);
+   }
+}
+
+void AnalyticCompositeSurface::ConvertParamCoordToPhys(const Vector &coord_t,
+                                                       Vector &coord_x) const
+{
+   coord_x = coord_t;
+   for (int s = 0; s < surfaces.Size(); s++)
+   {
+      surfaces[s]->ConvertParamCoordToPhys(coord_t, coord_x);
+   }
+}
+
+void AnalyticCompositeSurface::ConvertParamToPhys(const Array<int> &vdofs,
+                                                  const Vector &coord_t,
+                                                  Vector &coord_x) const
+{
+   coord_x = coord_t;
+   for (int s = 0; s < surfaces.Size(); s++)
+   {
+      surfaces[s]->ConvertParamToPhys(vdofs, coord_t, coord_x);
+   }
+}
+
+void Analytic2DCurve::ConvertPhysCoordToParam(const Vector &coord_x,
+                                              Vector &coord_t)
+{
+   const int ndof = coord_x.Size() / 2;
+   double t;
+   dists.SetSize(ndof);
+   dists = 0.0;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         t_of_xy(coord_x(i), coord_x(ndof + i), dists(i), t);
+         coord_t(i)        = t;
+         coord_t(ndof + i) = 0.0;
+      }
+   }
+}
+
+void Analytic2DCurve::ConvertParamCoordToPhys(const Vector &coord_t,
+                                              Vector &coord_x) const
+{
+   const int ndof = coord_x.Size() / 2;
+   double x, y;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         xy_of_t(coord_t(i), dists(i), x, y);
+         coord_x(i)        = x;
+         coord_x(ndof + i) = y;
+      }
+   }
+}
+
+void Analytic2DCurve::ConvertParamToPhys(const Array<int> &vdofs,
+                                         const Vector &coord_t,
+                                         Vector &coord_x) const
+{
+   const int ndof = vdofs.Size() / 2;
+   double x, y;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[vdofs[i]])
+      {
+         xy_of_t(coord_t(i), dists(vdofs[i]), x, y);
+         coord_x(i)        = x;
+         coord_x(ndof + i) = y;
+      }
+   }
+}
+
+void Analytic2DCurve::Deriv_1(const double *param, DenseMatrix &deriv) const
+{
+   deriv.SetSize(2, 1);
+   deriv(0, 0) = dx_dt(param[0]);
+   deriv(1, 0) = dy_dt(param[0]);
+}
+
+void Analytic2DCurve::Deriv_2(const double *param, DenseTensor &deriv) const
+{
+   deriv.SetSize(2, 1, 1);
+   deriv(0, 0, 0) = dx_dtdt(param[0]);
+   deriv(1, 0, 0) = dy_dtdt(param[0]);
+}
+
+void Analytic3DCurve::ConvertPhysCoordToParam(const Vector &coord_x,
+                                              Vector &coord_t)
+{
+   const int ndof = coord_x.Size() / 3;
+   double t;
+   dists.SetSize(2 * ndof);
+   dists = 0.0;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         t_of_xyz(coord_x(i), coord_x(ndof + i), coord_x(2 * ndof + i),
+                  dists(i), dists(ndof + i), t);
+         coord_t(i)            = t;
+         coord_t(ndof + i)     = 0.0;
+         coord_t(2 * ndof + i) = 0.0;
+      }
+   }
+}
+
+void Analytic3DCurve::ConvertParamCoordToPhys(const Vector &coord_t,
+                                              Vector &coord_x) const
+{
+   const int ndof = coord_x.Size() / 3;
+   double x, y, z;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         xyz_of_t(coord_t(i), dists(i), dists(ndof + i),
+                  x, y, z);
+         coord_x(i)            = x;
+         coord_x(ndof + i)     = y;
+         coord_x(2 * ndof + i) = z;
+      }
+   }
+}
+
+void Analytic3DCurve::ConvertParamToPhys(const Array<int> &vdofs,
+                                         const Vector &coord_t,
+                                         Vector &coord_x) const
+{
+   const int ndof = vdofs.Size() / 3;
+   double x, y, z;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[vdofs[i]])
+      {
+         xyz_of_t(coord_t(i), dists(vdofs[i]),
+                  dists(dof_marker.Size() + vdofs[i]), x, y, z);
+         coord_x(i)            = x;
+         coord_x(ndof + i)     = y;
+         coord_x(2 * ndof + i) = z;
+      }
+   }
+}
+
+void Analytic3DCurve::Deriv_1(const double *param, DenseMatrix &deriv) const
+{
+   deriv.SetSize(3, 1);
+   deriv(0, 0) = dx_dt(param[0]);
+   deriv(1, 0) = dy_dt(param[0]);
+   deriv(2, 0) = dz_dt(param[0]);
+}
+
+void Analytic3DCurve::Deriv_2(const double *param, DenseTensor &deriv) const
+{
+   deriv.SetSize(3, 1, 1);
+   deriv(0, 0, 0) = dx_dtdt(param[0]);
+   deriv(1, 0, 0) = dy_dtdt(param[0]);
+   deriv(2, 0, 0) = dz_dtdt(param[0]);
+}
+
+void Analytic3DSurface::ConvertPhysCoordToParam(const Vector &coord_x,
+                                                Vector &coord_t)
+{
+   const int ndof = coord_x.Size() / 3;
+   double u, v;
+   dists.SetSize(ndof);
+   dists = 0.0;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         uv_of_xyz(coord_x(i), coord_x(ndof + i), coord_x(2 * ndof + i),
+                   dists(i), u, v);
+         coord_t(i)            = u;
+         coord_t(ndof + i)     = v;
+         coord_t(2 * ndof + i) = 0.0;
+      }
+   }
+}
+
+void Analytic3DSurface::ConvertParamCoordToPhys(const Vector &coord_t,
+                                                Vector &coord_x) const
+{
+   const int ndof = coord_x.Size() / 3;
+   double x, y, z;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[i])
+      {
+         xyz_of_uv(coord_t(i), coord_t(ndof + i), dists(i), x, y, z);
+         coord_x(i)            = x;
+         coord_x(ndof + i)     = y;
+         coord_x(2 * ndof + i) = z;
+      }
+   }
+}
+
+void Analytic3DSurface::ConvertParamToPhys(const Array<int> &vdofs,
+                                           const Vector &coord_t,
+                                           Vector &coord_x) const
+{
+   const int ndof = vdofs.Size() / 3;
+   double x, y, z;
+   for (int i = 0; i < ndof; i++)
+   {
+      if (dof_marker[vdofs[i]])
+      {
+         xyz_of_uv(coord_t(i), coord_t(ndof + i), dists(vdofs[i]),
+                  x, y, z);
+         coord_x(i)            = x;
+         coord_x(ndof + i)     = y;
+         coord_x(2 * ndof + i) = z;
+      }
+   }
+}
+
+void Analytic3DSurface::Deriv_1(const double *param, DenseMatrix &deriv) const
+{
+   deriv.SetSize(3, 2);
+   deriv(0, 0) = dx_du(param[0], param[1]);
+   deriv(1, 0) = dy_du(param[0], param[1]);
+   deriv(2, 0) = dz_du(param[0], param[1]);
+   deriv(0, 1) = dx_dv(param[0], param[1]);
+   deriv(1, 1) = dy_dv(param[0], param[1]);
+   deriv(2, 1) = dz_dv(param[0], param[1]);
+}
+
+void Analytic3DSurface::Deriv_2(const double *param, DenseTensor &deriv) const
+{
+   deriv.SetSize(3, 2, 2);
+   deriv(0, 0, 0) = dx_dudu(param[0], param[1]);
+   deriv(1, 0, 0) = dy_dudu(param[0], param[1]);
+   deriv(2, 0, 0) = dz_dudu(param[0], param[1]);
+   deriv(0, 0, 1) = dx_dudv(param[0], param[1]);
+   deriv(1, 0, 1) = dy_dudv(param[0], param[1]);
+   deriv(2, 0, 1) = dz_dudv(param[0], param[1]);
+   deriv(0, 1, 0) = deriv(0, 0, 1);
+   deriv(1, 1, 0) = deriv(1, 0, 1);
+   deriv(2, 1, 0) = deriv(2, 0, 1);
+   deriv(0, 1, 1) = dx_dvdv(param[0], param[1]);
+   deriv(1, 1, 1) = dy_dvdv(param[0], param[1]);
+   deriv(2, 1, 1) = dz_dvdv(param[0], param[1]);
+}
+
+} // namespace mfem

fem/tmop/AnalyticalSurface.hpp

@@ -0,0 +1,255 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#ifndef MFEM_ANALYTICAL_SURFACE
+#define MFEM_ANALYTICAL_SURFACE
+
+#include "../pgridfunc.hpp"
+
+namespace mfem
+{
+
+class AnalyticSurface
+{
+protected:
+   friend class AnalyticCompositeSurface;
+
+   Array<bool> dof_marker;
+   Vector dists;
+
+public:
+   AnalyticSurface() : dof_marker(0), dists(0) { }
+   AnalyticSurface(const Array<bool> &marker);
+   virtual ~AnalyticSurface() = default;
+
+   const Array<bool> &GetMarker() const { return dof_marker; }
+
+   // Go from physical to parametric coordinates on the whole mesh.
+   // If the point is not on the curve, computes its offset.
+   // 2D curve: (x, y)    -> t.
+   // 3D curve: (x, y, z) -> t.
+   // 3D surf:  (x, y, z) -> (u, v).
+   virtual void ConvertPhysCoordToParam(const Vector &coord_x,
+                                        Vector &coord_t) = 0;
+
+   // Go from parametric to physical coordinates on the whole mesh:
+   // 2D curve: t      -> (x, y).
+   // 3D curve: t      -> (x, y, z).
+   // 3D surf:  (u, v) -> (x, y, z).
+   virtual void ConvertParamCoordToPhys(const Vector &coord_t,
+                                        Vector &coord_x) const = 0;
+
+   // Go from parametric to physical coordinates on a single element with given
+   // vdofs.
+   virtual void ConvertParamToPhys(const Array<int> &vdofs,
+                                   const Vector &coord_t,
+                                   Vector &coord_x) const = 0;
+
+   /// Number of free parameters stored for each constrained node.
+   virtual int NumParams() const = 0;
+
+   // First derivatives:
+   // 2D curve:  t     -> (dx/dt, dy/dt).
+   // 3D curve:  t     -> (dx/dt, dy/dt, dz/dt).
+   // 3D surf:  (u, v) -> (dx/du, dy/du, dz/du, dx/dv, dy/dv, dz/dv).
+   virtual void Deriv_1(const double *param, DenseMatrix &deriv) const = 0;
+
+   // Second derivatives:
+   // 2D curve:  t     -> (dx_dtdt, dy_dtdt).
+   // 3D curve:  t     -> (dx_dtdt, dy_dtdt, dz_dtdt).
+   // 3D surf:  (u, v) -> (dx/dudu, dx/dudv, dx/dvdv, ...).
+   virtual void Deriv_2(const double *param, DenseTensor &deriv) const = 0;
+};
+
+class AnalyticCompositeSurface : public AnalyticSurface
+{
+protected:
+   const Array<AnalyticSurface *> &surfaces;
+   Array<int> dof_to_surface;
+
+public:
+   AnalyticCompositeSurface(const Array<AnalyticSurface *> &surf);
+
+   int GetNumSurfaces() const { return surfaces.Size(); }
+
+   /// Must be called after the Array of surfaces is changed.
+   void UpdateDofToSurface();
+
+   const AnalyticSurface *GetSurfaceID(int s_id) const
+   {
+      MFEM_VERIFY(s_id < surfaces.Size(), "wrong index");
+      return surfaces[s_id];
+   }
+
+   /// Surface corresponding to dof_id.
+   const AnalyticSurface *GetSurface(int dof_id) const;
+
+   void ConvertPhysCoordToParam(const Vector &coord_x,
+                                Vector &coord_t) override;
+
+   void ConvertParamCoordToPhys(const Vector &coord_t,
+                                Vector &coord_x) const override;
+
+   void ConvertParamToPhys(const Array<int> &vdofs,
+                           const Vector &coord_t,
+                           Vector &coord_x) const override;
+
+   int NumParams() const override
+   { return 0; }
+
+   void Deriv_1(const double *param, DenseMatrix &deriv) const override
+   { MFEM_ABORT("Use GetSurface(surface_id)->Deriv_1(...);") }
+
+   void Deriv_2(const double *param, DenseTensor &deriv) const override
+   { MFEM_ABORT("Use GetSurface(surface_id)->Deriv_2(...);") }
+};
+
+class Analytic2DCurve : public AnalyticSurface
+{
+public:
+   Analytic2DCurve(const Array<bool> &marker)
+    : AnalyticSurface(marker) { }
+
+   // (x, y) -> t on the whole curve.
+   void ConvertPhysCoordToParam(const Vector &coord_x,
+                                Vector &coord_t) override;
+
+   // t -> (x, y) on the whole curve.
+   void ConvertParamCoordToPhys(const Vector &coord_t,
+                                Vector &coord_x) const override;
+
+   // t -> (x, y) on a single element given by vdofs.
+   void ConvertParamToPhys(const Array<int> &vdofs,
+                           const Vector &coord_t,
+                           Vector &coord_x) const override;
+
+   int NumParams() const override { return 1; }
+
+   // t -> (dx_dt, dy_dt).
+   void Deriv_1(const double *param, DenseMatrix &deriv) const override;
+
+   // t -> (dx_dtdt, dy_dtdt).
+   void Deriv_2(const double *param, DenseTensor &deriv) const override;
+
+   virtual void t_of_xy(double x, double y,
+                        double &dist, double &t) const = 0;
+   virtual void xy_of_t(double t, double dist,
+                        double &x, double &y) const = 0;
+
+   virtual double dx_dt(double t) const = 0;
+   virtual double dy_dt(double t) const = 0;
+
+   virtual double dx_dtdt(double t) const = 0;
+   virtual double dy_dtdt(double t) const = 0;
+};
+
+class Analytic3DCurve : public AnalyticSurface
+{
+public:
+   Analytic3DCurve(const Array<bool> &marker)
+    : AnalyticSurface(marker) { }
+
+   // (x,y,z) -> t on the whole curve.
+   void ConvertPhysCoordToParam(const Vector &coord_x,
+                                Vector &coord_t) override;
+
+   // t -> (x,y,z) on the whole curve.
+   void ConvertParamCoordToPhys(const Vector &coord_t,
+                                Vector &coord_x) const override;
+
+   // t -> (x,y,z) on a single element given by vdofs.
+   void ConvertParamToPhys(const Array<int> &vdofs,
+                           const Vector &coord_t,
+                           Vector &coord_x) const override;
+
+   int NumParams() const override { return 1; }
+
+   // t -> (dx_dt, dy_dt, dz_dt).
+   void Deriv_1(const double *param, DenseMatrix &deriv) const override;
+
+   // t -> (dx_dtdt, dy_dtdt, dz_dtdt).
+   void Deriv_2(const double *param, DenseTensor &deriv) const override;
+
+   virtual void t_of_xyz(double x, double y, double z,
+                         double &dist1, double &dist2, double &t) const = 0;
+   virtual void xyz_of_t(double t, double dist1, double dist2,
+                         double &x, double &y, double &z) const = 0;
+
+   virtual double dx_dt(double t) const = 0;
+   virtual double dy_dt(double t) const = 0;
+   virtual double dz_dt(double t) const = 0;
+
+   virtual double dx_dtdt(double t) const = 0;
+   virtual double dy_dtdt(double t) const = 0;
+   virtual double dz_dtdt(double t) const = 0;
+};
+
+class Analytic3DSurface : public AnalyticSurface
+{
+public:
+   Analytic3DSurface(const Array<bool> &marker)
+    : AnalyticSurface(marker) { }
+
+   // (x,y,z) -> (u,v) on whole surface.
+   void ConvertPhysCoordToParam(const Vector &coord_x,
+                                Vector &coord_t) override;
+
+   // (u,v) -> (x,y,z) on whole surface.
+   void ConvertParamCoordToPhys(const Vector &coord_t,
+                                Vector &coord_x) const override;
+
+   // (u,v) -> (x,y,z) on a single element given by vdofs.
+   void ConvertParamToPhys(const Array<int> &vdofs,
+                           const Vector &coord_t,
+                           Vector &coord_x) const override;
+
+   int NumParams() const override { return 2; }
+
+   // (u,v) -> (dx_du, dx_dv
+   //           dy_du, dy_dv,
+   //           dz_du, dz_dv)
+   void Deriv_1(const double *param, DenseMatrix &deriv) const override;
+
+   // (u,v) -> (dx_dudu, dx_dudv    (dx_dudv, dx_dvdv,
+   //           dy_dudu, dy_dudv,    dy_dudv, dy_dvdv,
+   //           dz_dudu, dz_dudv),   dz_dudv, dz_dvdv)
+   void Deriv_2(const double *param, DenseTensor &deriv) const override;
+
+   virtual void uv_of_xyz(double x, double y, double z,
+                          double &dist, double &u, double &v) const = 0;
+   virtual void xyz_of_uv(double u, double v, double dist,
+                          double &x, double &y, double &z) const = 0;
+
+   virtual double dx_du(double u, double v) const = 0;
+   virtual double dy_du(double u, double v) const = 0;
+   virtual double dz_du(double u, double v) const = 0;
+   virtual double dx_dv(double u, double v) const = 0;
+   virtual double dy_dv(double u, double v) const = 0;
+   virtual double dz_dv(double u, double v) const = 0;
+
+   virtual double dx_dudu(double u, double v) const = 0;
+   virtual double dy_dudu(double u, double v) const = 0;
+   virtual double dz_dudu(double u, double v) const = 0;
+   virtual double dx_dudv(double u, double v) const = 0;
+   virtual double dy_dudv(double u, double v) const = 0;
+   virtual double dz_dudv(double u, double v) const = 0;
+   virtual double dx_dvdv(double u, double v) const = 0;
+   virtual double dy_dvdv(double u, double v) const = 0;
+   virtual double dz_dvdv(double u, double v) const = 0;
+};
+
+}
+#endif

fem/tmop_tools.cpp

@@ -696,6 +696,13 @@ void TMOPNewtonSolver::Mult(const Vector &b, Vector &x) const
       if (co) { co->SetInitialMeshPos(&x_0); }
    }
 
+   // Check if there's a surface parametrization.
+   auto integ = dynamic_cast<TMOP_Integrator *>(integs[0]);
+   parametric = (integ->GetAnalyticSurface() == nullptr) ? false : true;
+
+   MFEM_VERIFY((periodic && parametric) == false,
+               "Periodic meshes + parametrization is not supported");
+
    // Solve for the displacement, which always starts from zero.
    Vector dx(height); dx = 0.0;
    if (solver_type == 0)      { NewtonSolver::Mult(b, dx); }
@@ -974,12 +981,22 @@ real_t TMOPNewtonSolver::ComputeMinDet(const Vector &d_loc,
 {
    real_t min_detJ = infinity();
    const int NE = fes.GetNE(), dim = fes.GetMesh()->Dimension();
-   Array<int> xdofs;
-   DenseMatrix Jpr(dim);
    const bool mixed_mesh = fes.GetMesh()->GetNumGeometries(dim) > 1;
-   if (dim == 1 || mixed_mesh ||
+   if (1 || dim == 1 || mixed_mesh ||
        UsesTensorBasis(fes) == false || fes.IsVariableOrder())
    {
+      const AnalyticCompositeSurface *surf = nullptr;
+      if (parametric)
+      {
+         const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
+         const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
+         auto ti = dynamic_cast<TMOP_Integrator *>(integs[0]);
+         surf = ti->GetAnalyticSurface();
+      }
+
+      Array<int> xdofs;
+      DenseMatrix Jpr(dim);
+
       for (int i = 0; i < NE; i++)
       {
          const int dof = fes.GetFE(i)->GetDof();
@@ -998,6 +1015,13 @@ real_t TMOPNewtonSolver::ComputeMinDet(const Vector &d_loc,
          d_loc.GetSubVector(xdofs, d_loc_el);
          posV += d_loc_el;
 
+         Vector convertedX(posV.Size());
+         if (parametric)
+         {
+            surf->ConvertParamToPhys(xdofs, posV, convertedX);
+            posV = convertedX;
+         }
+
          const IntegrationRule &irule = GetIntegrationRule(*fes.GetFE(i));
          const int nsp = irule.GetNPoints();
          for (int j = 0; j < nsp; j++)

fem/tmop_tools.hpp

@@ -178,6 +178,9 @@ protected:
    mutable GridFunction x_0;
    mutable bool periodic = false;
 
+   // Is there a parametrization (e.g. x(t), y(t) for 2D) in the TMOP problem.
+   mutable bool parametric = false;
+
    // Line search step is rejected if min(detJ) <= min_detJ_limit.
    real_t min_detJ_limit = 0.0;
 

miniapps/meshing/CMakeLists.txt

@@ -14,6 +14,7 @@ set(MESH_FILES
   blade.mesh
   cube.mesh
   cube-tet.mesh
+  cube01.mesh
   icf.mesh
   jagged.mesh
   square01.mesh
@@ -151,6 +152,12 @@ if (MFEM_USE_MPI)
     MAIN phpref.cpp
     LIBRARIES mfem)
 
+  add_mfem_miniapp(tmop-tangent
+    MAIN tmop-tangent.cpp
+    ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem-common)
+  add_dependencies(tmop-tangent copy_miniapps_meshing_data)
+
   # Add parallel tests.
   if (MFEM_ENABLE_TESTING)
     set(PARALLEL_TESTS
@@ -160,6 +167,7 @@ if (MFEM_USE_MPI)
       fit-node-position
       pminimal-surface
       phpref
+      tmop-tangent
     )
     # Meshing miniapps that return MFEM_SKIP_RETURN_VALUE in some cases:
     set(SKIP_TESTS)

miniapps/meshing/cube01.mesh

@@ -0,0 +1,312 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+64
+1 5 0 1 6 5 25 26 31 30
+1 5 25 26 31 30 50 51 56 55
+1 5 26 27 32 31 51 52 57 56
+1 5 1 2 7 6 26 27 32 31
+1 5 6 7 12 11 31 32 37 36
+1 5 31 32 37 36 56 57 62 61
+1 5 30 31 36 35 55 56 61 60
+1 5 5 6 11 10 30 31 36 35
+1 5 10 11 16 15 35 36 41 40
+1 5 11 12 17 16 36 37 42 41
+1 5 16 17 22 21 41 42 47 46
+1 5 15 16 21 20 40 41 46 45
+1 5 40 41 46 45 65 66 71 70
+1 5 41 42 47 46 66 67 72 71
+1 5 36 37 42 41 61 62 67 66
+1 5 35 36 41 40 60 61 66 65
+1 5 60 61 66 65 85 86 91 90
+1 5 61 62 67 66 86 87 92 91
+1 5 66 67 72 71 91 92 97 96
+1 5 65 66 71 70 90 91 96 95
+1 5 90 91 96 95 115 116 121 120
+1 5 91 92 97 96 116 117 122 121
+1 5 86 87 92 91 111 112 117 116
+1 5 85 86 91 90 110 111 116 115
+1 5 80 81 86 85 105 106 111 110
+1 5 75 76 81 80 100 101 106 105
+1 5 50 51 56 55 75 76 81 80
+1 5 55 56 61 60 80 81 86 85
+1 5 56 57 62 61 81 82 87 86
+1 5 51 52 57 56 76 77 82 81
+1 5 76 77 82 81 101 102 107 106
+1 5 81 82 87 86 106 107 112 111
+1 5 82 83 88 87 107 108 113 112
+1 5 77 78 83 82 102 103 108 107
+1 5 52 53 58 57 77 78 83 82
+1 5 57 58 63 62 82 83 88 87
+1 5 58 59 64 63 83 84 89 88
+1 5 53 54 59 58 78 79 84 83
+1 5 78 79 84 83 103 104 109 108
+1 5 83 84 89 88 108 109 114 113
+1 5 88 89 94 93 113 114 119 118
+1 5 87 88 93 92 112 113 118 117
+1 5 92 93 98 97 117 118 123 122
+1 5 93 94 99 98 118 119 124 123
+1 5 68 69 74 73 93 94 99 98
+1 5 67 68 73 72 92 93 98 97
+1 5 62 63 68 67 87 88 93 92
+1 5 63 64 69 68 88 89 94 93
+1 5 38 39 44 43 63 64 69 68
+1 5 37 38 43 42 62 63 68 67
+1 5 42 43 48 47 67 68 73 72
+1 5 43 44 49 48 68 69 74 73
+1 5 18 19 24 23 43 44 49 48
+1 5 17 18 23 22 42 43 48 47
+1 5 12 13 18 17 37 38 43 42
+1 5 13 14 19 18 38 39 44 43
+1 5 8 9 14 13 33 34 39 38
+1 5 33 34 39 38 58 59 64 63
+1 5 32 33 38 37 57 58 63 62
+1 5 7 8 13 12 32 33 38 37
+1 5 2 3 8 7 27 28 33 32
+1 5 27 28 33 32 52 53 58 57
+1 5 28 29 34 33 53 54 59 58
+1 5 3 4 9 8 28 29 34 33
+
+boundary
+96
+1 3 0 5 6 1
+1 3 1 6 7 2
+1 3 2 7 8 3
+1 3 3 8 9 4
+1 3 5 10 11 6
+1 3 6 11 12 7
+1 3 7 12 13 8
+1 3 8 13 14 9
+1 3 10 15 16 11
+1 3 11 16 17 12
+1 3 12 17 18 13
+1 3 13 18 19 14
+1 3 15 20 21 16
+1 3 16 21 22 17
+1 3 17 22 23 18
+1 3 18 23 24 19
+6 3 100 101 106 105
+6 3 101 102 107 106
+6 3 102 103 108 107
+6 3 103 104 109 108
+6 3 105 106 111 110
+6 3 106 107 112 111
+6 3 107 108 113 112
+6 3 108 109 114 113
+6 3 110 111 116 115
+6 3 111 112 117 116
+6 3 112 113 118 117
+6 3 113 114 119 118
+6 3 115 116 121 120
+6 3 116 117 122 121
+6 3 117 118 123 122
+6 3 118 119 124 123
+5 3 0 25 30 5
+5 3 5 30 35 10
+5 3 10 35 40 15
+5 3 15 40 45 20
+5 3 25 50 55 30
+5 3 30 55 60 35
+5 3 35 60 65 40
+5 3 40 65 70 45
+5 3 50 75 80 55
+5 3 55 80 85 60
+5 3 60 85 90 65
+5 3 65 90 95 70
+5 3 75 100 105 80
+5 3 80 105 110 85
+5 3 85 110 115 90
+5 3 90 115 120 95
+3 3 4 9 34 29
+3 3 9 14 39 34
+3 3 14 19 44 39
+3 3 19 24 49 44
+3 3 29 34 59 54
+3 3 34 39 64 59
+3 3 39 44 69 64
+3 3 44 49 74 69
+3 3 54 59 84 79
+3 3 59 64 89 84
+3 3 64 69 94 89
+3 3 69 74 99 94
+3 3 79 84 109 104
+3 3 84 89 114 109
+3 3 89 94 119 114
+3 3 94 99 124 119
+2 3 0 1 26 25
+2 3 25 26 51 50
+2 3 50 51 76 75
+2 3 75 76 101 100
+2 3 1 2 27 26
+2 3 26 27 52 51
+2 3 51 52 77 76
+2 3 76 77 102 101
+2 3 2 3 28 27
+2 3 27 28 53 52
+2 3 52 53 78 77
+2 3 77 78 103 102
+2 3 3 4 29 28
+2 3 28 29 54 53
+2 3 53 54 79 78
+2 3 78 79 104 103
+4 3 20 45 46 21
+4 3 45 70 71 46
+4 3 70 95 96 71
+4 3 95 120 121 96
+4 3 21 46 47 22
+4 3 46 71 72 47
+4 3 71 96 97 72
+4 3 96 121 122 97
+4 3 22 47 48 23
+4 3 47 72 73 48
+4 3 72 97 98 73
+4 3 97 122 123 98
+4 3 23 48 49 24
+4 3 48 73 74 49
+4 3 73 98 99 74
+4 3 98 123 124 99
+
+vertices
+125
+3
+0 0 0
+0.25 0 0
+0.5 0 0
+0.75 0 0
+1 0 0
+0 0.25 0
+0.25 0.25 0
+0.5 0.25 0
+0.75 0.25 0
+1 0.25 0
+0 0.5 0
+0.25 0.5 0
+0.5 0.5 0
+0.75 0.5 0
+1 0.5 0
+0 0.75 0
+0.25 0.75 0
+0.5 0.75 0
+0.75 0.75 0
+1 0.75 0
+0 1 0
+0.25 1 0
+0.5 1 0
+0.75 1 0
+1 1 0
+0 0 0.25
+0.25 0 0.25
+0.5 0 0.25
+0.75 0 0.25
+1 0 0.25
+0 0.25 0.25
+0.25 0.25 0.25
+0.5 0.25 0.25
+0.75 0.25 0.25
+1 0.25 0.25
+0 0.5 0.25
+0.25 0.5 0.25
+0.5 0.5 0.25
+0.75 0.5 0.25
+1 0.5 0.25
+0 0.75 0.25
+0.25 0.75 0.25
+0.5 0.75 0.25
+0.75 0.75 0.25
+1 0.75 0.25
+0 1 0.25
+0.25 1 0.25
+0.5 1 0.25
+0.75 1 0.25
+1 1 0.25
+0 0 0.5
+0.25 0 0.5
+0.5 0 0.5
+0.75 0 0.5
+1 0 0.5
+0 0.25 0.5
+0.25 0.25 0.5
+0.5 0.25 0.5
+0.75 0.25 0.5
+1 0.25 0.5
+0 0.5 0.5
+0.25 0.5 0.5
+0.5 0.5 0.5
+0.75 0.5 0.5
+1 0.5 0.5
+0 0.75 0.5
+0.25 0.75 0.5
+0.5 0.75 0.5
+0.75 0.75 0.5
+1 0.75 0.5
+0 1 0.5
+0.25 1 0.5
+0.5 1 0.5
+0.75 1 0.5
+1 1 0.5
+0 0 0.75
+0.25 0 0.75
+0.5 0 0.75
+0.75 0 0.75
+1 0 0.75
+0 0.25 0.75
+0.25 0.25 0.75
+0.5 0.25 0.75
+0.75 0.25 0.75
+1 0.25 0.75
+0 0.5 0.75
+0.25 0.5 0.75
+0.5 0.5 0.75
+0.75 0.5 0.75
+1 0.5 0.75
+0 0.75 0.75
+0.25 0.75 0.75
+0.5 0.75 0.75
+0.75 0.75 0.75
+1 0.75 0.75
+0 1 0.75
+0.25 1 0.75
+0.5 1 0.75
+0.75 1 0.75
+1 1 0.75
+0 0 1
+0.25 0 1
+0.5 0 1
+0.75 0 1
+1 0 1
+0 0.25 1
+0.25 0.25 1
+0.5 0.25 1
+0.75 0.25 1
+1 0.25 1
+0 0.5 1
+0.25 0.5 1
+0.5 0.5 1
+0.75 0.5 1
+1 0.5 1
+0 0.75 1
+0.25 0.75 1
+0.5 0.75 1
+0.75 0.75 1
+1 0.75 1
+0 1 1
+0.25 1 1
+0.5 1 1
+0.75 1 1
+1 1 1

miniapps/meshing/makefile

@@ -24,7 +24,7 @@ SEQ_MINIAPPS = mobius-strip klein-bottle toroid trimmer twist mesh-explorer\
 	shaper extruder mesh-optimizer minimal-surface polar-nc reflector\
 	ref321 mesh-quality hpref
 PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting fit-node-position\
-	phpref mesh-bounding-boxes
+	phpref mesh-bounding-boxes tmop-tangent
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -57,8 +57,9 @@ lib-common:
 # Rules to copy the *.mesh files - needed for running the sample runs when
 # building out-of-source:
 ifneq ($(SRC),)
-MESH_FILES = amr-quad-q2.mesh blade.mesh cube.mesh cube-tet.mesh icf.mesh\
-jagged.mesh square01.mesh square01-tri.mesh stretched2D.mesh stretched3D.mesh
+MESH_FILES = amr-quad-q2.mesh blade.mesh cube.mesh cube-tet.mesh cube01.mesh\
+icf.mesh jagged.mesh square01.mesh square01-tri.mesh stretched2D.mesh\
+stretched3D.mesh
 $(MESH_FILES): %: $(SRC)%
 	ln -sf $(<) .
 mesh-optimizer pmesh-optimizer pmesh-fitting fit-node-position: | $(MESH_FILES)
@@ -123,6 +124,7 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
+	@rm -f tmop-tangent.mesh tmop-tangent.gf
 	@rm -f mobius-strip.mesh klein-bottle.mesh mesh-explorer.mesh*
 	@rm -f toroid-*.mesh twist-*.mesh trimmer.mesh reflected.mesh
 	@rm -f partitioning.txt shaper.mesh extruder.mesh ref321.mesh sol.gf

miniapps/meshing/rectangle01.mesh

@@ -0,0 +1,114 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+16
+1 3 0 1 6 5
+1 3 1 2 7 6
+1 3 2 3 8 7
+1 3 3 4 9 8
+1 3 5 6 11 10
+1 3 6 7 12 11
+1 3 7 8 13 12
+1 3 8 9 14 13
+1 3 10 11 16 15
+1 3 11 12 17 16
+1 3 12 13 18 17
+1 3 13 14 19 18
+1 3 15 16 21 20
+1 3 16 17 22 21
+1 3 17 18 23 22
+1 3 18 19 24 23
+
+boundary
+16
+3 1 0 1
+3 1 1 2
+3 1 2 3
+3 1 3 4
+1 1 21 20
+1 1 22 21
+1 1 23 22
+1 1 24 23
+4 1 5 0
+4 1 10 5
+4 1 15 10
+4 1 20 15
+2 1 4 9
+2 1 9 14
+2 1 14 19
+2 1 19 24
+
+vertices
+25
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: H1_2D_P1
+VDim: 2
+Ordering: 0
+
+0
+0.5
+1.0
+1.5
+2
+0
+0.5
+1.0
+1.5
+2
+0
+0.5
+1.0
+1.5
+2
+0
+0.5
+1.0
+1.5
+2
+0
+0.5
+1.0
+1.5
+2
+0
+0
+0
+0
+0
+0.25
+0.25
+0.25
+0.25
+0.25
+0.5
+0.5
+0.5
+0.5
+0.5
+0.75
+0.75
+0.75
+0.75
+0.75
+1
+1
+1
+1
+1

miniapps/meshing/square01.mesh

@@ -36,22 +36,22 @@ elements
 
 boundary
 16
-2 1 0 1
-2 1 1 2
-2 1 2 3
-2 1 3 4
-2 1 21 20
-2 1 22 21
-2 1 23 22
-2 1 24 23
-1 1 5 0
-1 1 10 5
-1 1 15 10
-1 1 20 15
-1 1 4 9
-1 1 9 14
-1 1 14 19
-1 1 19 24
+3 1 0 1
+3 1 1 2
+3 1 2 3
+3 1 3 4
+1 1 21 20
+1 1 22 21
+1 1 23 22
+1 1 24 23
+4 1 5 0
+4 1 10 5
+4 1 15 10
+4 1 20 15
+2 1 4 9
+2 1 9 14
+2 1 14 19
+2 1 19 24
 
 vertices
 25

miniapps/meshing/tmop-tangent.cpp

@@ -0,0 +1,489 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//    -------------------------------------------------------------------
+//    Tangential relaxation of boundary nodes to a given analytic surface
+//    -------------------------------------------------------------------
+//
+//
+// Sample runs:
+// mpirun -np 4 tmop-tangent -rs 2 -m square01.mesh -o 2 -qo 6
+// mpirun -np 4 tmop-tangent -rs 2 -m cube01.mesh -o 2 -qo 6
+
+#include "../common/mfem-common.hpp"
+#include "tmop-tangent.hpp"
+#include <memory>
+#include <vector>
+
+using namespace mfem;
+using namespace std;
+
+char vishost[] = "localhost";
+int  visport   = 19916;
+int  wsize     = 350;
+
+int main (int argc, char *argv[])
+{
+   // Initialize MPI.
+   Mpi::Init();
+   int myid = Mpi::WorldRank();
+   const char *mesh_file = "square01.mesh";
+   int rs_levels     = 1;
+   int mesh_poly_deg = 2;
+   int quad_order    = 5;
+   real_t d          = 100;
+   bool glvis        = true;
+   // Parse command-line options.
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+		 "Mesh file to use.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+       "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&mesh_poly_deg, "-o", "--order",
+       "Polynomial degree of mesh finite element space.");
+   args.AddOption(&quad_order, "-qo", "--quad_order",
+		 "Order of the quadrature rule.");
+   args.AddOption(&d, "-dist", "--distance",
+                  "Physical distance for limiting.");
+   args.AddOption(&glvis, "-vis", "--visualization", "-no-vis",
+		 "--no-visualization",
+		 "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (myid == 0) { args.PrintOptions(cout); }
+
+   // Read and refine the mesh.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
+   for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   const int dim = pmesh.Dimension();
+
+   // Setup mesh curvature and GridFunction that stores the coordinates.
+   H1_FECollection fec_mesh(mesh_poly_deg, dim);
+   ParFiniteElementSpace pfes_mesh(&pmesh, &fec_mesh, dim);
+   pmesh.SetNodalFESpace(&pfes_mesh);
+   ParGridFunction coord_x(&pfes_mesh), coord_t(&pfes_mesh);
+   pmesh.SetNodalGridFunction(&coord_x);
+
+   // Move the mesh nodes to have non-trivial problem.
+   const int N = coord_x.Size() / dim;
+   double a, b, c;
+   for (int i = 0; i < N; i++)
+   {
+      double x = coord_x(i);
+      double y = coord_x(i + N);
+
+      // Displace all x and y, so that the spacing is non-uniform.
+      x = x + x * (1 - x) * 0.4;
+      y = y + y * (1 - y) * 0.4;
+
+      if (dim == 2)
+      {
+         // a adds deformation inside.
+         // b pulls the top-right corner out.
+         // c adds boundary deformation.
+         // a = 0.0, b = 0.5, c = 0.0; // linear.
+         a = 0.2, b = 0.5, c = 1.5; // curved.
+         coord_x(i)     = x + a * sin(0.5 * M_PI * x) * sin(c * M_PI * y)   + b * x * y;
+         coord_x(i + N) = y + a * sin(c * M_PI * x)   * sin(0.5 * M_PI * y) + b * x * y;
+      }
+
+      else if (dim == 3)
+      {
+         double z = coord_x(i + 2*N);
+         z = z + z * (1 - z) * 0.4;
+
+         // a, b, and c have similar actions to 2D case.
+         a = 0.2, b = 0.1, c = 1.1;
+
+         coord_x(i)       = x + a * sin(0.5 * M_PI * x) * sin(c * M_PI * y) * sin(c * M_PI * z) + b * x * y * z;
+         coord_x(i + N)   = y + a * sin(c * M_PI * x) * sin(0.5 * M_PI * y) * sin(c * M_PI * z) + b * x * y * z;
+         coord_x(i + 2*N) = z + a * sin(c * M_PI * x) * sin(c * M_PI * y) * sin(0.5 * M_PI * z) + b * x * y * z;
+      }
+      
+   }
+
+   ParGridFunction x0(coord_x);
+
+   // Compute the minimum det(J) of the starting mesh.
+   double min_detJ = infinity();
+   const int NE = pmesh.GetNE();
+   for (int e = 0; e < NE; e++)
+   {
+      const IntegrationRule &ir =
+          IntRules.Get(pfes_mesh.GetFE(0)->GetGeomType(), quad_order);
+      ElementTransformation *transf = pmesh.GetElementTransformation(e);
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         transf->SetIntPoint(&ir.IntPoint(j));
+         min_detJ = min(min_detJ, transf->Jacobian().Det());
+      }
+   }
+   MPI_Allreduce(MPI_IN_PLACE, &min_detJ, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+   if (myid == 0)
+   { cout << "Minimum det(J) of the original mesh is " << min_detJ << endl; }
+   MFEM_VERIFY(min_detJ > 0.0, "Inverted initial meshes are not supported.");
+
+   // Mark which nodes to move tangentially.
+   std::vector<Array<bool>> fit_markers;
+   Array<int> ess_vdofs, ess_vdofs_marker;
+   ess_vdofs_marker.SetSize(pfes_mesh.GetVSize());
+   ess_vdofs_marker = 0;
+   ParFiniteElementSpace pfes_scalar(&pmesh, &fec_mesh, 1);
+   ParGridFunction fit_marker_vis_gf(&pfes_scalar);
+   fit_marker_vis_gf = 0.0;
+
+   if (dim == 2)
+   {
+      enum FitMarkerId
+      {
+         TOP_MARKER = 0,
+         RIGHT_MARKER,
+         BOTTOM_MARKER,
+         LEFT_MARKER,
+         FIXED_MARKER
+      };
+      fit_markers.reserve(5);
+      for (int i = 0; i < 5; i++)
+      {
+         fit_markers.emplace_back(pfes_mesh.GetNDofs());
+         fit_markers.back() = false;
+      }
+      Array<bool> &fit_marker_top = fit_markers[TOP_MARKER];
+      Array<bool> &fit_marker_right = fit_markers[RIGHT_MARKER];
+      Array<bool> &fit_marker_bottom = fit_markers[BOTTOM_MARKER];
+      Array<bool> &fit_marker_left = fit_markers[LEFT_MARKER];
+      Array<bool> &fit_marker_fixed = fit_markers[FIXED_MARKER];
+
+      Array<int> top_bdr(pmesh.bdr_attributes.Max());
+      Array<int> right_bdr(pmesh.bdr_attributes.Max());
+      Array<int> bottom_bdr(pmesh.bdr_attributes.Max());
+      Array<int> left_bdr(pmesh.bdr_attributes.Max());
+      top_bdr = 0;
+      right_bdr = 0;
+      bottom_bdr = 0;
+      left_bdr = 0;
+
+      // For the square meshes used here, boundary attributes map to sides as
+      // follows: 1 -> top, 2 -> right, 3 -> bottom, 4 -> left.
+      top_bdr[0] = 1;
+      right_bdr[1] = 1;
+      bottom_bdr[2] = 1;
+      left_bdr[3] = 1;
+
+      Array<int> top_x_marker, right_x_marker, bottom_x_marker, left_x_marker;
+      pfes_mesh.GetEssentialVDofs(top_bdr, top_x_marker, 0);
+      pfes_mesh.GetEssentialVDofs(right_bdr, right_x_marker, 0);
+      pfes_mesh.GetEssentialVDofs(bottom_bdr, bottom_x_marker, 0);
+      pfes_mesh.GetEssentialVDofs(left_bdr, left_x_marker, 0);
+
+      for (int dof = 0; dof < pfes_mesh.GetNDofs(); dof++)
+      {
+         const int vdof = pfes_mesh.DofToVDof(dof, 0);
+         const bool on_top = top_x_marker[vdof];
+         const bool on_right = right_x_marker[vdof];
+         const bool on_bottom = bottom_x_marker[vdof];
+         const bool on_left = left_x_marker[vdof];
+         const int boundary_count =
+            int(on_top) + int(on_right) + int(on_bottom) + int(on_left);
+
+         if (boundary_count == 2)
+         {
+            fit_marker_fixed[dof] = true;
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 0)] = 1;
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 1)] = 1;
+         }
+         else if (boundary_count == 1)
+         {
+            if (on_top) { fit_marker_top[dof] = true; }
+            if (on_right) { fit_marker_right[dof] = true; }
+            if (on_bottom) { fit_marker_bottom[dof] = true; }
+            if (on_left) { fit_marker_left[dof] = true; }
+
+            // The parameter value is stored in the x-component of coord_t for
+            // all 2D analytic curves, so the y-component is constrained.
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 1)] = 1;
+         }
+      }
+   }
+
+   else if (dim == 3)
+   {
+      enum FitMarkerId
+      {
+         FACE_X0_MARKER = 0,
+         FACE_X1_MARKER,
+         FACE_Y0_MARKER,
+         FACE_Y1_MARKER,
+         FACE_Z0_MARKER,
+         FACE_Z1_MARKER,
+         EDGE_X0Y0_MARKER,
+         EDGE_X0Y1_MARKER,
+         EDGE_X1Y0_MARKER,
+         EDGE_X1Y1_MARKER,
+         EDGE_X0Z0_MARKER,
+         EDGE_X0Z1_MARKER,
+         EDGE_X1Z0_MARKER,
+         EDGE_X1Z1_MARKER,
+         EDGE_Y0Z0_MARKER,
+         EDGE_Y0Z1_MARKER,
+         EDGE_Y1Z0_MARKER,
+         EDGE_Y1Z1_MARKER,
+         FIXED_MARKER
+      };
+
+      fit_markers.reserve(19);
+      for (int i = 0; i < 19; i++)
+      {
+         fit_markers.emplace_back(pfes_mesh.GetNDofs());
+         fit_markers.back() = false;
+      }
+
+      // x0 <--> face aligned with plane x = 0.
+      // x0y0 <--> edge between x0 and y0.
+      Array<bool> &fit_marker_face_x0 = fit_markers[FACE_X0_MARKER];
+      Array<bool> &fit_marker_face_x1 = fit_markers[FACE_X1_MARKER];
+      Array<bool> &fit_marker_face_y0 = fit_markers[FACE_Y0_MARKER];
+      Array<bool> &fit_marker_face_y1 = fit_markers[FACE_Y1_MARKER];
+      Array<bool> &fit_marker_face_z0 = fit_markers[FACE_Z0_MARKER];
+      Array<bool> &fit_marker_face_z1 = fit_markers[FACE_Z1_MARKER];
+      Array<bool> &fit_marker_edge_x0y0 = fit_markers[EDGE_X0Y0_MARKER];
+      Array<bool> &fit_marker_edge_x0y1 = fit_markers[EDGE_X0Y1_MARKER];
+      Array<bool> &fit_marker_edge_x1y0 = fit_markers[EDGE_X1Y0_MARKER];
+      Array<bool> &fit_marker_edge_x1y1 = fit_markers[EDGE_X1Y1_MARKER];
+      Array<bool> &fit_marker_edge_x0z0 = fit_markers[EDGE_X0Z0_MARKER];
+      Array<bool> &fit_marker_edge_x0z1 = fit_markers[EDGE_X0Z1_MARKER];
+      Array<bool> &fit_marker_edge_x1z0 = fit_markers[EDGE_X1Z0_MARKER];
+      Array<bool> &fit_marker_edge_x1z1 = fit_markers[EDGE_X1Z1_MARKER];
+      Array<bool> &fit_marker_edge_y0z0 = fit_markers[EDGE_Y0Z0_MARKER];
+      Array<bool> &fit_marker_edge_y0z1 = fit_markers[EDGE_Y0Z1_MARKER];
+      Array<bool> &fit_marker_edge_y1z0 = fit_markers[EDGE_Y1Z0_MARKER];
+      Array<bool> &fit_marker_edge_y1z1 = fit_markers[EDGE_Y1Z1_MARKER];
+      Array<bool> &fit_marker_fixed = fit_markers[FIXED_MARKER];
+
+      Array<int> x0_bdr(pmesh.bdr_attributes.Max());
+      Array<int> x1_bdr(pmesh.bdr_attributes.Max());
+      Array<int> y0_bdr(pmesh.bdr_attributes.Max());
+      Array<int> y1_bdr(pmesh.bdr_attributes.Max());
+      Array<int> z0_bdr(pmesh.bdr_attributes.Max());
+      Array<int> z1_bdr(pmesh.bdr_attributes.Max());
+      x0_bdr = 0;
+      x1_bdr = 0;
+      y0_bdr = 0;
+      y1_bdr = 0;
+      z0_bdr = 0;
+      z1_bdr = 0;
+
+      // For cube01.mesh the boundary attributes map as:
+      // 1 -> z = 0, 2 -> y = 0, 3 -> x = 1,
+      // 4 -> y = 1, 5 -> x = 0, 6 -> z = 1.
+      z0_bdr[0] = 1;
+      y0_bdr[1] = 1;
+      x1_bdr[2] = 1;
+      y1_bdr[3] = 1;
+      x0_bdr[4] = 1;
+      z1_bdr[5] = 1;
+
+      Array<int> x0_marker, x1_marker, y0_marker, y1_marker, z0_marker, z1_marker;
+      pfes_mesh.GetEssentialVDofs(x0_bdr, x0_marker, 0);
+      pfes_mesh.GetEssentialVDofs(x1_bdr, x1_marker, 0);
+      pfes_mesh.GetEssentialVDofs(y0_bdr, y0_marker, 0);
+      pfes_mesh.GetEssentialVDofs(y1_bdr, y1_marker, 0);
+      pfes_mesh.GetEssentialVDofs(z0_bdr, z0_marker, 0);
+      pfes_mesh.GetEssentialVDofs(z1_bdr, z1_marker, 0);
+
+      for (int dof = 0; dof < pfes_mesh.GetNDofs(); dof++)
+      {
+         const int vdof = pfes_mesh.DofToVDof(dof, 0);
+         const bool on_x0 = x0_marker[vdof];
+         const bool on_x1 = x1_marker[vdof];
+         const bool on_y0 = y0_marker[vdof];
+         const bool on_y1 = y1_marker[vdof];
+         const bool on_z0 = z0_marker[vdof];
+         const bool on_z1 = z1_marker[vdof];
+
+         const int boundary_face_count =
+            int(on_x0) + int(on_x1) + int(on_y0) +
+            int(on_y1) + int(on_z0) + int(on_z1);
+         if (boundary_face_count == 3)
+         {
+            fit_marker_fixed[dof] = true;
+            for (int c = 0; c < dim; c++)
+            {
+               ess_vdofs_marker[pfes_mesh.DofToVDof(dof, c)] = 1;
+            }
+         }
+         else if (boundary_face_count == 2)
+         {
+            if (on_x0 && on_y0) { fit_marker_edge_x0y0[dof] = true; }
+            if (on_x0 && on_y1) { fit_marker_edge_x0y1[dof] = true; }
+            if (on_x1 && on_y0) { fit_marker_edge_x1y0[dof] = true; }
+            if (on_x1 && on_y1) { fit_marker_edge_x1y1[dof] = true; }
+            if (on_x0 && on_z0) { fit_marker_edge_x0z0[dof] = true; }
+            if (on_x0 && on_z1) { fit_marker_edge_x0z1[dof] = true; }
+            if (on_x1 && on_z0) { fit_marker_edge_x1z0[dof] = true; }
+            if (on_x1 && on_z1) { fit_marker_edge_x1z1[dof] = true; }
+            if (on_y0 && on_z0) { fit_marker_edge_y0z0[dof] = true; }
+            if (on_y0 && on_z1) { fit_marker_edge_y0z1[dof] = true; }
+            if (on_y1 && on_z0) { fit_marker_edge_y1z0[dof] = true; }
+            if (on_y1 && on_z1) { fit_marker_edge_y1z1[dof] = true; }
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 1)] = 1;
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 2)] = 1;
+         }
+         else if (boundary_face_count == 1)
+         {
+            if (on_x0) { fit_marker_face_x0[dof] = true; }
+            if (on_x1) { fit_marker_face_x1[dof] = true; }
+            if (on_y0) { fit_marker_face_y0[dof] = true; }
+            if (on_y1) { fit_marker_face_y1[dof] = true; }
+            if (on_z0) { fit_marker_face_z0[dof] = true; }
+            if (on_z1) { fit_marker_face_z1[dof] = true; }
+            ess_vdofs_marker[pfes_mesh.DofToVDof(dof, 2)] = 1;
+         }
+      }
+   }
+
+   for (int dof = 0; dof < pfes_mesh.GetNDofs(); dof++)
+   {
+      int nconstrained = 0;
+      for (int c = 0; c < dim; c++)
+      {
+         nconstrained += ess_vdofs_marker[pfes_mesh.DofToVDof(dof, c)];
+      }
+      fit_marker_vis_gf(dof) = nconstrained;
+   }
+   FiniteElementSpace::MarkerToList(ess_vdofs_marker, ess_vdofs);
+
+   // Visualize the selected nodes and their target positions.
+   if (glvis)
+   {
+      socketstream vis1, vis2, vis3;
+      common::VisualizeField(vis1, "localhost", 19916, fit_marker_vis_gf,
+			     "Target positions (DOFS with value 1)",
+			     0, 0, 400, 400, (dim == 2) ? "Rjm" : "");
+      common::VisualizeMesh(vis2, "localhost", 19916, pmesh, "Initial mesh",
+                            400, 0, 400, 400, (dim == 2) ? "me" : "M");
+   }
+
+   Array<AnalyticSurface *> surf_array;
+   std::vector<std::unique_ptr<AnalyticSurface>> owned_surfaces;
+
+   if (dim == 2)
+   {
+      owned_surfaces.emplace_back(new Line_Bottom(fit_markers[2]));
+      owned_surfaces.emplace_back(new Line_Left(fit_markers[3]));
+      owned_surfaces.emplace_back(new Curve_Sine_Top(fit_markers[0], a, b, c));
+      owned_surfaces.emplace_back(new Curve_Sine_Right(fit_markers[1], a, b, c));
+   }
+   else if (dim == 3)
+   {
+      // Each 2D surface and 1D edge of cube has its own AnalyticSurface.
+      owned_surfaces.emplace_back(new AxisAlignedPlane(fit_markers[0], 0, 0.0));
+      owned_surfaces.emplace_back(new CubeFace_X(fit_markers[1], a, b, c));
+      owned_surfaces.emplace_back(new AxisAlignedPlane(fit_markers[2], 1, 0.0));
+      owned_surfaces.emplace_back(new CubeFace_Y(fit_markers[3], a, b, c));
+      owned_surfaces.emplace_back(new AxisAlignedPlane(fit_markers[4], 2, 0.0));
+      owned_surfaces.emplace_back(new CubeFace_Z(fit_markers[5], a, b, c));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[6], 2, 0.0, 0.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[7], 2, 0.0, 1.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[8], 2, 1.0, 0.0));
+      owned_surfaces.emplace_back(new CubeEdge_XY(fit_markers[9], a, b, c));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[10], 1, 0.0, 0.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[11], 1, 0.0, 1.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[12], 1, 1.0, 0.0));
+      owned_surfaces.emplace_back(new CubeEdge_XZ(fit_markers[13], a, b, c));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[14], 0, 0.0, 0.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[15], 0, 0.0, 1.0));
+      owned_surfaces.emplace_back(new AxisAlignedEdge(fit_markers[16], 0, 1.0, 0.0));
+      owned_surfaces.emplace_back(new CubeEdge_YZ(fit_markers[17], a, b, c));
+   }
+   for (auto &surface : owned_surfaces)
+   {
+      surf_array.Append(surface.get());
+   }
+
+   AnalyticCompositeSurface surfaces(surf_array);
+   surfaces.ConvertPhysCoordToParam(coord_x, coord_t);
+
+   // std::ostringstream mesh_name;
+   // mesh_name << "mesh_a02_b05_c15.mesh";
+   // std::ofstream mesh_ofs(mesh_name.str().c_str());
+   // mesh_ofs.precision(8);
+   // pmesh.Print(mesh_ofs);
+   // mesh_ofs.close();
+   // return 0;
+
+   if (glvis)
+   {
+      surfaces.ConvertParamCoordToPhys(coord_t, coord_x);
+      socketstream vis1;
+      common::VisualizeMesh(vis1, "localhost", 19916, pmesh, "Mesh x->t->x",
+                            400, 0, 400, 400, (dim == 2) ? "me" : "M");
+   }
+
+   // TMOP setup.
+   TMOP_QualityMetric *metric;
+   if (dim == 2) { metric = new TMOP_Metric_002; }
+   else          { metric = new TMOP_Metric_302; }
+   metric->use_old_invariants_code = true;
+   TargetConstructor target(TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
+                            pfes_mesh.GetComm());
+   target.SetNodes(coord_x);
+   auto integ = new TMOP_Integrator(metric, &target, nullptr);
+   integ->EnableTangentialMovement(surfaces, pfes_mesh);
+
+   ParFiniteElementSpace pfes_dist(&pmesh, pfes_mesh.FEColl(), 1);
+   ParGridFunction dist(&pfes_dist);
+   dist = d;
+   ConstantCoefficient limit_coeff(1.0);
+   integ->EnableLimiting(x0, dist, limit_coeff);
+
+   // Linear solver.
+   MINRESSolver minres(pfes_mesh.GetComm());
+   minres.SetMaxIter(100);
+   minres.SetRelTol(1e-8);
+   minres.SetAbsTol(0.0);
+
+   // Nonlinear solver.
+   ParNonlinearForm nlf(&pfes_mesh);
+   nlf.SetEssentialVDofs(ess_vdofs);
+   nlf.AddDomainIntegrator(integ);
+   const IntegrationRule &ir =
+    IntRules.Get(pfes_mesh.GetFE(0)->GetGeomType(), quad_order);
+   TMOPNewtonSolver solver(pfes_mesh.GetComm(), ir, 0);
+   solver.SetOperator(nlf);
+   solver.SetPreconditioner(minres);
+   solver.SetPrintLevel(1);
+   solver.SetMaxIter(50);
+   solver.SetRelTol(1e-6);
+   solver.SetAbsTol(0.0);
+
+   // Solve.
+   Vector zero(0);
+   coord_t.SetTrueVector();
+   solver.Mult(zero, coord_t.GetTrueVector());
+   coord_t.SetFromTrueVector();
+   surfaces.ConvertParamCoordToPhys(coord_t, coord_x);
+   if(glvis)
+   {
+      socketstream vis2;
+      common::VisualizeMesh(vis2, "localhost", 19916, pmesh, "Final mesh",
+                            800, 0, 400, 400, (dim == 2) ? "me" : "M");
+   }
+
+   delete metric;
+   return 0;
+}

miniapps/meshing/tmop-tangent.hpp

@@ -0,0 +1,772 @@
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// MFEM Mesh Optimizer Miniapp - Serial/Parallel Shared Code
+
+#include "mfem.hpp"
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+using namespace std;
+
+// x = t.
+// y = 0.
+// The distance is the error in y.
+class Line_Bottom : public Analytic2DCurve
+{
+   public:
+   Line_Bottom(const Array<bool> &marker) : Analytic2DCurve(marker) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = x;
+      dist = y - 0.0;
+   }
+   void xy_of_t(double t, double dist, double &x, double &y) const override
+   {
+      x = t;
+      y = 0.0 + dist;
+   }
+
+   virtual double dx_dt(double t) const override { return 1.0; }
+   virtual double dy_dt(double t) const override { return 0.0; }
+   virtual double dx_dtdt(double t) const override { return 0.0; }
+   virtual double dy_dtdt(double t) const override { return 0.0; }
+};
+
+// x = 0.
+// y = t.
+// The distance is the error in x.
+class Line_Left : public Analytic2DCurve
+{
+   public:
+   Line_Left(const Array<bool> &marker) : Analytic2DCurve(marker) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = y;
+      dist = x - 0.0;
+   }
+   void xy_of_t(double t, double dist, double &x, double &y) const override
+   {
+      x = 0.0 + dist;
+      y = t;
+   }
+
+   virtual double dx_dt(double t) const override { return 0.0; }
+   virtual double dy_dt(double t) const override { return 1.0; }
+   virtual double dx_dtdt(double t) const override { return 0.0; }
+   virtual double dy_dtdt(double t) const override { return 0.0; }
+};
+
+// x = 1.5 t.
+// y = 1 + 0.5 t.
+// The distance is the error in y.
+class Line_Top : public Analytic2DCurve
+{
+   public:
+   Line_Top(const Array<bool> &marker) : Analytic2DCurve(marker) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = x / 1.5;
+      dist = y - (1 + 0.5 * t);
+   }
+   void xy_of_t(double t, double dist, double &x, double &y) const override
+   {
+      x = 1.5 * t;
+      y = dist + 1.0 + 0.5 * t;
+   }
+
+   virtual double dx_dt(double t) const override { return 1.5; }
+   virtual double dy_dt(double t) const override { return 0.5; }
+   virtual double dx_dtdt(double t) const override { return 0.0; }
+   virtual double dy_dtdt(double t) const override { return 0.0; }
+};
+
+// x = 1 + 0.5 t.
+// y = 1.5 t.
+// The distance is the error in x.
+class Line_Right : public Analytic2DCurve
+{
+   public:
+   Line_Right(const Array<bool> &marker) : Analytic2DCurve(marker) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = y / 1.5;
+      dist = x - (1 + 0.5 * t);
+   }
+   void xy_of_t(double t, double dist, double &x, double &y) const override
+   {
+      x = dist + 1.0 + 0.5 * t;
+      y = 1.5 * t;
+   }
+
+   virtual double dx_dt(double t) const override { return 0.5; }
+   virtual double dy_dt(double t) const override { return 1.5; }
+   virtual double dx_dtdt(double t) const override { return 0.0; }
+   virtual double dy_dtdt(double t) const override { return 0.0; }
+};
+
+// x = [1.0 + a sin(c pi) + b] t.
+// y = 1 + a sin(c pi t) + b t.
+// The distance is the error in y.
+class Curve_Sine_Top : public Analytic2DCurve
+{
+private:
+   const double a, b, c, x_scale;
+
+public:
+   Curve_Sine_Top(const Array<bool> &marker, double a_, double b_, double c_)
+    : Analytic2DCurve(marker),
+    a(a_), b(b_), c(c_), x_scale(1.0 + a * sin(c * M_PI) + b) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = x / x_scale;
+      dist = y - (1.0 + a * sin(c * M_PI * t) + b * t);
+   }
+   void xy_of_t(double t, double dist, double &x, double &y) const override
+   {
+      x = x_scale * t;
+      y = dist + 1.0 + a * sin(c * M_PI * t) + b * t;
+   }
+
+   virtual double dx_dt(double t) const override
+   { return x_scale; }
+   virtual double dy_dt(double t) const override
+   { return a * c * M_PI * cos(c * M_PI * t) + b; }
+
+   virtual double dx_dtdt(double t) const override { return 0.0; }
+   virtual double dy_dtdt(double t) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+};
+
+// x = 1 + a sin(s pi t) + b t.
+// y = [1.0 + a sin(s pi) + b] t.
+// The distance is the error in x.
+class Curve_Sine_Right : public Analytic2DCurve
+{
+private:
+   const double a, b, c, y_scale;
+
+public:
+   Curve_Sine_Right(const Array<bool> &marker, double a_, double b_, double c_)
+    : Analytic2DCurve(marker),
+      a(a_), b(b_), c(c_), y_scale(1.0 + a * sin(c * M_PI) + b) { }
+
+   void t_of_xy(double x, double y, double &dist, double &t) const override
+   {
+      t    = y / y_scale;
+      dist = x - (1.0 + a * sin(c * M_PI * t) + b * t);
+   }
+   void xy_of_t(double t,  double dist, double &x, double &y) const override
+   {
+      x = dist + 1.0 + a * sin(c * M_PI * t) + b * t;
+      y = y_scale * t;
+   }
+
+   virtual double dx_dt(double t) const override
+   { return a * c * M_PI * cos(c * M_PI * t) + b; }
+   virtual double dy_dt(double t) const override
+   { return y_scale; }
+
+
+   virtual double dx_dtdt(double t) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+   virtual double dy_dtdt(double t) const override
+   { return 0.0; }
+};
+
+class AxisAlignedEdge : public Analytic3DCurve
+{
+private:
+   const int axis;
+   const double fixed_1;
+   const double fixed_2;
+
+public:
+   // Creates straight edge aligned with an axis (0 <-> x, 1 <-> y, 2 <-> z).
+   // fixed_1 and fixed_2 give the coordinates on the normal plane to the axis.
+   AxisAlignedEdge(const Array<bool> &marker, int axis_,
+                   double fixed_1_, double fixed_2_)
+      : Analytic3DCurve(marker),
+        axis(axis_),
+        fixed_1(fixed_1_),
+        fixed_2(fixed_2_) { }
+
+   void t_of_xyz(double x, double y, double z,
+                 double &dist1, double &dist2, double &t) const override
+   {
+      if (axis == 0)
+      {
+         t = x;
+         dist1 = y - fixed_1;
+         dist2 = z - fixed_2;
+      }
+      else if (axis == 1)
+      {
+         t = y;
+         dist1 = x - fixed_1;
+         dist2 = z - fixed_2;
+      }
+      else
+      {
+         t = z;
+         dist1 = x - fixed_1;
+         dist2 = y - fixed_2;
+      }
+   }
+
+   void xyz_of_t(double t, double dist1, double dist2,
+                 double &x, double &y, double &z) const override
+   {
+      if (axis == 0)
+      {
+         x = t;
+         y = fixed_1 + dist1;
+         z = fixed_2 + dist2;
+      }
+      else if (axis == 1)
+      {
+         x = fixed_1 + dist1;
+         y = t;
+         z = fixed_2 + dist2;
+      }
+      else
+      {
+         x = fixed_1 + dist1;
+         y = fixed_2 + dist2;
+         z = t;
+      }
+   }
+
+   double dx_dt(double) const override { return axis == 0 ? 1.0 : 0.0; }
+   double dy_dt(double) const override { return axis == 1 ? 1.0 : 0.0; }
+   double dz_dt(double) const override { return axis == 2 ? 1.0 : 0.0; }
+
+   double dx_dtdt(double) const override { return 0.0; }
+   double dy_dtdt(double) const override { return 0.0; }
+   double dz_dtdt(double) const override { return 0.0; }
+};
+
+class AxisAlignedPlane : public Analytic3DSurface
+{
+private:
+   const int normal_axis;
+   const double fixed_value;
+
+public:
+   AxisAlignedPlane(const Array<bool> &marker, int normal_axis_,
+                    double fixed_value_)
+      : Analytic3DSurface(marker),
+        normal_axis(normal_axis_),
+        fixed_value(fixed_value_) { }
+
+   void uv_of_xyz(double x, double y, double z,
+                  double &dist, double &u, double &v) const override
+   {
+      if (normal_axis == 0)
+      {
+         dist = x - fixed_value;
+         u = y;
+         v = z;
+      }
+      else if (normal_axis == 1)
+      {
+         dist = y - fixed_value;
+         u = x;
+         v = z;
+      }
+      else
+      {
+         dist = z - fixed_value;
+         u = x;
+         v = y;
+      }
+   }
+
+   void xyz_of_uv(double u, double v, double dist,
+                  double &x, double &y, double &z) const override
+   {
+      if (normal_axis == 0)
+      {
+         x = fixed_value + dist;
+         y = u;
+         z = v;
+      }
+      else if (normal_axis == 1)
+      {
+         x = u;
+         y = fixed_value + dist;
+         z = v;
+      }
+      else
+      {
+         x = u;
+         y = v;
+         z = fixed_value + dist;
+      }
+   }
+
+   double dx_du(double, double) const override { return normal_axis == 1 || normal_axis == 2 ? 1.0 : 0.0; }
+   double dy_du(double, double) const override { return normal_axis == 0 ? 1.0 : 0.0; }
+   double dz_du(double, double) const override { return 0.0; }
+   double dx_dv(double, double) const override { return 0.0; }
+   double dy_dv(double, double) const override { return normal_axis == 2 ? 1.0 : 0.0; }
+   double dz_dv(double, double) const override { return normal_axis == 0 || normal_axis == 1 ? 1.0 : 0.0; }
+
+   double dx_dudu(double, double) const override { return 0.0; }
+   double dy_dudu(double, double) const override { return 0.0; }
+   double dz_dudu(double, double) const override { return 0.0; }
+   double dx_dudv(double, double) const override { return 0.0; }
+   double dy_dudv(double, double) const override { return 0.0; }
+   double dz_dudv(double, double) const override { return 0.0; }
+   double dx_dvdv(double, double) const override { return 0.0; }
+   double dy_dvdv(double, double) const override { return 0.0; }
+   double dz_dvdv(double, double) const override { return 0.0; }
+};
+
+class CubeEdge_XY : public Analytic3DCurve
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeEdge_XY(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DCurve(marker), a(a_), b(b_), c(c_) { }
+
+   void t_of_xyz(double x, double y, double z,
+                 double &dist1, double &dist2, double &t) const override
+   {
+      t = std::max(0.0, std::min(1.0, z));
+      for (int it = 0; it < 12; it++)
+      {
+         const double f = z_of_t(t) - z;
+         const double df = dz_dt(t);
+         if (std::abs(df) < 1e-12) { break; }
+         const double dt = -f / df;
+         t = std::max(0.0, std::min(1.0, t + dt));
+         if (std::abs(dt) < 1e-12) { break; }
+      }
+      dist1 = x - x_of_t(t);
+      dist2 = y - y_of_t(t);
+   }
+
+   void xyz_of_t(double t, double dist1, double dist2,
+                 double &x, double &y, double &z) const override
+   {
+      x = x_of_t(t) + dist1;
+      y = y_of_t(t) + dist2;
+      z = z_of_t(t);
+   }
+
+   double dx_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+   double dy_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+   double dz_dt(double t) const override
+   { return 1.0 + 0.5 * a * s * s * M_PI * cos(0.5 * M_PI * t) + b; }
+
+   double dx_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+   double dy_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+   double dz_dtdt(double t) const override
+   { return -0.25 * a * s * s * M_PI * M_PI * sin(0.5 * M_PI * t); }
+
+private:
+   double x_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+   double y_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+   double z_of_t(double t) const
+   { return t + a * s * s * sin(0.5 * M_PI * t) + b * t; }
+};
+
+class CubeEdge_XZ : public Analytic3DCurve
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeEdge_XZ(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DCurve(marker), a(a_), b(b_), c(c_) { }
+
+   void t_of_xyz(double x, double y, double z,
+                 double &dist1, double &dist2, double &t) const override
+   {
+      // Becuase the transformation is nonlinear, the problem to find t
+      // to parameterize the curve is solved using Newton's method.
+      t = std::max(0.0, std::min(1.0, y));
+      for (int it = 0; it < 12; it++)
+      {
+         const double f = y_of_t(t) - y;
+         const double df = dy_dt(t);
+         if (std::abs(df) < 1e-12) { break; }
+         const double dt = -f / df;
+         t = std::max(0.0, std::min(1.0, t + dt));
+         if (std::abs(dt) < 1e-12) { break; }
+      }
+      dist1 = x - x_of_t(t);
+      dist2 = z - z_of_t(t);
+   }
+
+   void xyz_of_t(double t, double dist1, double dist2,
+                 double &x, double &y, double &z) const override
+   {
+      x = x_of_t(t) + dist1;
+      y = y_of_t(t);
+      z = z_of_t(t) + dist2;
+   }
+
+   double dx_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+   double dy_dt(double t) const override
+   { return 1.0 + 0.5 * a * s * s * M_PI * cos(0.5 * M_PI * t) + b; }
+   double dz_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+
+   double dx_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+   double dy_dtdt(double t) const override
+   { return -0.25 * a * s * s * M_PI * M_PI * sin(0.5 * M_PI * t); }
+   double dz_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+
+private:
+   double x_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+   double y_of_t(double t) const
+   { return t + a * s * s * sin(0.5 * M_PI * t) + b * t; }
+   double z_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+};
+
+class CubeEdge_YZ : public Analytic3DCurve
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeEdge_YZ(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DCurve(marker), a(a_), b(b_), c(c_) { }
+
+   void t_of_xyz(double x, double y, double z,
+                 double &dist1, double &dist2, double &t) const override
+   {
+      // Becuase the transformation is nonlinear, the problem to find t
+      // to parameterize the curve is solved using Newton's method.
+      t = std::max(0.0, std::min(1.0, x));
+      for (int it = 0; it < 12; it++)
+      {
+         const double f = x_of_t(t) - x;
+         const double df = dx_dt(t);
+         if (std::abs(df) < 1e-12) { break; }
+         const double dt = -f / df;
+         t = std::max(0.0, std::min(1.0, t + dt));
+         if (std::abs(dt) < 1e-12) { break; }
+      }
+      dist1 = y - y_of_t(t);
+      dist2 = z - z_of_t(t);
+   }
+
+   void xyz_of_t(double t, double dist1, double dist2,
+                 double &x, double &y, double &z) const override
+   {
+      x = x_of_t(t);
+      y = y_of_t(t) + dist1;
+      z = z_of_t(t) + dist2;
+   }
+
+   double dx_dt(double t) const override
+   { return 1.0 + 0.5 * a * s * s * M_PI * cos(0.5 * M_PI * t) + b; }
+   double dy_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+   double dz_dt(double t) const override
+   { return a * s * c * M_PI * cos(c * M_PI * t) + b; }
+
+   double dx_dtdt(double t) const override
+   { return -0.25 * a * s * s * M_PI * M_PI * sin(0.5 * M_PI * t); }
+   double dy_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+   double dz_dtdt(double t) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * t); }
+
+private:
+   double x_of_t(double t) const
+   { return t + a * s * s * sin(0.5 * M_PI * t) + b * t; }
+   double y_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+   double z_of_t(double t) const
+   { return 1.0 + a * s * sin(c * M_PI * t) + b * t; }
+};
+
+class CubeFace_X : public Analytic3DSurface
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeFace_X(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DSurface(marker), a(a_), b(b_), c(c_) { }
+
+   void uv_of_xyz(double x, double y, double z,
+                  double &dist, double &u, double &v) const override
+   {
+      // Becuase the transformation is nonlinear, the problem to find u and v 
+      // to parameterize the surface is solved using Newton's method.
+      u = std::max(0.0, std::min(1.0, y));
+      v = std::max(0.0, std::min(1.0, z));
+      for (int it = 0; it < 15; it++)
+      {
+         const double r1 = y_of_uv(u, v) - y;
+         const double r2 = z_of_uv(u, v) - z;
+         const double j11 = dy_du(u, v);
+         const double j12 = dy_dv(u, v);
+         const double j21 = dz_du(u, v);
+         const double j22 = dz_dv(u, v);
+         const double det = j11 * j22 - j12 * j21;
+         if (std::abs(det) < 1e-12) { break; }
+         const double du = (-r1 * j22 + r2 * j12) / det;
+         const double dv = (-j11 * r2 + j21 * r1) / det;
+         u = std::max(0.0, std::min(1.0, u + du));
+         v = std::max(0.0, std::min(1.0, v + dv));
+         if (std::max(std::abs(du), std::abs(dv)) < 1e-12) { break; }
+      }
+      dist = x - x_of_uv(u, v);
+   }
+
+   void xyz_of_uv(double u, double v, double dist,
+                  double &x, double &y, double &z) const override
+   {
+      x = x_of_uv(u, v) + dist;
+      y = y_of_uv(u, v);
+      z = z_of_uv(u, v);
+   }
+
+   double dx_du(double u, double v) const override
+   { return a * c * M_PI * cos(c * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dy_du(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * cos(0.5 * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dz_du(double u, double v) const override
+   { return a * s * c * M_PI * cos(c * M_PI * u) * sin(0.5 * M_PI * v) + b * v; }
+   double dx_dv(double u, double v) const override
+   { return a * c * M_PI * sin(c * M_PI * u) * cos(c * M_PI * v) + b * u; }
+   double dy_dv(double u, double v) const override
+   { return a * s * c * M_PI * sin(0.5 * M_PI * u) * cos(c * M_PI * v) + b * u; }
+   double dz_dv(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * sin(c * M_PI * u) * cos(0.5 * M_PI * v) + b * u; }
+
+   double dx_dudu(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dudu(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dz_dudu(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+   double dx_dudv(double u, double v) const override
+   { return a * c * c * M_PI * M_PI * cos(c * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dy_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(0.5 * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dz_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(c * M_PI * u) * cos(0.5 * M_PI * v) + b; }
+   double dx_dvdv(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dvdv(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dz_dvdv(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+
+private:
+   double x_of_uv(double u, double v) const
+   { return 1.0 + a * sin(c * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+   double y_of_uv(double u, double v) const
+   { return u + a * s * sin(0.5 * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+   double z_of_uv(double u, double v) const
+   { return v + a * s * sin(c * M_PI * u) * sin(0.5 * M_PI * v) + b * u * v; }
+};
+
+class CubeFace_Y : public Analytic3DSurface
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeFace_Y(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DSurface(marker), a(a_), b(b_), c(c_) { }
+
+   void uv_of_xyz(double x, double y, double z,
+                  double &dist, double &u, double &v) const override
+   {
+      // Becuase the transformation is nonlinear, the problem to find u and v 
+      // to parameterize the surface is solved using Newton's method.
+      u = std::max(0.0, std::min(1.0, x));
+      v = std::max(0.0, std::min(1.0, z));
+      for (int it = 0; it < 15; it++)
+      {
+         const double r1 = x_of_uv(u, v) - x;
+         const double r2 = z_of_uv(u, v) - z;
+         const double j11 = dx_du(u, v);
+         const double j12 = dx_dv(u, v);
+         const double j21 = dz_du(u, v);
+         const double j22 = dz_dv(u, v);
+         const double det = j11 * j22 - j12 * j21;
+         if (std::abs(det) < 1e-12) { break; }
+         const double du = (-r1 * j22 + r2 * j12) / det;
+         const double dv = (-j11 * r2 + j21 * r1) / det;
+         u = std::max(0.0, std::min(1.0, u + du));
+         v = std::max(0.0, std::min(1.0, v + dv));
+         if (std::max(std::abs(du), std::abs(dv)) < 1e-12) { break; }
+      }
+      dist = y - y_of_uv(u, v);
+   }
+
+   void xyz_of_uv(double u, double v, double dist,
+                  double &x, double &y, double &z) const override
+   {
+      x = x_of_uv(u, v);
+      y = y_of_uv(u, v) + dist;
+      z = z_of_uv(u, v);
+   }
+
+   double dx_du(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * cos(0.5 * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dy_du(double u, double v) const override
+   { return a * c * M_PI * cos(c * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dz_du(double u, double v) const override
+   { return a * s * c * M_PI * cos(c * M_PI * u) * sin(0.5 * M_PI * v) + b * v; }
+   double dx_dv(double u, double v) const override
+   { return a * s * c * M_PI * sin(0.5 * M_PI * u) * cos(c * M_PI * v) + b * u; }
+   double dy_dv(double u, double v) const override
+   { return a * c * M_PI * sin(c * M_PI * u) * cos(c * M_PI * v) + b * u; }
+   double dz_dv(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * sin(c * M_PI * u) * cos(0.5 * M_PI * v) + b * u; }
+
+   double dx_dudu(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dudu(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+   double dz_dudu(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+   double dx_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(0.5 * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dy_dudv(double u, double v) const override
+   { return a * c * c * M_PI * M_PI * cos(c * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dz_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(c * M_PI * u) * cos(0.5 * M_PI * v) + b; }
+   double dx_dvdv(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dvdv(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+   double dz_dvdv(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+
+private:
+   double x_of_uv(double u, double v) const
+   { return u + a * s * sin(0.5 * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+   double y_of_uv(double u, double v) const
+   { return 1.0 + a * sin(c * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+   double z_of_uv(double u, double v) const
+   { return v + a * s * sin(c * M_PI * u) * sin(0.5 * M_PI * v) + b * u * v; }
+};
+
+class CubeFace_Z : public Analytic3DSurface
+{
+private:
+   const double a, b, c;
+   const double s = sin(c * M_PI);
+
+public:
+   CubeFace_Z(const Array<bool> &marker, double a_, double b_, double c_)
+      : Analytic3DSurface(marker), a(a_), b(b_), c(c_) { }
+
+   void uv_of_xyz(double x, double y, double z,
+                  double &dist, double &u, double &v) const override
+   {
+      // Becuase the transformation is nonlinear, the problem to find u and v 
+      // to parameterize the surface is solved using Newton's method.
+      u = std::max(0.0, std::min(1.0, x));
+      v = std::max(0.0, std::min(1.0, y));
+      for (int it = 0; it < 15; it++)
+      {
+         const double r1 = x_of_uv(u, v) - x;
+         const double r2 = y_of_uv(u, v) - y;
+         const double j11 = dx_du(u, v);
+         const double j12 = dx_dv(u, v);
+         const double j21 = dy_du(u, v);
+         const double j22 = dy_dv(u, v);
+         const double det = j11 * j22 - j12 * j21;
+         if (std::abs(det) < 1e-12) { break; }
+         const double du = (-r1 * j22 + r2 * j12) / det;
+         const double dv = (-j11 * r2 + j21 * r1) / det;
+         u = std::max(0.0, std::min(1.0, u + du));
+         v = std::max(0.0, std::min(1.0, v + dv));
+         if (std::max(std::abs(du), std::abs(dv)) < 1e-12) { break; }
+      }
+      dist = z - z_of_uv(u, v);
+   }
+
+   void xyz_of_uv(double u, double v, double dist,
+                  double &x, double &y, double &z) const override
+   {
+      x = x_of_uv(u, v);
+      y = y_of_uv(u, v);
+      z = z_of_uv(u, v) + dist;
+   }
+
+   double dx_du(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * cos(0.5 * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dy_du(double u, double v) const override
+   { return a * s * c * M_PI * cos(c * M_PI * u) * sin(0.5 * M_PI * v) + b * v; }
+   double dz_du(double u, double v) const override
+   { return a * c * M_PI * cos(c * M_PI * u) * sin(c * M_PI * v) + b * v; }
+   double dx_dv(double u, double v) const override
+   { return a * s * c * M_PI * sin(0.5 * M_PI * u) * cos(c * M_PI * v) + b * u; }
+   double dy_dv(double u, double v) const override
+   { return 1.0 + 0.5 * a * s * M_PI * sin(c * M_PI * u) * cos(0.5 * M_PI * v) + b * u; }
+   double dz_dv(double u, double v) const override
+   { return a * c * M_PI * sin(c * M_PI * u) * cos(c * M_PI * v) + b * u; }
+
+   double dx_dudu(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dudu(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+   double dz_dudu(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+   double dx_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(0.5 * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dy_dudv(double u, double v) const override
+   { return 0.5 * a * s * c * M_PI * M_PI * cos(c * M_PI * u) * cos(0.5 * M_PI * v) + b; }
+   double dz_dudv(double u, double v) const override
+   { return a * c * c * M_PI * M_PI * cos(c * M_PI * u) * cos(c * M_PI * v) + b; }
+   double dx_dvdv(double u, double v) const override
+   { return -a * s * c * c * M_PI * M_PI * sin(0.5 * M_PI * u) * sin(c * M_PI * v); }
+   double dy_dvdv(double u, double v) const override
+   { return -0.25 * a * s * M_PI * M_PI * sin(c * M_PI * u) * sin(0.5 * M_PI * v); }
+   double dz_dvdv(double u, double v) const override
+   { return -a * c * c * M_PI * M_PI * sin(c * M_PI * u) * sin(c * M_PI * v); }
+
+private:
+   double x_of_uv(double u, double v) const
+   { return u + a * s * sin(0.5 * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+   double y_of_uv(double u, double v) const
+   { return v + a * s * sin(c * M_PI * u) * sin(0.5 * M_PI * v) + b * u * v; }
+   double z_of_uv(double u, double v) const
+   { return 1.0 + a * sin(c * M_PI * u) * sin(c * M_PI * v) + b * u * v; }
+};

miniapps/tools/tmop-check-metric.cpp

@@ -32,6 +32,7 @@ int main(int argc, char *argv[])
    bool a_metric_version = false;
    int convergence_iter = 10;
    bool verbose = false;
+   bool old_inv = false;
 
    // Choose metric.
    OptionsParser args(argc, argv);
@@ -40,6 +41,8 @@ int main(int argc, char *argv[])
                   "Use the A-version of the metric, if available.");
    args.AddOption(&verbose, "-v", "-verbose", "-no-v", "--no-verbose",
                   "Enable extra screen output.");
+   args.AddOption(&old_inv, "-oi", "--old-invar", "-ni", "--new-invar",
+                  "Use old or new code or invariants.");
    args.AddOption(&convergence_iter, "-i", "--iterations",
                   "Number of iterations to check convergence of derivatives.");
 
@@ -109,6 +112,8 @@ int main(int argc, char *argv[])
       default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
    }
 
+   if (old_inv) { metric->use_old_invariants_code = true; }
+
    const int dim = (metric_id < 300) ? 2 : 3;
    Mesh *mesh;
    if (dim == 2)

tests/unit/fem/test_pa_kernels.cpp

@@ -228,6 +228,25 @@ void gradf1(const Vector &x, Vector &u)
    if (x.Size() >= 3) { u(2) = 4*pow(x(2), 3); }
 }
 
+void f_matrix(const Vector &x, DenseMatrix &M)
+{
+   M(0,0) = x(0);
+   if (x.Size() >= 2)
+   {
+      M(1,1) = x(1);
+      M(0,1) = x(0) * x(1);
+      M(1,0) = x(0) * x(1) * 2;
+   }
+   if (x.Size() >= 3)
+   {
+      M(2,2) = x(2);
+      M(0,2) = x(0) * x(1) * x(2);
+      M(2,0) = x(0) * x(1) * x(2) * 3;
+      M(1,2) = x(1) * x(2);
+      M(2,1) = x(1) * x(2) * 2;
+   }
+}
+
 void pa_gradient_testnd(int dim, FECType fec_type,
                         real_t (*f1)(const Vector &),
                         void (*gradf1)(const Vector &, Vector &))
@@ -784,6 +803,43 @@ TEST_CASE("PA Boundary Mass", "[PartialAssembly], [GPU]")
    REQUIRE(y_fa.Normlinf() == MFEM_Approx(0.0));
 }
 
+TEST_CASE("PA BoundaryVectorMassIntegrator",
+          "[PartialAssembly] [BoundaryVectorMassIntegrator] [CUDA]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   auto fname = GENERATE("../../data/star.mesh", "../../data/star-q3.mesh",
+                         "../../data/fichera.mesh", "../../data/fichera-q3.mesh");
+   auto order = !all_tests ? 2 : GENERATE(1, 2, 3);
+
+   Mesh mesh(fname);
+   int dim = mesh.Dimension();
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec, dim);
+
+   GridFunction x(&fes), y_fa(&fes), y_pa(&fes);
+   x.Randomize(1);
+
+   MatrixFunctionCoefficient mat_coeff(dim, f_matrix);
+
+   BilinearForm blf_fa(&fes);
+   blf_fa.AddBdrFaceIntegrator(new BoundaryVectorMassIntegrator(mat_coeff));
+   blf_fa.Assemble();
+   blf_fa.Finalize();
+   blf_fa.Mult(x, y_fa);
+
+   BilinearForm blf_pa(&fes);
+   blf_pa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   blf_pa.AddBdrFaceIntegrator(new BoundaryVectorMassIntegrator(mat_coeff));
+
+   blf_pa.Assemble();
+   blf_pa.Mult(x, y_pa);
+
+   y_fa -= y_pa;
+
+   REQUIRE(y_fa.Normlinf() == MFEM_Approx(0.0));
+}
+
 namespace
 {
 template <typename T> struct ParTypeHelper { };