Axisymmetric geometry support in NLFormIntegrators

kernel/include/Coefficients/AxiCylindricalCoefficient.hpp

@@ -0,0 +1,43 @@
+/**
+ * @file AxiCylindricalCoefficient.hpp
+ * @author Clément Introïni (clement.introini@cea.fr)
+ * @brief Built a MFEM Coefficient for AxiCylindricalCoefficient
+ * @version 0.1
+ * @date 2026-06-16
+ *
+ * @copyright CEA (C) 2026
+ *
+ * This file is part of SLOTH.
+ *
+ * SLOTH is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * SLOTH is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+#include <cmath>
+#include <cstdlib>
+#include <span>
+#include <vector>
+
+#include "Coefficients/Coefficients.hpp"
+#include "mfem.hpp"  // NOLINT [no include the directory when naming mfem include file]
+
+class AxiCylindricalCoefficient : public mfem::Coefficient {
+ private:
+  mutable mfem::Vector transip;
+
+ public:
+  AxiCylindricalCoefficient();
+
+  mfem::real_t Eval(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip) override;
+};

kernel/include/Coefficients/Coefficients.hpp

@@ -62,8 +62,9 @@ class Coefficients {
 
   void add(Coefficient coef);
   std::vector<Coefficient> getCoefficients() const;
-  size_t size() noexcept;
-  Coefficient operator[](size_t i);
+  size_t size() const noexcept;
+  Coefficient& operator[](size_t i);
+  const Coefficient& operator[](size_t i) const;
 
   std::vector<GlossaryType> get_types() const;
 

kernel/include/Coefficients/ListCoefficients.hpp

@@ -25,6 +25,7 @@
  */
 #pragma once
 
+#include "Coefficients/AxiCylindricalCoefficient.hpp"
 #include "Coefficients/CommonCoefficients.hpp"
 #include "Coefficients/ProductFunctionCoefficients.hpp"
 #include "Coefficients/SumFunctionCoefficients.hpp"

kernel/include/Coefficients/SlothBaseCoefficient.hpp

@@ -117,5 +117,5 @@ class SlothBaseCoefficient {
 
   void set_time(double time);
   void set_bdr_index_coef(std::vector<int> ids);
-  std::vector<int> get_bdr_index_coef();
+  std::vector<int> get_bdr_index_coef() const;
 };

kernel/include/Integrators/AllenCahnNLFormIntegrator.hpp

@@ -72,7 +72,7 @@ class AllenCahnNLFormIntegrator : public SlothNLFormIntegrator<VARS> {
 
  public:
   void init() override;
-  AllenCahnNLFormIntegrator(const std::vector<mfem::ParGridFunction>& u_old,
+  AllenCahnNLFormIntegrator(Geometry geometry, const std::vector<mfem::ParGridFunction>& u_old,
                             const std::vector<mfem::ParGridFunction>& aux_old,
                             const Parameters& params, std::vector<VARS*> auxvars,
                             const std::vector<Coefficients>& coefficients);

kernel/include/Integrators/AllenCahnNLFormIntegrator.tpp

@@ -33,6 +33,7 @@
 #include <utility>
 #include <vector>
 
+#include "Coefficients/AxiCylindricalCoefficient.hpp"
 #include "Coefficients/SlothBaseCoefficient.hpp"
 #include "Integrators/SlothNLFormIntegrator.hpp"
 #include "MAToolsProfiling/MATimersAPI.hxx"
@@ -60,10 +61,10 @@
  */
 template <class VARS>
 AllenCahnNLFormIntegrator<VARS>::AllenCahnNLFormIntegrator(
-    const std::vector<mfem::ParGridFunction>& u_old,
+    Geometry geometry, const std::vector<mfem::ParGridFunction>& u_old,
     const std::vector<mfem::ParGridFunction>& aux_old, const Parameters& params,
     std::vector<VARS*> auxvars, const std::vector<Coefficients>& coefficients)
-    : SlothNLFormIntegrator<VARS>(u_old, aux_old, params, auxvars, coefficients) {
+    : SlothNLFormIntegrator<VARS>(geometry, u_old, aux_old, params, auxvars, coefficients) {
   this->integrator_name_ = "AllenCahn";
 
   this->check_variables_consistency();
@@ -167,15 +168,21 @@ void AllenCahnNLFormIntegrator<VARS>::AssembleElementVector(
                                              std::span<const double>(vaux_gf_at_ip));
       const double lamb = this->compute_coefficient(lambda[blk], std::span<const double>(u_values),
                                                     std::span<const double>(vaux_gf_at_ip));
-      gradU *= coef_mobi * lamb * ip.weight * Tr.Weight();
-      ;
+
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
+      gradU *= coef_mobi * lamb * weight_coef;
+
       gradPsi.AddMult(gradU, *elvect[blk]);
 
       // Given u, compute (w'(u), psi), psi is shape function
       double ww = coef_mobi * this->compute_gradient_energy_coefficient(
                                   double_well_energy[blk], blk, std::span<const double>(u_values),
                                   std::span<const double>(vaux_gf_at_ip));
-      ww *= ip.weight * Tr.Weight();
+      ww *= weight_coef;
       add(*elvect[blk], ww, Psi, *elvect[blk]);
     }
   }
@@ -246,15 +253,20 @@ void AllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
       const double lamb = this->compute_coefficient(lambda[blk], std::span<const double>(u_values),
                                                     std::span<const double>(vaux_gf_at_ip));
 
-      const double mob_lamb = coef_mobi * lamb * ip.weight * Tr.Weight();
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
+      const double mob_lamb = coef_mobi * lamb * weight_coef;
       AddMult_a_AAt(mob_lamb, gradPsi, *elmats(blk, blk));
 
       // Compute w'(u)*(du,psi), psi is shape function ( // w''(u))
       double fun_val =
           coef_mobi * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk,
                                                         std::span<const double>(u_values),
                                                         std::span<const double>(vaux_gf_at_ip));
-      fun_val *= ip.weight * Tr.Weight();
+      fun_val *= weight_coef;
 
       AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk));  // w'(u)*(du, psi)
     }
@@ -286,7 +298,12 @@ void AllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
             coef_mobi * this->compute_hessian_coefficient(double_well_energy[jblk], blk, jblk,
                                                           std::span<const double>(u_values),
                                                           std::span<const double>(vaux_gf_at_ip));
-        fun_val *= ip.weight * Tr.Weight();
+        double weight_coef = ip.weight * Tr.Weight();
+        if (this->isAxisymmetric()) {
+          weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+        }
+
+        fun_val *= weight_coef;
         AddMult_a_VVt(fun_val, Psi, *elmats(blk, jblk));
       }
     }
@@ -316,7 +333,13 @@ void AllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
             coef_mobi * this->compute_hessian_coefficient(double_well_energy[jblk], blk, jblk,
                                                           std::span<const double>(u_values),
                                                           std::span<const double>(vaux_gf_at_ip));
-        fun_val *= ip.weight * Tr.Weight();
+
+        double weight_coef = ip.weight * Tr.Weight();
+        if (this->isAxisymmetric()) {
+          weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+        }
+
+        fun_val *= weight_coef;
         AddMult_a_VVt(fun_val, Psi, *elmats(blk, jblk));
       }
     }

kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.hpp

@@ -65,7 +65,7 @@ class BlockAllenCahnNLFormIntegrator : public SlothNLFormIntegrator<VARS> {
 
  public:
   void init() override;
-  BlockAllenCahnNLFormIntegrator(const std::vector<mfem::ParGridFunction>& u_old,
+  BlockAllenCahnNLFormIntegrator(Geometry geometry, const std::vector<mfem::ParGridFunction>& u_old,
                                  const std::vector<mfem::ParGridFunction>& aux_old,
                                  const Parameters& params, std::vector<VARS*> auxvars,
                                  const std::vector<Coefficients>& coefficients);

kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.tpp

@@ -32,6 +32,7 @@
 #include <utility>
 #include <vector>
 
+#include "Coefficients/AxiCylindricalCoefficient.hpp"
 #include "Coefficients/SlothBaseCoefficient.hpp"
 #include "Integrators/SlothGridFunction.hpp"
 #include "Integrators/SlothNLFormIntegrator.hpp"
@@ -59,10 +60,10 @@
  */
 template <class VARS>
 BlockAllenCahnNLFormIntegrator<VARS>::BlockAllenCahnNLFormIntegrator(
-    const std::vector<mfem::ParGridFunction>& u_old,
+    Geometry geometry, const std::vector<mfem::ParGridFunction>& u_old,
     const std::vector<mfem::ParGridFunction>& aux_old, const Parameters& params,
     std::vector<VARS*> auxvars, const std::vector<Coefficients>& coefficients)
-    : SlothNLFormIntegrator<VARS>(u_old, aux_old, params, auxvars, coefficients) {
+    : SlothNLFormIntegrator<VARS>(geometry, u_old, aux_old, params, auxvars, coefficients) {
   this->check_variables_consistency();
 }
 
@@ -163,18 +164,23 @@ void BlockAllenCahnNLFormIntegrator<VARS>::AssembleElementVector(
         u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip);
       }
 
-      const double xx = ip.weight * Tr.Weight();
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
       el[blk]->CalcPhysDShape(Tr, gradPsi);
       gradPsi.MultTranspose(*elfun[blk], gradU);
 
-      gradU *= xx * this->compute_coefficient(lambda[blk], std::span<const double>(u_values),
-                                              std::span<const double>(vaux_gf_at_ip));
+      gradU *=
+          weight_coef * this->compute_coefficient(lambda[blk], std::span<const double>(u_values),
+                                                  std::span<const double>(vaux_gf_at_ip));
       gradPsi.AddMult(gradU, *elvect[blk]);
 
       const double ww =
-          xx * (eta + this->compute_gradient_energy_coefficient(
-                          double_well_energy[blk], blk, std::span<const double>(u_values),
-                          std::span<const double>(vaux_gf_at_ip)));
+          weight_coef * (eta + this->compute_gradient_energy_coefficient(
+                                   double_well_energy[blk], blk, std::span<const double>(u_values),
+                                   std::span<const double>(vaux_gf_at_ip)));
 
       add(*elvect[blk], ww, Psi, *elvect[blk]);
     }
@@ -212,10 +218,15 @@ void BlockAllenCahnNLFormIntegrator<VARS>::AssembleElementVector(
         u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip);
       }
 
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
       const double coef_mob =
           this->compute_coefficient(mobility[blk], std::span<const double>(u_values),
                                     std::span<const double>(vaux_gf_at_ip)) *
-          ip.weight * Tr.Weight();
+          weight_coef;
 
       const double ww = -coef_mob * eta;
 
@@ -281,18 +292,22 @@ void BlockAllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
         u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip);
       }
 
-      const double xx = ip.weight * Tr.Weight();
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
       el[blk]->CalcPhysDShape(Tr, gradPsi);
 
       const double coef_lambda = this->compute_coefficient(
           lambda[blk], std::span<const double>(u_values), std::span<const double>(vaux_gf_at_ip));
 
-      AddMult_a_AAt(xx * coef_lambda, gradPsi, *elmats(blk, blk));
+      AddMult_a_AAt(weight_coef * coef_lambda, gradPsi, *elmats(blk, blk));
 
       double fun_val =
-          xx * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk,
-                                                 std::span<const double>(u_values),
-                                                 std::span<const double>(vaux_gf_at_ip));
+          weight_coef * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk,
+                                                          std::span<const double>(u_values),
+                                                          std::span<const double>(vaux_gf_at_ip));
       AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk));
     }
   }
@@ -318,8 +333,11 @@ void BlockAllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
       el[blk]->CalcShape(ip, Psi);
       Tr.SetIntPoint(&ip);
 
-      double ww = ip.weight * Tr.Weight();
-      AddMult_a_VVt(ww, Psi, *elmats(blk, off_blk));
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+      AddMult_a_VVt(weight_coef, Psi, *elmats(blk, off_blk));
     }
   }
   // Block 1 0  dR(eta)dphi=d(M eta )/dphi
@@ -368,10 +386,16 @@ void BlockAllenCahnNLFormIntegrator<VARS>::AssembleElementGrad(
         u_values[off_blk] = (*elfun[off_blk]) * Psi;
         u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip);
       }
+
+      double weight_coef = ip.weight * Tr.Weight();
+      if (this->isAxisymmetric()) {
+        weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip);
+      }
+
       const double coef_mob =
           -this->compute_coefficient(mobility[blk], std::span<const double>(u_values),
                                      std::span<const double>(vaux_gf_at_ip)) *
-          ip.weight * Tr.Weight();
+          weight_coef;
 
       AddMult_a_VVt(coef_mob, Psi, *elmats(blk, blk));
     }

kernel/include/Integrators/CahnHilliardNLFormIntegrator.hpp

@@ -66,7 +66,7 @@ class CahnHilliardNLFormIntegrator : public SlothNLFormIntegrator<VARS> {
 
  public:
   void init() override;
-  CahnHilliardNLFormIntegrator(const std::vector<mfem::ParGridFunction>& u_old,
+  CahnHilliardNLFormIntegrator(Geometry geometry, const std::vector<mfem::ParGridFunction>& u_old,
                                const std::vector<mfem::ParGridFunction>& aux_old,
                                const Parameters& params, std::vector<VARS*> auxvars,
                                const std::vector<Coefficients>& coefficients);