Density.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "Density.h"

registerMooseObjectDeprecated("MiscApp", Density, "12/31/2025 24:00");
registerMooseObjectDeprecated("MiscApp", ADDensity, "12/31/2025 24:00");

template <bool is_ad>
InputParameters
DensityTempl<is_ad>::validParams()
{
  InputParameters params = Material::validParams();

  params.addClassDescription(
      "Creates density material property. This class is deprecated, and its functionality"
      "is replaced by StrainAdjustedDensity for cases when the density should be adjusted"
      "to account for material deformation. If it is not desired to adjust the density for"
      "deformation, a variety of general-purpose Materials, such as GenericConstantMaterial"
      "or ParsedMaterial can be used to define the density. A StrainAdjustedDensity can also "
      "be used with '0 0 0' for the displacements.");

  params.addCoupledVar(
      "displacements",
      "The displacements appropriate for the simulation geometry and coordinate system");

  params.addParam<std::string>("base_name",
                               "Optional parameter that allows the user to define "
                               "multiple material systems on the same block, "
                               "e.g. for multiple phases");
  params.addRequiredParam<Real>("density", "Density");

  return params;
}

template <bool is_ad>
DensityTempl<is_ad>::DensityTempl(const InputParameters & parameters)
  : Material(parameters),
    _is_coupled(isCoupled("displacements")),
    _coord_system(getBlockCoordSystem()),
    _disp_r(_is_coupled ? this->template coupledGenericValue<is_ad>("displacements", 0)
                        : genericZeroValue<is_ad>()),
    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
    _initial_density(getParam<Real>("density")),
    _grad_disp(_is_coupled ? this->template coupledGenericGradients<is_ad>("displacements")
                           : std::vector<const GenericVariableGradient<is_ad> *>()),
    _density(declareGenericProperty<Real, is_ad>(_base_name + "density"))
{
  if (getParam<bool>("use_displaced_mesh"))
    paramError("use_displaced_mesh",
               "Density needs to act on an undisplaced mesh. Use of a displaced mesh leads to "
               "incorrect gradient values");

  // get coupled gradients
  const unsigned int ndisp = coupledComponents("displacements");

  if (ndisp == 0 && _fe_problem.getDisplacedProblem())
    paramError(
        "displacements",
        "The system uses a displaced problem but 'displacements' are not provided in Density.");

  // fill remaining components with zero
  _grad_disp.resize(3, &genericZeroGradient<is_ad>());
}

template <bool is_ad>
void
DensityTempl<is_ad>::initQpStatefulProperties()
{
  _density[_qp] = _initial_density;
}

template <bool is_ad>
void
DensityTempl<is_ad>::computeQpProperties()
{
  GenericReal<is_ad> density = _initial_density;

  // TODO: We should deprecate the !_is_coupled case and have the
  // user use a GenericConstantMaterial
  if (_is_coupled)
  {
    // rho * V = rho0 * V0
    // rho = rho0 * V0 / V
    // rho = rho0 / det(F)
    // rho = rho0 / det(grad(u) + 1)

    const auto Axx = (*_grad_disp[0])[_qp](0) + 1.0;
    const auto & Axy = (*_grad_disp[0])[_qp](1);
    const auto & Axz = (*_grad_disp[0])[_qp](2);
    const auto & Ayx = (*_grad_disp[1])[_qp](0);
    auto Ayy = (*_grad_disp[1])[_qp](1) + 1.0;
    const auto & Ayz = (*_grad_disp[1])[_qp](2);
    const auto & Azx = (*_grad_disp[2])[_qp](0);
    const auto & Azy = (*_grad_disp[2])[_qp](1);
    auto Azz = (*_grad_disp[2])[_qp](2) + 1.0;

    switch (_coord_system)
    {
      case Moose::COORD_XYZ:
        Azz = (*_grad_disp[2])[_qp](2) + 1.0;
        break;

      case Moose::COORD_RZ:
        if (_q_point[_qp](0) != 0.0)
          Azz = _disp_r[_qp] / _q_point[_qp](0) + 1.0;
        break;

      case Moose::COORD_RSPHERICAL:
        if (_q_point[_qp](0) != 0.0)
          Ayy = Azz = _disp_r[_qp] / _q_point[_qp](0) + 1.0;
        break;
    }

    const auto detF = Axx * Ayy * Azz + Axy * Ayz * Azx + Axz * Ayx * Azy - Azx * Ayy * Axz -
                      Azy * Ayz * Axx - Azz * Ayx * Axy;
    density /= detF;
  }

  _density[_qp] = density;
}

template class DensityTempl<false>;
template class DensityTempl<true>;