Clean up code and comments, add test for InterfaceNormalCurvatures. Closes #32436

Author: Larry Aagesen

modules/phase_field/doc/content/source/materials/InterfaceNormalCurvatures.md

@@ -6,7 +6,7 @@
 
 `InterfaceNormalCurvatures` is a `Material` object that computes the two **normal curvatures** of a diffuse interface defined by an order parameter $\eta$. The interface geometry is encoded entirely in the gradient and Hessian of $\eta$; no explicit surface mesh is required.
 
-The two normal curvatures characterise how the interface bends in two orthogonal tangent directions:
+The two normal curvatures characterize how the interface bends in two orthogonal tangent directions:
 
 | Property | Symbol | Direction |
 |---|---|---|
@@ -121,9 +121,6 @@ The mean curvature $\kappa = \nabla \cdot \hat{n}$ is also declared as a materia
 | `kappa1` | `Real` | Normal curvature $\kappa_1$ along the in-plane tangent $\hat{t}_1$ |
 | `kappa2` | `Real` | Normal curvature $\kappa_2$ along the out-of-plane tangent $\hat{t}_2$ |
 | `kappa_mean` | `Real` | Mean curvature $\kappa = \kappa_1 + \kappa_2 = \nabla\cdot\hat{n}$ |
-| `interface_normal` | `RealVectorValue` | Unit interface normal $\hat{n}$ |
-| `tangent1` | `RealVectorValue` | First tangent vector $\hat{t}_1$ (in-plane) |
-| `tangent2` | `RealVectorValue` | Second tangent vector $\hat{t}_2$ (binormal) |
 
 All property names can be prefixed by setting `base_name`.
 

modules/phase_field/include/materials/InterfaceNormalCurvatures.h

@@ -13,25 +13,25 @@
 
 /**
  * Computes the two normal curvatures of a diffuse interface defined
- * by an order parameter η using n̂ = ∇η / |∇η|.
+ * by an order parameter eta using n = grad(eta) / | grad(eta) |.
  *
  * At every quadrature point the local frame is:
  *
- *   n̂  — unit interface normal   (= ∇η / |∇η|)
- *   t̂₁ — tangent in the xy-plane (= ẑ × n̂, normalised)
- *   t̂₂ — binormal                (= n̂ × t̂₁)
+ *   n  — unit interface normal   (grad(eta) / | grad(eta) |)
+ *   t_1 — tangent in the xy-plane (= z x n, normalized)
+ *   t_2 — binormal                (= n x t_1)
  *
  * The normal curvatures are then:
  *
- *   κ₁ = t̂₁ · S · t̂₁    (curvature along t̂₁, the in-plane tangent)
- *   κ₂ = t̂₂ · S · t̂₂    (curvature along t̂₂, the out-of-plane tangent)
+ *   kappa_1 = t_1 · S · t_1    (curvature along t_1, the in-plane tangent)
+ *   kappa_2 = t_2 · S · t_2    (curvature along t_2, the out-of-plane tangent)
  *
- * where  S = −∇n̂  is the shape operator (Weingarten map).
+ * where  S is the shape operator (Weingarten map).
  *
- * Note: κ₁ + κ₂ = mean curvature κ = ∇·n̂  (as a consistency check).
+ * Note: kappa_1 + kappa_2 = mean curvature kappa = div(n)  (as a consistency check).
  *
  * Requires second-order Lagrange (or higher) elements so that the
- * second derivatives of η are available.
+ * second derivatives of eta are available.
  */
 class InterfaceNormalCurvatures : public Material
 {
@@ -44,20 +44,17 @@ class InterfaceNormalCurvatures : public Material
 
 private:
   // ── order parameter ────────────────────────────────────────
-  const VariableValue  & _eta;        ///< η
-  const VariableGradient & _grad_eta; ///< ∇η
-  const VariableSecond & _second_eta; ///< ∇∇η  (Hessian)
+  const VariableValue & _eta;
+  const VariableGradient & _grad_eta;
+  const VariableSecond & _second_eta;
 
   // ── user parameters ────────────────────────────────────────
-  const Real _eps;   ///< regularisation floor for |∇η| (avoids /0)
+  const Real _eps; /// regularization floor for grad(eta) (avoids /0)
 
   // ── material properties produced ───────────────────────────
-  MaterialProperty<Real> & _kappa1; ///< normal curvature along t̂₁
-  MaterialProperty<Real> & _kappa2; ///< normal curvature along t̂₂
+  MaterialProperty<Real> & _kappa1; /// normal curvature along t_1
+  MaterialProperty<Real> & _kappa2; /// normal curvature along t_2
 
   // ── optional diagnostics ───────────────────────────────────
-  MaterialProperty<Real>         & _kappa_mean;  ///< κ₁+κ₂  (= ∇·n̂)
-  MaterialProperty<RealVectorValue> & _normal;   ///< n̂
-  MaterialProperty<RealVectorValue> & _tangent1; ///< t̂₁
-  MaterialProperty<RealVectorValue> & _tangent2; ///< t̂₂
+  MaterialProperty<Real> & _kappa_mean;
 };

modules/phase_field/src/materials/InterfaceNormalCurvatures.C

@@ -17,19 +17,14 @@ InterfaceNormalCurvatures::validParams()
 {
   InputParameters params = Material::validParams();
   params.addClassDescription(
-      "Computes the two normal curvatures κ₁ and κ₂ of a diffuse interface "
-      "from an order parameter η via n̂ = ∇η/|∇η| and the shape operator "
-      "S = −∇n̂.  t̂₁ lies in the xy-plane (ẑ×n̂); t̂₂ = n̂×t̂₁.");
-
-  params.addRequiredCoupledVar("eta", "Order parameter η that defines the interface");
-
-  params.addParam<Real>("regularization", 1e-8,
-      "Floor added to |∇η| before division to prevent singularities in "
-      "bulk regions where ∇η ≈ 0");
-
-  params.addParam<std::string>("base_name", "",
-      "Optional prefix for all material property names");
-
+      "Computes the two normal curvatures kappa_1 and kappa_2 of a diffuse interface "
+      "from an order parameter eta.");
+  params.addRequiredCoupledVar("eta", "Order parameter that defines the interface");
+  params.addParam<Real>("regularization",
+                        1e-8,
+                        "Floor added to |grad(eta)| before division to prevent singularities in "
+                        "bulk regions");
+  params.addParam<std::string>("base_name", "", "Optional prefix for all material property names");
   return params;
 }
 
@@ -39,13 +34,9 @@ InterfaceNormalCurvatures::InterfaceNormalCurvatures(const InputParameters & par
     _grad_eta(coupledGradient("eta")),
     _second_eta(coupledSecond("eta")),
     _eps(getParam<Real>("regularization")),
-
-    _kappa1   (declareProperty<Real>            (getParam<std::string>("base_name") + "kappa1")),
-    _kappa2   (declareProperty<Real>            (getParam<std::string>("base_name") + "kappa2")),
-    _kappa_mean(declareProperty<Real>           (getParam<std::string>("base_name") + "kappa_mean")),
-    _normal   (declareProperty<RealVectorValue> (getParam<std::string>("base_name") + "interface_normal")),
-    _tangent1 (declareProperty<RealVectorValue> (getParam<std::string>("base_name") + "tangent1")),
-    _tangent2 (declareProperty<RealVectorValue> (getParam<std::string>("base_name") + "tangent2"))
+    _kappa1(declareProperty<Real>(getParam<std::string>("base_name") + "kappa1")),
+    _kappa2(declareProperty<Real>(getParam<std::string>("base_name") + "kappa2")),
+    _kappa_mean(declareProperty<Real>(getParam<std::string>("base_name") + "kappa_mean"))
 {
 }
 
@@ -54,58 +45,34 @@ InterfaceNormalCurvatures::computeQpProperties()
 {
   // ── 1.  Interface normal ─────────────────────────────────────────────────
   const RealVectorValue & g = _grad_eta[_qp];
-  const Real              g_mag = g.norm();
-  const Real              g_reg = g_mag + _eps;
+  const Real g_mag = g.norm();
+  const Real g_reg = g_mag + _eps;
 
-  const RealVectorValue nhat = g / g_reg;   // n̂  (unit normal)
-  _normal[_qp] = nhat;
+  const RealVectorValue nhat = g / g_reg; // n  (unit normal)
 
   // ── 2.  Tangent frame ────────────────────────────────────────────────────
-  // t̂₁ = ẑ × n̂  projected into the xy-plane.
-  // For a normal that is purely along ẑ the cross product vanishes;
-  // in that case we fall back to x̂ so the frame is always well-defined.
   static const RealVectorValue ZHAT(0., 0., 1.);
   static const RealVectorValue XHAT(1., 0., 0.);
 
   RealVectorValue t1 = ZHAT.cross(nhat);
   const Real t1_mag = t1.norm();
   if (t1_mag < 1e-12)
-    t1 = XHAT;   // n̂ ≈ ±ẑ: choose x̂ as in-plane tangent
+    t1 = XHAT; // n ≈ +/- z: choose x as in-plane tangent
   else
     t1 /= t1_mag;
 
-  const RealVectorValue t2 = nhat.cross(t1);   // already unit length
-
-  _tangent1[_qp] = t1;
-  _tangent2[_qp] = t2;
+  const RealVectorValue t2 = nhat.cross(t1); // already unit length
 
-  // ── 3.  Shape operator  S = −∇n̂ ─────────────────────────────────────────
-  //
-  // n̂ = ∇η / |∇η|,  so its Jacobian is:
-  //
-  //   ∂n̂ᵢ/∂xⱼ = ( Hᵢⱼ  −  nᵢ (∇η · ∇∂η/∂xⱼ) / |∇η| ) / |∇η|
-  //
-  // where H = ∇∇η is the Hessian of η.
-  // Written compactly:  ∇n̂ = (H − n̂ ⊗ (H·n̂)) / |∇η|
-  // and the shape operator  S = −∇n̂.
-  //
-  // The normal curvature in direction v̂ is  κ(v̂) = v̂ · S · v̂.
+  // ── 3.  Shape operator  S ─────────────────────────────────────────
 
-  const RankTwoTensor & H = _second_eta[_qp];  // Hessian ∇∇η  (3×3)
+  const RankTwoTensor & H = _second_eta[_qp]; // Hessian (3×3)
 
-  // H·n̂  (matrix-vector product)
   RealVectorValue Hn;
   for (unsigned i = 0; i < 3; ++i)
     for (unsigned j = 0; j < 3; ++j)
       Hn(i) += H(i, j) * nhat(j);
 
-  // Normal curvature along a unit tangent v̂:
-  //   κ(v̂) = v̂ · S · v̂
-  //         = −v̂ · (∇n̂) · v̂
-  //         = −[ (v̂·H·v̂) − (v̂·n̂)(n̂·H·n̂) ] / |∇η|
-  // but since v̂ ⊥ n̂  →  v̂·n̂ = 0, this simplifies to:
-  //   κ(v̂) = −(v̂·H·v̂) / |∇η|                               (*)
-
+  // Normal curvature along a unit tangent v:
   auto vHv = [&](const RealVectorValue & v) -> Real
   {
     Real val = 0.;
@@ -115,11 +82,11 @@ InterfaceNormalCurvatures::computeQpProperties()
     return val;
   };
 
-  _kappa1[_qp] = -vHv(t1) / g_reg;   // κ₁  (in-plane tangent)
-  _kappa2[_qp] = -vHv(t2) / g_reg;   // κ₂  (out-of-plane tangent)
+  _kappa1[_qp] = -vHv(t1) / g_reg; // in-plane tangent
+  _kappa2[_qp] = -vHv(t2) / g_reg; // out-of-plane tangent
+
+  // ── 4.  Mean curvature ───────────────────────────────────────────
 
-  // ── 4.  Mean curvature (diagnostic)  κ = κ₁ + κ₂ = ∇·n̂ ─────────────────
-  // Trace of shape operator = H.trace()/|∇η| − (n̂·H·n̂)/|∇η|
-  const Real nHn = nhat * Hn;          // n̂·H·n̂
+  const Real nHn = nhat * Hn;
   _kappa_mean[_qp] = -(H.tr() - nHn) / g_reg;
 }

modules/phase_field/test/tests/misc/gold/interface_normal_curvatures_out.e

@@ -0,0 +1,81 @@
+#
+# Test calculation of normal curvatures of an interface defined by order parameter c
+#
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  nx = 8
+  ny = 8
+  nz = 8
+  xmin = 0
+  xmax = 20
+  ymin = -10
+  ymax = 10
+  zmin = 0
+  zmax = 20
+[]
+
+[Modules]
+  [PhaseField]
+    [Conserved]
+      [./c]
+        solve_type = direct
+        free_energy = F
+        kappa = 2.0
+        mobility = 1.0
+      []
+    []
+  []
+[]
+
+[ICs]
+  [InitialCondition]
+    type = FunctionIC
+    function = ic_func_c
+    variable = c
+  []
+[]
+
+[Functions]
+  [ic_func_c]
+    type = ParsedFunction
+    symbol_names = ' A   lambda '
+    symbol_values = '0.5 40     '
+    expression = '0.5*(1.0-tanh((y - A * sin(2*pi*x/lambda) * sin(2*pi*z/lambda)) / 2))'
+  []
+[]
+
+[Materials]
+  [free_energy]
+    type = DerivativeParsedMaterial
+    property_name = F
+    coupled_variables = 'c'
+    expression = 'c^2 * (1 - c)^2'
+  []
+  [curvatures]
+    type = InterfaceNormalCurvatures
+    eta = c
+    outputs = 'exodus'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  scheme = 'bdf2'
+
+  solve_type = 'NEWTON'
+
+  l_max_its = 30
+  l_tol = 1.0e-4
+  nl_max_its = 10
+  nl_rel_tol = 1.0e-9
+  nl_abs_tol = 1.0e-11
+
+  start_time = 0.0
+  end_time = 0.001
+  dt = 0.001
+[]
+
+[Outputs]
+  exodus = true
+[]

modules/phase_field/test/tests/misc/interface_normal_curvatures.i

@@ -56,4 +56,15 @@
                   'the value of a function evaluated over a set of up to four coupled variables'
     design = 'CoupledValueFunctionIC.md'
   []
+
+  [interface_normal_curvatures]
+    type = 'Exodiff'
+    input = 'interface_normal_curvatures.i'
+    exodiff = 'interface_normal_curvatures_out.e'
+    issues = '32436'
+    requirement = 'An Material shall be implemented that can calculate the normal curvatures of'
+                  'an interface defined by an order parameter'
+    design = 'InterfaceNormalCurvatures.md'
+  []
+
 []