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