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 "Density.h"
      11             : 
      12             : registerMooseObjectDeprecated("MiscApp", Density, "12/31/2025 24:00");
      13             : registerMooseObjectDeprecated("MiscApp", ADDensity, "12/31/2025 24:00");
      14             : 
      15             : template <bool is_ad>
      16             : InputParameters
      17           0 : DensityTempl<is_ad>::validParams()
      18             : {
      19           0 :   InputParameters params = Material::validParams();
      20             : 
      21           0 :   params.addClassDescription(
      22             :       "Creates density material property. This class is deprecated, and its functionality"
      23             :       "is replaced by StrainAdjustedDensity for cases when the density should be adjusted"
      24             :       "to account for material deformation. If it is not desired to adjust the density for"
      25             :       "deformation, a variety of general-purpose Materials, such as GenericConstantMaterial"
      26             :       "or ParsedMaterial can be used to define the density. A StrainAdjustedDensity can also "
      27             :       "be used with '0 0 0' for the displacements.");
      28             : 
      29           0 :   params.addCoupledVar(
      30             :       "displacements",
      31             :       "The displacements appropriate for the simulation geometry and coordinate system");
      32             : 
      33           0 :   params.addParam<std::string>("base_name",
      34             :                                "Optional parameter that allows the user to define "
      35             :                                "multiple material systems on the same block, "
      36             :                                "e.g. for multiple phases");
      37           0 :   params.addRequiredParam<Real>("density", "Density");
      38             : 
      39           0 :   return params;
      40           0 : }
      41             : 
      42             : template <bool is_ad>
      43           0 : DensityTempl<is_ad>::DensityTempl(const InputParameters & parameters)
      44             :   : Material(parameters),
      45           0 :     _is_coupled(isCoupled("displacements")),
      46           0 :     _coord_system(getBlockCoordSystem()),
      47           0 :     _disp_r(_is_coupled ? this->template coupledGenericValue<is_ad>("displacements", 0)
      48           0 :                         : genericZeroValue<is_ad>()),
      49           0 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      50           0 :     _initial_density(getParam<Real>("density")),
      51           0 :     _grad_disp(_is_coupled ? this->template coupledGenericGradients<is_ad>("displacements")
      52             :                            : std::vector<const GenericVariableGradient<is_ad> *>()),
      53           0 :     _density(declareGenericProperty<Real, is_ad>(_base_name + "density"))
      54             : {
      55           0 :   if (getParam<bool>("use_displaced_mesh"))
      56           0 :     paramError("use_displaced_mesh",
      57             :                "Density needs to act on an undisplaced mesh. Use of a displaced mesh leads to "
      58             :                "incorrect gradient values");
      59             : 
      60             :   // get coupled gradients
      61           0 :   const unsigned int ndisp = coupledComponents("displacements");
      62             : 
      63           0 :   if (ndisp == 0 && _fe_problem.getDisplacedProblem())
      64           0 :     paramError(
      65             :         "displacements",
      66             :         "The system uses a displaced problem but 'displacements' are not provided in Density.");
      67             : 
      68             :   // fill remaining components with zero
      69           0 :   _grad_disp.resize(3, &genericZeroGradient<is_ad>());
      70           0 : }
      71             : 
      72             : template <bool is_ad>
      73             : void
      74           0 : DensityTempl<is_ad>::initQpStatefulProperties()
      75             : {
      76           0 :   _density[_qp] = _initial_density;
      77           0 : }
      78             : 
      79             : template <bool is_ad>
      80             : void
      81           0 : DensityTempl<is_ad>::computeQpProperties()
      82             : {
      83           0 :   GenericReal<is_ad> density = _initial_density;
      84             : 
      85             :   // TODO: We should deprecate the !_is_coupled case and have the
      86             :   // user use a GenericConstantMaterial
      87           0 :   if (_is_coupled)
      88             :   {
      89             :     // rho * V = rho0 * V0
      90             :     // rho = rho0 * V0 / V
      91             :     // rho = rho0 / det(F)
      92             :     // rho = rho0 / det(grad(u) + 1)
      93             : 
      94           0 :     const auto Axx = (*_grad_disp[0])[_qp](0) + 1.0;
      95             :     const auto & Axy = (*_grad_disp[0])[_qp](1);
      96             :     const auto & Axz = (*_grad_disp[0])[_qp](2);
      97           0 :     const auto & Ayx = (*_grad_disp[1])[_qp](0);
      98           0 :     auto Ayy = (*_grad_disp[1])[_qp](1) + 1.0;
      99             :     const auto & Ayz = (*_grad_disp[1])[_qp](2);
     100           0 :     const auto & Azx = (*_grad_disp[2])[_qp](0);
     101             :     const auto & Azy = (*_grad_disp[2])[_qp](1);
     102           0 :     auto Azz = (*_grad_disp[2])[_qp](2) + 1.0;
     103             : 
     104           0 :     switch (_coord_system)
     105             :     {
     106           0 :       case Moose::COORD_XYZ:
     107           0 :         Azz = (*_grad_disp[2])[_qp](2) + 1.0;
     108           0 :         break;
     109             : 
     110           0 :       case Moose::COORD_RZ:
     111           0 :         if (_q_point[_qp](0) != 0.0)
     112           0 :           Azz = _disp_r[_qp] / _q_point[_qp](0) + 1.0;
     113             :         break;
     114             : 
     115           0 :       case Moose::COORD_RSPHERICAL:
     116           0 :         if (_q_point[_qp](0) != 0.0)
     117           0 :           Ayy = Azz = _disp_r[_qp] / _q_point[_qp](0) + 1.0;
     118             :         break;
     119             :     }
     120             : 
     121           0 :     const auto detF = Axx * Ayy * Azz + Axy * Ayz * Azx + Axz * Ayx * Azy - Azx * Ayy * Axz -
     122           0 :                       Azy * Ayz * Axx - Azz * Ayx * Axy;
     123           0 :     density /= detF;
     124             :   }
     125             : 
     126           0 :   _density[_qp] = density;
     127           0 : }
     128             : 
     129             : template class DensityTempl<false>;
     130             : template class DensityTempl<true>;