mfem · kmittal2 · Apr 27, 2025 · Apr 29, 2025 · Apr 29, 2025 · Apr 29, 2025
diff --git a/fem/tmop.cpp b/fem/tmop.cpp
@@ -140,6 +140,28 @@ void add_3D(const scalartype &scalar, const std::vector<type> &u,
 
 /* Metric definitions */
 
+// W = ||T||^2 - 2*det(T).
+template <typename type>
+type mu4_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   auto fnorm2 = fnorm2_2D(T);
+   auto det = det_2D(T);
+   return fnorm2 - 2*det;
+};
+
+// W = ||T-I||^2.
+template <typename type>
+type mu14_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   DenseMatrix Id(2,2); Id = 0.0;
+   Id(0,0) = 1; Id(1,1) = 1;
+
+   std::vector<type> Mat;
+   add_2D(real_t{-1.0}, T, &Id, Mat);
+
+   return fnorm2_2D(Mat);
+};
+
 // W = |T-T'|^2, where T'= |T|*I/sqrt(2).
 template <typename type>
 type mu85_ad(const std::vector<type> &T, const std::vector<type> &W)
@@ -163,6 +185,63 @@ type mu98_ad(const std::vector<type> &T, const std::vector<type> &W)
    return fnorm2_2D(Mat)/det_2D(T);
 };
 
+template <typename type>
+type make_one_type()
+{
+   return 1.0;
+}
+// add specialization for AD1Type
+template <>
+AD1Type make_one_type<AD1Type>()
+{
+   return AD1Type{1.0, 0.0};
+}
+// add specialization for AD2Type
+template <>
+AD2Type make_one_type<AD2Type>()
+{
+   return AD2Type{AD1Type{1.0, 0.0}, AD1Type{0.0, 0.0}};
+}
+
+using TWCUO = TMOP_WorstCaseUntangleOptimizer_Metric;
+template <typename type>
+type wcuo_ad(type mu,
+             const std::vector<type> &T, const std::vector<type> &W,
+             real_t alpha, real_t min_detT, real_t detT_ep,
+             int exponent, real_t max_muT, real_t muT_ep,
+             TWCUO::BarrierType bt,
+             TWCUO::WorstCaseType wct)
+{
+   type one = make_one_type<type>();
+   type zero = 0.0*one;
+   type denom = one;
+   if (bt == TWCUO::BarrierType::Shifted)
+   {
+      auto val1 = alpha*min_detT-detT_ep < 0.0 ?
+                  (alpha*min_detT-detT_ep)*one :
+                  zero;
+      denom = 2.0*(det_2D(T)-val1);
+   }
+   else if (bt == TWCUO::BarrierType::Pseudo)
+   {
+      auto detT = det_2D(T);
+      denom = detT + sqrt(detT*detT + detT_ep*detT_ep);
+   }
+   mu = mu/denom;
+
+   if (wct == TWCUO::WorstCaseType::PMean)
+   {
+      auto exp = exponent*one;
+      mu = pow(mu, exp);
+   }
+   else if (wct == TWCUO::WorstCaseType::Beta)
+   {
+      auto beta = (max_muT+muT_ep)*one;
+      mu = mu/(beta-mu);
+   }
+   return mu;
+}
+
 // W = 1/(tau^0.5) |T-I|^2.
 template <typename type>
 type mu342_ad(const std::vector<type> &T, const std::vector<type> &W)
@@ -421,7 +500,7 @@ void TMOP_QualityMetric::DefaultAssembleH(const DenseTensor &H,
          {
             for (int cc = 0; cc < dim; cc++)
             {
-               const double entry_rr_cc = Hrc(rr, cc);
+               const real_t entry_rr_cc = Hrc(rr, cc);
 
                for (int i = 0; i < dof; i++)
                {
@@ -481,6 +560,30 @@ void TMOP_Combo_QualityMetric::EvalPW(const DenseMatrix &Jpt,
    }
 }
 
+AD1Type TMOP_Combo_QualityMetric::EvalW_AD1(const std::vector<AD1Type> &T,
+                                            const std::vector<AD1Type> &W)
+const
+{
+   AD1Type metric = {0., 0.};
+   for (int i = 0; i < tmop_q_arr.Size(); i++)
+   {
+      metric += wt_arr[i]*tmop_q_arr[i]->EvalW_AD1(T, W);
+   }
+   return metric;
+}
+
+AD2Type TMOP_Combo_QualityMetric::EvalW_AD2(const std::vector<AD2Type> &T,
+                                            const std::vector<AD2Type> &W)
+const
+{
+   AD2Type metric = {{0., 0.},{0., 0.}};
+   for (int i = 0; i < tmop_q_arr.Size(); i++)
+   {
+      metric += wt_arr[i]*tmop_q_arr[i]->EvalW_AD2(T, W);
+   }
+   return metric;
+}
+
 void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
                                          const DenseMatrix &DS,
                                          const real_t weight,
@@ -645,6 +748,64 @@ real_t TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier(
    return tmop_metric.EvalW(Jpt)/denominator;
 }
 
+AD1Type TMOP_WorstCaseUntangleOptimizer_Metric::EvalW_AD1(
+   const std::vector<AD1Type> &T,
+   const std::vector<AD1Type> &W) const
+{
+   return wcuo_ad(tmop_metric.EvalW_AD1(T,W), T, W, alpha, min_detT, detT_ep,
+                  exponent, max_muT, muT_ep, btype, wctype);
+}
+
+AD2Type TMOP_WorstCaseUntangleOptimizer_Metric::EvalW_AD2(
+   const std::vector<AD2Type> &T,
+   const std::vector<AD2Type> &W) const
+{
+   return wcuo_ad(tmop_metric.EvalW_AD2(T,W), T, W, alpha, min_detT, detT_ep,
+                  exponent, max_muT, muT_ep, btype, wctype);
+}
+
+void TMOP_WorstCaseUntangleOptimizer_Metric::EvalP(const DenseMatrix &Jpt,
+                                                   DenseMatrix &P) const
+{
+   auto mu_ad_fn = [this](
+                      std::vector<AD1Type> &T, std::vector<AD1Type> &W)
+   {
+      return EvalW_AD1(T,W);
+   };
+   if (tmop_metric.Id() == 4 || tmop_metric.Id() == 14 || tmop_metric.Id() == 66)
+   {
+      ADGrad(mu_ad_fn, P, Jpt);
+      return;
+   }
+   MFEM_ABORT("EvalW_AD1 not implemented with this metric for "
+              "TMOP_WorstCaseUntangleOptimizer_Metric. "
+              "Please use metric 4/14/66.");
+}
+
+void TMOP_WorstCaseUntangleOptimizer_Metric::AssembleH(
+   const DenseMatrix &Jpt,
+   const DenseMatrix &DS,
+   const real_t weight,
+   DenseMatrix &A) const
+{
+   DenseTensor H(Jpt.Height(), Jpt.Height(), Jpt.TotalSize());
+   H = 0.0;
+   auto mu_ad_fn = [this](
+                      std::vector<AD2Type> &T, std::vector<AD2Type> &W)
+   {
+      return EvalW_AD2(T,W);
+   };
+   if (tmop_metric.Id() == 4 || tmop_metric.Id() == 14 ||tmop_metric.Id() == 66)
+   {
+      ADHessian(mu_ad_fn, H, Jpt);
+      this->DefaultAssembleH(H,DS,weight,A);
+      return;
+   }
+   MFEM_ABORT("EvalW_AD1 not implemented with this metric for "
+              "TMOP_WorstCaseUntangleOptimizer_Metric. "
+              "Please use metric 4/14/66.");
+}
+
 real_t TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
 {
    ie.SetJacobian(Jpt.GetData());
@@ -850,6 +1011,25 @@ void TMOP_Metric_004::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI2b(-2.0*weight, A.GetData());
 }
 
+template <typename type>
+type TMOP_Metric_004::EvalW_AD_impl(const std::vector<type> &T,
+                                    const std::vector<type> &W) const
+{
+   return mu4_ad(T, W);
+}
+
+AD1Type TMOP_Metric_004::EvalW_AD1(const std::vector<AD1Type> &T,
+                                   const std::vector<AD1Type> &W) const
+{
+   return EvalW_AD_impl<AD1Type>(T,W);
+}
+
+AD2Type TMOP_Metric_004::EvalW_AD2(const std::vector<AD2Type> &T,
+                                   const std::vector<AD2Type> &W) const
+{
+   return EvalW_AD_impl<AD2Type>(T,W);
+}
+
 real_t TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
 {
    // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
@@ -971,6 +1151,25 @@ void TMOP_Metric_014::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI1(weight, A.GetData());
 }
 
+template <typename type>
+type TMOP_Metric_014::EvalW_AD_impl(const std::vector<type> &T,
+                                    const std::vector<type> &W) const
+{
+   return mu14_ad(T, W);
+}
+
+AD1Type TMOP_Metric_014::EvalW_AD1(const std::vector<AD1Type> &T,
+                                   const std::vector<AD1Type> &W) const
+{
+   return EvalW_AD_impl<AD1Type>(T,W);
+}
+
+AD2Type TMOP_Metric_014::EvalW_AD2(const std::vector<AD2Type> &T,
+                                   const std::vector<AD2Type> &W) const
+{
+   return EvalW_AD_impl<AD2Type>(T,W);
+}
+
 real_t TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
 {
    // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
@@ -4096,7 +4295,7 @@ real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
          const IntegrationPoint &ip_s = ir_s->IntPoint(s);
          Tpr->SetIntPoint(&ip_s);
-         double w = surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
+         real_t w = surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
                     1.0 / surf_fit_dof_count[scalar_dof_id];
 
          if (surf_fit_gf)

diff --git a/fem/tmop.hpp b/fem/tmop.hpp
@@ -14,10 +14,14 @@
 
 #include "../linalg/invariants.hpp"
 #include "nonlininteg.hpp"
+#include "../linalg/dual.hpp"
 
 namespace mfem
 {
 
+using AD1Type = internal::dual<real_t, real_t>;
+using AD2Type = internal::dual<AD1Type, AD1Type>;
+
 /** @brief Abstract class for local mesh quality metrics in the target-matrix
     optimization paradigm (TMOP) by P. Knupp et al. */
 class TMOP_QualityMetric : public HyperelasticModel
@@ -69,6 +73,20 @@ class TMOP_QualityMetric : public HyperelasticModel
    virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
    { PW = 0.0;}
 
+   // First-derivative hook for AD-based computations:
+   virtual AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                             const std::vector<AD1Type> &W) const
+   {
+      return AD1Type{0.0, 0.0};
+   }
+
+   // Second-derivative hook for AD-based computations:
+   virtual AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                             const std::vector<AD2Type> &W) const
+   {
+      return AD2Type{{0.0, 0.0}, {0.0,0.0}};
+   }
+
    /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
        and assemble its contribution to the local gradient matrix 'A'.
        @param[in] Jpt     Represents the target->physical transformation
@@ -124,6 +142,12 @@ class TMOP_Combo_QualityMetric : public TMOP_QualityMetric
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
    /// Computes the averages of all metrics (integral of metric / volume).
    /// Works in parallel when called with a ParGridFunction.
    void ComputeAvgMetrics(const GridFunction &nodes,
@@ -221,12 +245,16 @@ class TMOP_WorstCaseUntangleOptimizer_Metric : public TMOP_QualityMetric
 
    real_t EvalW(const DenseMatrix &Jpt) const override;
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
 
    // Compute mu_hat.
    real_t EvalWBarrier(const DenseMatrix &Jpt) const;
@@ -368,6 +396,10 @@ class TMOP_Metric_004 : public TMOP_QualityMetric
 protected:
    mutable InvariantsEvaluator2D<real_t> ie;
 
+   template<typename type>
+   type EvalW_AD_impl(const std::vector<type> &T,
+                      const std::vector<type> &W) const;
+
 public:
    // W = |J|^2 - 2*det(J)
    real_t EvalW(const DenseMatrix &Jpt) const override;
@@ -377,6 +409,12 @@ class TMOP_Metric_004 : public TMOP_QualityMetric
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
    int Id() const override { return 4; }
 };
 
@@ -420,6 +458,10 @@ class TMOP_Metric_014 : public TMOP_QualityMetric
 protected:
    mutable InvariantsEvaluator2D<real_t> ie;
 
+   template <typename type>
+   type EvalW_AD_impl(const std::vector<type> &T,
+                      const std::vector<type> &W) const;
+
 public:
    // W = |J - I|^2.
    real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
@@ -431,6 +473,14 @@ class TMOP_Metric_014 : public TMOP_QualityMetric
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
+
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
+   int Id() const override { return 14; }
 };
 
 /// 2D Shifted barrier form of shape metric (mu_2).

diff --git a/miniapps/meshing/mesh-optimizer.cpp b/miniapps/meshing/mesh-optimizer.cpp
@@ -105,7 +105,7 @@
 //   2D untangling:
 //     mesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
 //   2D untangling with shifted barrier metric:
-//     mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
+//     mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -vl 1 -btype 1
 //   3D untangling (the mesh is in the mfem/data GitHub repository):
 //   * mesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
 

diff --git a/miniapps/meshing/pmesh-optimizer.cpp b/miniapps/meshing/pmesh-optimizer.cpp
@@ -107,7 +107,7 @@
 //   2D untangling:
 //     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
 //   2D untangling with shifted barrier metric:
-//     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
+//     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -vl 1 -btype 1
 //   3D untangling (the mesh is in the mfem/data GitHub repository):
 //   * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
 //   Shape optimization for a Kershaw transformed mesh using partial assembly: