Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "InterfaceNormalCurvatures.h"
      11             : #include "Assembly.h"
      12             : 
      13             : registerMooseObject("PhaseFieldApp", InterfaceNormalCurvatures);
      14             : 
      15             : InputParameters
      16          31 : InterfaceNormalCurvatures::validParams()
      17             : {
      18          31 :   InputParameters params = Material::validParams();
      19          31 :   params.addClassDescription(
      20             :       "Computes the two normal curvatures kappa_1 and kappa_2 of a diffuse interface "
      21             :       "from an order parameter eta.");
      22          62 :   params.addRequiredCoupledVar("eta", "Order parameter that defines the interface");
      23          62 :   params.addParam<Real>(
      24             :       "gradient_threshold",
      25          62 :       1e-6,
      26             :       "If |grad(eta)| is less than this threshold at a point, curvatures are set to zero there.");
      27          62 :   params.addParam<std::string>("base_name", "", "Optional prefix for all material property names");
      28          31 :   return params;
      29           0 : }
      30             : 
      31          24 : InterfaceNormalCurvatures::InterfaceNormalCurvatures(const InputParameters & params)
      32             :   : Material(params),
      33          24 :     _eta(coupledValue("eta")),
      34          24 :     _grad_eta(coupledGradient("eta")),
      35          24 :     _second_eta(coupledSecond("eta")),
      36          48 :     _grad_threshold(getParam<Real>("gradient_threshold")),
      37          48 :     _kappa1(declareProperty<Real>(getParam<std::string>("base_name") + "kappa1")),
      38          48 :     _kappa2(declareProperty<Real>(getParam<std::string>("base_name") + "kappa2")),
      39          72 :     _kappa_mean(declareProperty<Real>(getParam<std::string>("base_name") + "kappa_mean"))
      40             : {
      41          24 : }
      42             : 
      43             : void
      44      327680 : InterfaceNormalCurvatures::computeQpProperties()
      45             : {
      46             :   // -- 1.  Interface normal -------------------------------------------------
      47      327680 :   const RealVectorValue & g = _grad_eta[_qp];
      48      327680 :   const Real g_mag = g.norm();
      49      327680 :   if (g_mag < _grad_threshold)
      50           0 :     _kappa1[_qp] = _kappa2[_qp] = _kappa_mean[_qp] = 0;
      51             :   else
      52             :   {
      53             :     const RealVectorValue nhat = g / g_mag; // n  (unit normal)
      54             : 
      55             :     // -- 2.  Tangent frame ----------------------------------------------------
      56      327680 :     static const RealVectorValue ZHAT(0., 0., 1.);
      57      327680 :     static const RealVectorValue XHAT(1., 0., 0.);
      58             : 
      59      327680 :     RealVectorValue t1 = ZHAT.cross(nhat); // chosen to lie in xy-plane
      60      327680 :     const Real t1_mag = t1.norm();
      61      327680 :     if (t1_mag < 1e-12)
      62           0 :       t1 = XHAT; // if n = +/- z, choose x to be in-plane tangent t_1
      63             :     else
      64             :       t1 /= t1_mag;
      65             : 
      66      327680 :     const RealVectorValue t2 = nhat.cross(t1); // already unit length
      67             : 
      68             :     // -- 3.  Shape operator  S ------------------------------------------------
      69             : 
      70      327680 :     const RankTwoTensor & H = _second_eta[_qp]; // Hessian (3 x 3)
      71             : 
      72             :     RealVectorValue Hn;
      73     1310720 :     for (unsigned i = 0; i < 3; ++i)
      74     3932160 :       for (unsigned j = 0; j < 3; ++j)
      75     2949120 :         Hn(i) += H(i, j) * nhat(j);
      76             : 
      77             :     // Normal curvature along a unit tangent v:
      78      655360 :     auto vHv = [&](const RealVectorValue & v) -> Real
      79             :     {
      80             :       Real val = 0.;
      81     2621440 :       for (unsigned i = 0; i < 3; ++i)
      82     7864320 :         for (unsigned j = 0; j < 3; ++j)
      83     5898240 :           val += v(i) * H(i, j) * v(j);
      84      655360 :       return val;
      85      327680 :     };
      86             : 
      87      327680 :     _kappa1[_qp] = -vHv(t1) / g_mag; // curvature in t_1 direction
      88      327680 :     _kappa2[_qp] = -vHv(t2) / g_mag; // curvature in t_2 direction
      89             : 
      90             :     // -- 4.  Mean curvature ---------------------------------------------------
      91             : 
      92             :     const Real nHn = nhat * Hn;
      93      327680 :     _kappa_mean[_qp] = -(H.tr() - nHn) / g_mag / 2;
      94             :   }
      95      327680 : }