MooseVariableData.h

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

#pragma once

#include "MooseArray.h"
#include "MooseTypes.h"
#include "MooseVariableDataBase.h"
#include "Conversion.h"

#include "libmesh/tensor_tools.h"
#include "libmesh/vector_value.h"
#include "libmesh/tensor_value.h"
#include "libmesh/type_n_tensor.h"
#include "libmesh/enum_fe_family.h"
#include "libmesh/fe_type.h"
#include "ADUtils.h"

#include <functional>
#include <vector>

namespace libMesh
{
class QBase;
class DofMap;
template <typename>
class DenseVector;
}

using libMesh::DenseVector;
using libMesh::DofMap;
using libMesh::QBase;
using libMesh::VectorValue;

class TimeIntegrator;
class Assembly;
class SubProblem;
template <typename>
class MooseVariableFE;
class SystemBase;

template <typename OutputType>
class MooseVariableData : public MooseVariableDataBase<OutputType>
{
public:
  // type for gradient, second and divergence of template class OutputType
  typedef typename libMesh::TensorTools::IncrementRank<OutputType>::type OutputGradient;
  typedef typename libMesh::TensorTools::IncrementRank<OutputGradient>::type OutputSecond;
  typedef typename libMesh::TensorTools::DecrementRank<OutputType>::type OutputDivergence;

  // shortcut for types storing values on quadrature points
  typedef MooseArray<OutputType> FieldVariableValue;
  typedef MooseArray<OutputGradient> FieldVariableGradient;
  typedef MooseArray<OutputSecond> FieldVariableSecond;
  typedef MooseArray<OutputType> FieldVariableCurl;
  typedef MooseArray<OutputDivergence> FieldVariableDivergence;

  // shape function type for the template class OutputType
  typedef typename Moose::ShapeType<OutputType>::type OutputShape;

  // type for gradient, second and divergence of shape functions of template class OutputType
  typedef typename libMesh::TensorTools::IncrementRank<OutputShape>::type OutputShapeGradient;
  typedef typename libMesh::TensorTools::IncrementRank<OutputShapeGradient>::type OutputShapeSecond;
  typedef typename libMesh::TensorTools::DecrementRank<OutputShape>::type OutputShapeDivergence;

  // shortcut for types storing shape function values on quadrature points
  typedef MooseArray<std::vector<OutputShape>> FieldVariablePhiValue;
  typedef MooseArray<std::vector<OutputShapeGradient>> FieldVariablePhiGradient;
  typedef MooseArray<std::vector<OutputShapeSecond>> FieldVariablePhiSecond;
  typedef MooseArray<std::vector<OutputShape>> FieldVariablePhiCurl;
  typedef MooseArray<std::vector<OutputShapeDivergence>> FieldVariablePhiDivergence;

  // shortcut for types storing test function values on quadrature points
  // Note: here we assume the types are the same as of shape functions.
  typedef MooseArray<std::vector<OutputShape>> FieldVariableTestValue;
  typedef MooseArray<std::vector<OutputShapeGradient>> FieldVariableTestGradient;
  typedef MooseArray<std::vector<OutputShapeSecond>> FieldVariableTestSecond;
  typedef MooseArray<std::vector<OutputShape>> FieldVariableTestCurl;
  typedef MooseArray<std::vector<OutputShapeDivergence>> FieldVariableTestDivergence;

  // DoF value type for the template class OutputType
  typedef typename Moose::DOFType<OutputType>::type OutputData;
  typedef MooseArray<OutputData> DoFValue;

  MooseVariableData(const MooseVariableField<OutputType> & var,
                    SystemBase & sys,
                    THREAD_ID tid,
                    Moose::ElementType element_type,
                    const QBase * const & qrule_in,
                    const QBase * const & qrule_face_in,
                    const Node * const & node,
                    const Elem * const & elem);

  bool needsAD() const { return _need_ad; }

  void setGeometry(Moose::GeometryType gm_type);


  void computeValues();

  void computeMonomialValues();

  void computeAD(const unsigned int num_dofs, const unsigned int nqp);

  void computeNodalValues();


  const FieldVariablePhiValue & phi() const { return *_phi; }

  const FieldVariablePhiValue & phiFace() const { return *_phi_face; }

  const FieldVariablePhiGradient & gradPhi() const { return *_grad_phi; }

  const MappedArrayVariablePhiGradient & arrayGradPhi() const
  {
    mooseAssert(var().fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY, "Not an array variable");
    return _mapped_grad_phi;
  }

  const FieldVariablePhiGradient & gradPhiFace() const { return *_grad_phi_face; }

  const MappedArrayVariablePhiGradient & arrayGradPhiFace() const
  {
    mooseAssert(var().fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY, "Not an array variable");
    return _mapped_grad_phi_face;
  }

  const FieldVariablePhiSecond & secondPhi() const;

  const FieldVariablePhiSecond & secondPhiFace() const;

  const FieldVariablePhiCurl & curlPhi() const;

  const FieldVariablePhiCurl & curlPhiFace() const;

  const FieldVariablePhiDivergence & divPhi() const;

  const FieldVariablePhiDivergence & divPhiFace() const;

  const ADTemplateVariablePhiGradient<OutputShape> & adGradPhi() const;

  const ADTemplateVariablePhiGradient<OutputShape> & adGradPhiFace() const;

  std::size_t phiSize() const { return _phi->size(); }

  std::size_t phiFaceSize() const { return _phi_face->size(); }

  bool usesSecondPhi() const
  {
    return _need_second || _need_second_old || _need_second_older || _need_second_previous_nl;
  }

  bool computingCurl() const { return _need_curl || _need_curl_old; }

  bool computingDiv() const { return _need_div || _need_div_old; }


  bool isNodal() const override { return _is_nodal; }
  bool hasDoFsOnNodes() const override { return _continuity != libMesh::DISCONTINUOUS; }
  libMesh::FEContinuity getContinuity() const override { return _continuity; };
  const Node * const & node() const { return _node; }
  const dof_id_type & nodalDofIndex() const { return _nodal_dof_index; }
  bool isNodalDefined() const { return _has_dof_indices; }

  const Elem * const & currentElem() const { return _elem; }

  const unsigned int & currentSide() const { return _current_side; }

  void prepareIC();


  const FieldVariableGradient & gradSlnDot() const;

  const FieldVariableGradient & gradSlnDotDot() const;

  const FieldVariableSecond & secondSln(Moose::SolutionState state) const;

  const FieldVariableCurl & curlSln(Moose::SolutionState state) const;

  const FieldVariableDivergence & divSln(Moose::SolutionState state) const;

  const ADTemplateVariableValue<OutputType> & adSln() const
  {
    _need_ad = _need_ad_u = true;
    return _ad_u;
  }

  const ADTemplateVariableGradient<OutputType> & adGradSln() const
  {
    _need_ad = _need_ad_grad_u = true;
    return _ad_grad_u;
  }

  const ADTemplateVariableGradient<OutputType> & adGradSlnDot() const
  {
    _need_ad = _need_ad_grad_u_dot = true;

    if (!_time_integrator)
      // If we don't have a time integrator (this will be the case for variables that are a part of
      // the AuxiliarySystem) then we have no way to calculate _ad_grad_u_dot and we are just going
      // to copy the values from _grad_u_dot. Of course in order to be able to do that we need to
      // calculate _grad_u_dot
      _need_grad_dot = true;

    return _ad_grad_u_dot;
  }

  const ADTemplateVariableSecond<OutputType> & adSecondSln() const
  {
    _need_ad = _need_ad_second_u = true;
    secondPhi();
    secondPhiFace();
    return _ad_second_u;
  }

  const ADTemplateVariableValue<OutputType> & adUDot() const;

  const ADTemplateVariableValue<OutputType> & adUDotDot() const;

  const FieldVariableValue & uDot() const;

  const FieldVariableValue & uDotDot() const;

  const FieldVariableValue & uDotOld() const;

  const FieldVariableValue & uDotDotOld() const;

  const VariableValue & duDotDu() const
  {
    _need_du_dot_du = true;
    return _du_dot_du;
  }

  const VariableValue & duDotDotDu() const
  {
    _need_du_dotdot_du = true;
    return _du_dotdot_du;
  }

  const ADTemplateVariableCurl<OutputType> & adCurlSln() const
  {
    _need_ad = _need_ad_curl_u = true;
    curlPhi();
    curlPhiFace();
    return _ad_curl_u;
  }


  const OutputType & nodalValueDot() const;
  const OutputType & nodalValueDotDot() const;
  const OutputType & nodalValueDotOld() const;
  const OutputType & nodalValueDotDotOld() const;
  const OutputType & nodalValueDuDotDu() const;
  const OutputType & nodalValueDuDotDotDu() const;

  const typename Moose::ADType<OutputType>::type & adNodalValue() const;

  void setDofValues(const DenseVector<OutputData> & values);


  void setDofValue(const OutputData & value, unsigned int index);

  void insertNodalValue(libMesh::NumericVector<libMesh::Number> & residual, const OutputData & v);
  OutputData getNodalValue(const Node & node, Moose::SolutionState state) const;
  OutputData
  getElementalValue(const Elem * elem, Moose::SolutionState state, unsigned int idx = 0) const;


  void getDofIndices(const Elem * elem, std::vector<dof_id_type> & dof_indices) const;
  const std::vector<dof_id_type> & dofIndices() const { return _dof_indices; }
  unsigned int numberOfDofs() const { return _dof_indices.size(); }
  void clearDofIndices() { _dof_indices.clear(); }

  void prepare();

  void reinitNode();

  void reinitAux();

  void reinitNodes(const std::vector<dof_id_type> & nodes);

  void addSolution(libMesh::NumericVector<libMesh::Number> & sol,
                   const DenseVector<libMesh::Number> & v) const;


  const DoFValue & dofValuesDot() const;
  const DoFValue & dofValuesDotOld() const;
  const DoFValue & dofValuesDotDot() const;
  const DoFValue & dofValuesDotDotOld() const;
  const MooseArray<libMesh::Number> & dofValuesDuDotDu() const;
  const MooseArray<libMesh::Number> & dofValuesDuDotDotDu() const;

  const MooseArray<ADReal> & adDofValues() const;

  const MooseArray<ADReal> & adDofValuesDot() const;


  const FieldVariableValue & increment() const { return _increment; }

  void computeIncrementAtQps(const libMesh::NumericVector<libMesh::Number> & increment_vec);

  void computeIncrementAtNode(const libMesh::NumericVector<libMesh::Number> & increment_vec);

private:
  void assignADNodalValue(const ADReal & value, const unsigned int & component);
  void fetchADNodalValues();

  template <bool monomial>
  void computeValuesInternal();

  const libMesh::FEType & _fe_type;

  const unsigned int _var_num;

  const Assembly & _assembly;

  Moose::ElementType _element_type;

  bool _is_nodal;

  dof_id_type _nodal_dof_index;

  libMesh::FEContinuity _continuity;

  FieldVariableValue _increment;

  typename Moose::ADType<OutputType>::type _ad_nodal_value;

  ADReal _ad_zero;

  mutable bool _need_ad_u_dot;
  mutable bool _need_ad_u_dotdot;

  mutable bool _need_second;
  mutable bool _need_second_old;
  mutable bool _need_second_older;
  mutable bool _need_second_previous_nl;

  mutable bool _need_curl;
  mutable bool _need_curl_old;
  mutable bool _need_curl_older;

  mutable bool _need_div;
  mutable bool _need_div_old;
  mutable bool _need_div_older;

  mutable bool _need_ad;
  mutable bool _need_ad_u;
  mutable bool _need_ad_grad_u;
  mutable bool _need_ad_grad_u_dot;
  mutable bool _need_ad_second_u;
  mutable bool _need_ad_curl_u;

  bool _has_dof_indices;

  FieldVariableGradient _grad_u_dot;
  FieldVariableGradient _grad_u_dotdot;

  FieldVariableSecond _second_u;
  FieldVariableSecond _second_u_old;
  FieldVariableSecond _second_u_older;
  FieldVariableSecond _second_u_previous_nl;

  FieldVariableCurl _curl_u;
  FieldVariableCurl _curl_u_old;
  FieldVariableCurl _curl_u_older;

  FieldVariableDivergence _div_u;
  FieldVariableDivergence _div_u_old;
  FieldVariableDivergence _div_u_older;

  ADTemplateVariableValue<OutputType> _ad_u;
  ADTemplateVariableGradient<OutputType> _ad_grad_u;
  ADTemplateVariableSecond<OutputType> _ad_second_u;
  MooseArray<ADReal> _ad_dof_values;
  MooseArray<ADReal> _ad_dofs_dot;
  MooseArray<ADReal> _ad_dofs_dotdot;
  ADTemplateVariableValue<OutputType> _ad_u_dot;
  ADTemplateVariableValue<OutputType> _ad_u_dotdot;
  ADTemplateVariableGradient<OutputType> _ad_grad_u_dot;
  ADTemplateVariableCurl<OutputType> _ad_curl_u;

  // time derivatives

  FieldVariableValue _u_dot;

  FieldVariableValue _u_dotdot;

  FieldVariableValue _u_dot_old;

  FieldVariableValue _u_dotdot_old;

  VariableValue _du_dot_du;

  VariableValue _du_dotdot_du;

  const QBase * const & _qrule;
  const QBase * const & _qrule_face;

  // Shape function values, gradients, second derivatives
  const FieldVariablePhiValue * _phi;
  const FieldVariablePhiGradient * _grad_phi;
  mutable const FieldVariablePhiSecond * _second_phi;
  mutable const FieldVariablePhiCurl * _curl_phi;
  mutable const FieldVariablePhiDivergence * _div_phi;

  // Mapped array phi
  MappedArrayVariablePhiGradient _mapped_grad_phi;
  MappedArrayVariablePhiGradient _mapped_grad_phi_face;
  MappedArrayVariablePhiGradient _mapped_grad_phi_neighbor;
  MappedArrayVariablePhiGradient _mapped_grad_phi_face_neighbor;

  // Values, gradients and second derivatives of shape function on faces
  const FieldVariablePhiValue * _phi_face;
  const FieldVariablePhiGradient * _grad_phi_face;
  mutable const FieldVariablePhiSecond * _second_phi_face;
  mutable const FieldVariablePhiCurl * _curl_phi_face;
  mutable const FieldVariablePhiDivergence * _div_phi_face;

  const ADTemplateVariablePhiGradient<OutputShape> * _ad_grad_phi;
  const ADTemplateVariablePhiGradient<OutputShape> * _ad_grad_phi_face;

  const QBase * _current_qrule;
  const FieldVariablePhiValue * _current_phi;
  const FieldVariablePhiGradient * _current_grad_phi;
  const FieldVariablePhiSecond * _current_second_phi;
  const FieldVariablePhiCurl * _current_curl_phi;
  const FieldVariablePhiDivergence * _current_div_phi;
  const ADTemplateVariablePhiGradient<OutputShape> * _current_ad_grad_phi;

  // dual mortar
  const bool _use_dual;

  std::function<const typename OutputTools<OutputType>::VariablePhiValue &(const Assembly &,
                                                                           libMesh::FEType)>
      _phi_assembly_method;
  std::function<const typename OutputTools<OutputShape>::VariablePhiValue &(const Assembly &,
                                                                            libMesh::FEType)>
      _phi_face_assembly_method;

  std::function<const typename OutputTools<OutputShape>::VariablePhiGradient &(const Assembly &,
                                                                               libMesh::FEType)>
      _grad_phi_assembly_method;
  std::function<const typename OutputTools<OutputShape>::VariablePhiGradient &(const Assembly &,
                                                                               libMesh::FEType)>
      _grad_phi_face_assembly_method;

  std::function<const typename OutputTools<OutputShape>::VariablePhiSecond &(const Assembly &,
                                                                             libMesh::FEType)>
      _second_phi_assembly_method;
  std::function<const typename OutputTools<OutputShape>::VariablePhiSecond &(const Assembly &,
                                                                             libMesh::FEType)>
      _second_phi_face_assembly_method;

  std::function<const typename OutputTools<OutputShape>::VariablePhiCurl &(const Assembly &,
                                                                           libMesh::FEType)>
      _curl_phi_assembly_method;
  std::function<const typename OutputTools<OutputShape>::VariablePhiCurl &(const Assembly &,
                                                                           libMesh::FEType)>
      _curl_phi_face_assembly_method;

  std::function<const typename OutputTools<OutputShape>::VariablePhiDivergence &(const Assembly &,
                                                                                 libMesh::FEType)>
      _div_phi_assembly_method;

  std::function<const typename OutputTools<OutputShape>::VariablePhiDivergence &(const Assembly &,
                                                                                 libMesh::FEType)>
      _div_phi_face_assembly_method;

  std::function<const ADTemplateVariablePhiGradient<OutputShape> &(const Assembly &,
                                                                   libMesh::FEType)>
      _ad_grad_phi_assembly_method;
  std::function<const ADTemplateVariablePhiGradient<OutputShape> &(const Assembly &,
                                                                   libMesh::FEType)>
      _ad_grad_phi_face_assembly_method;

  const TimeIntegrator * _time_integrator;

  const Node * const & _node;

  const Elem * const & _elem;

  const bool _displaced;

  const unsigned int & _current_side;

  ADReal _ad_real_dummy = 0;

  using MooseVariableDataBase<OutputType>::var;
  using MooseVariableDataBase<OutputType>::_sys;
  using MooseVariableDataBase<OutputType>::_subproblem;
  using MooseVariableDataBase<OutputType>::_need_vector_tag_dof_u;
  using MooseVariableDataBase<OutputType>::_need_matrix_tag_dof_u;
  using MooseVariableDataBase<OutputType>::_vector_tags_dof_u;
  using MooseVariableDataBase<OutputType>::_matrix_tags_dof_u;
  using MooseVariableDataBase<OutputType>::_vector_tag_u;
  using MooseVariableDataBase<OutputType>::_need_vector_tag_u;
  using MooseVariableDataBase<OutputType>::_vector_tag_grad;
  using MooseVariableDataBase<OutputType>::_need_vector_tag_grad;
  using MooseVariableDataBase<OutputType>::_matrix_tag_u;
  using MooseVariableDataBase<OutputType>::_need_matrix_tag_u;
  using MooseVariableDataBase<OutputType>::_dof_indices;
  using MooseVariableDataBase<OutputType>::_has_dof_values;
  using MooseVariableDataBase<OutputType>::fetchDoFValues;
  using MooseVariableDataBase<OutputType>::assignNodalValue;
  using MooseVariableDataBase<OutputType>::zeroSizeDofValues;
  using MooseVariableDataBase<OutputType>::_solution_tag;
  using MooseVariableDataBase<OutputType>::_old_solution_tag;
  using MooseVariableDataBase<OutputType>::_older_solution_tag;
  using MooseVariableDataBase<OutputType>::_previous_nl_solution_tag;
  using MooseVariableDataBase<OutputType>::_dof_map;
  using MooseVariableDataBase<OutputType>::_need_u_dot;
  using MooseVariableDataBase<OutputType>::_need_u_dotdot;
  using MooseVariableDataBase<OutputType>::_need_u_dot_old;
  using MooseVariableDataBase<OutputType>::_need_u_dotdot_old;
  using MooseVariableDataBase<OutputType>::_need_du_dot_du;
  using MooseVariableDataBase<OutputType>::_need_du_dotdot_du;
  using MooseVariableDataBase<OutputType>::_need_grad_dot;
  using MooseVariableDataBase<OutputType>::_need_grad_dotdot;
  using MooseVariableDataBase<OutputType>::_need_dof_values_dot;
  using MooseVariableDataBase<OutputType>::_need_dof_values_dotdot;
  using MooseVariableDataBase<OutputType>::_need_dof_values_dot_old;
  using MooseVariableDataBase<OutputType>::_need_dof_values_dotdot_old;
  using MooseVariableDataBase<OutputType>::_need_dof_du_dot_du;
  using MooseVariableDataBase<OutputType>::_need_dof_du_dotdot_du;
  using MooseVariableDataBase<OutputType>::_dof_values_dot;
  using MooseVariableDataBase<OutputType>::_dof_values_dotdot;
  using MooseVariableDataBase<OutputType>::_dof_values_dot_old;
  using MooseVariableDataBase<OutputType>::_dof_values_dotdot_old;
  using MooseVariableDataBase<OutputType>::_dof_du_dot_du;
  using MooseVariableDataBase<OutputType>::_dof_du_dotdot_du;
  using MooseVariableDataBase<OutputType>::_tid;
  using MooseVariableDataBase<OutputType>::_nodal_value_dot;
  using MooseVariableDataBase<OutputType>::_nodal_value_dotdot;
  using MooseVariableDataBase<OutputType>::_nodal_value_dot_old;
  using MooseVariableDataBase<OutputType>::_nodal_value_dotdot_old;
  using MooseVariableDataBase<OutputType>::_required_vector_tags;
};


template <typename OutputType>
const MooseArray<ADReal> &
MooseVariableData<OutputType>::adDofValues() const
{
  _need_ad = true;
  return _ad_dof_values;
}

template <typename OutputType>
const MooseArray<ADReal> &
MooseVariableData<OutputType>::adDofValuesDot() const
{
  _need_ad = _need_ad_u_dot = true;
  if (!_time_integrator)
    // See explanation in adUDot() body
    _need_u_dot = true;
  return _ad_dofs_dot;
}

template <typename OutputType>
const typename Moose::ADType<OutputType>::type &
MooseVariableData<OutputType>::adNodalValue() const
{
  _need_ad = true;
  return _ad_nodal_value;
}

template <typename OutputType>
const ADTemplateVariableValue<OutputType> &
MooseVariableData<OutputType>::adUDot() const
{
  _need_ad = _need_ad_u_dot = true;

  if (!_time_integrator)
    // If we don't have a time integrator (this will be the case for variables that are a part of
    // the AuxiliarySystem) then we have no way to calculate _ad_u_dot and we are just going to
    // copy the values from _u_dot. Of course in order to be able to do that we need to calculate
    // _u_dot
    _need_u_dot = true;

  return _ad_u_dot;
}

template <typename OutputType>
const ADTemplateVariableValue<OutputType> &
MooseVariableData<OutputType>::adUDotDot() const
{
  _need_ad = _need_ad_u_dotdot = true;

  if (!_time_integrator)
    // If we don't have a time integrator (this will be the case for variables that are a part
    // of the AuxiliarySystem) then we have no way to calculate _ad_u_dotdot and we are just
    // going to copy the values from _u_dotdot. Of course in order to be able to do that we need
    // to calculate _u_dotdot
    _need_u_dotdot = true;

  return _ad_u_dotdot;
}

template <typename OutputType>
const ADTemplateVariablePhiGradient<typename MooseVariableData<OutputType>::OutputShape> &
MooseVariableData<OutputType>::adGradPhi() const
{
  if (_element_type == Moose::ElementType::Neighbor || _element_type == Moose::ElementType::Lower)
    mooseError("Unsupported element type: ", Moose::stringify(_element_type));
  mooseAssert(_ad_grad_phi, "this should be non-null");
  return *_ad_grad_phi;
}

template <typename OutputType>
const ADTemplateVariablePhiGradient<typename MooseVariableData<OutputType>::OutputShape> &
MooseVariableData<OutputType>::adGradPhiFace() const
{
  if (_element_type == Moose::ElementType::Neighbor || _element_type == Moose::ElementType::Lower)
    mooseError("Unsupported element type: ", Moose::stringify(_element_type));
  mooseAssert(_ad_grad_phi_face, "this should be non-null");
  return *_ad_grad_phi_face;
}