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             : #pragma once
      11             : 
      12             : #include "InitialConditionBase.h"
      13             : 
      14             : // libMesh
      15             : #include "libmesh/point.h"
      16             : #include "libmesh/vector_value.h"
      17             : #include "libmesh/elem.h"
      18             : 
      19             : // forward declarations
      20             : class FEProblemBase;
      21             : class Assembly;
      22             : 
      23             : /**
      24             :  * This is a template class that implements the workhorse `compute` and `computeNodal` methods. The
      25             :  * former method is used for setting block initial conditions. It first projects the initial
      26             :  * condition field to nodes, then to edges, then to faces, then to interior dofs. The latter
      27             :  * `computeNodal` method sets dof values for boundary restricted initial conditions
      28             :  */
      29             : template <typename T>
      30             : class InitialConditionTempl : public InitialConditionBase
      31             : {
      32             : public:
      33             :   typedef typename OutputTools<T>::OutputShape ValueType;
      34             :   typedef typename OutputTools<T>::OutputGradient GradientType;
      35             :   typedef typename OutputTools<T>::OutputShapeGradient GradientShapeType;
      36             :   typedef typename OutputTools<T>::OutputData DataType;
      37             :   typedef libMesh::FEGenericBase<ValueType> FEBaseType;
      38             : 
      39             :   /**
      40             :    * Constructor
      41             :    *
      42             :    * @param parameters The parameters object holding data for the class to use.
      43             :    */
      44             :   InitialConditionTempl(const InputParameters & parameters);
      45             : 
      46             :   virtual ~InitialConditionTempl();
      47             : 
      48             :   static InputParameters validParams();
      49             : 
      50     4792709 :   virtual MooseVariableFEBase & variable() override { return _var; }
      51             : 
      52             :   virtual void compute() override;
      53             : 
      54             :   virtual void computeNodal(const Point & p) override;
      55             : 
      56             :   /**
      57             :    * The value of the variable at a point.
      58             :    *
      59             :    * This must be overridden by derived classes.
      60             :    */
      61             :   virtual T value(const Point & p) = 0;
      62             : 
      63             :   /**
      64             :    * The gradient of the variable at a point.
      65             :    *
      66             :    * This is optional.  Note that if you are using C1 continuous elements you will
      67             :    * want to use an initial condition that defines this!
      68             :    */
      69        9216 :   virtual GradientType gradient(const Point & /*p*/) { return GradientType(); };
      70             : 
      71             :   T gradientComponent(GradientType grad, unsigned int i);
      72             : 
      73             :   /**
      74             :    * set the temporary solution vector for node projections of C0 variables
      75             :    */
      76             :   void setCZeroVertices();
      77             :   /**
      78             :    * set the temporary solution vector for node projections of Hermite variables
      79             :    */
      80             :   void setHermiteVertices();
      81             :   /**
      82             :    * set the temporary solution vector for node projections of non-Hermitian C1 variables
      83             :    */
      84             :   void setOtherCOneVertices();
      85             : 
      86             :   /**
      87             :    * Helps perform multiplication of GradientTypes: a normal dot product for vectors and a
      88             :    * contraction for tensors
      89             :    */
      90           0 :   Real dotHelper(const libMesh::RealGradient & op1, const libMesh::RealGradient & op2)
      91             :   {
      92           0 :     return op1 * op2;
      93             :   }
      94           0 :   Real dotHelper(const RealTensor & op1, const RealTensor & op2) { return op1.contract(op2); }
      95           0 :   RealEigenVector dotHelper(const RealVectorArrayValue & op1, const libMesh::RealGradient & op2)
      96             :   {
      97           0 :     RealEigenVector v = op1.col(0) * op2(0);
      98           0 :     for (const auto i : make_range(Moose::dim))
      99           0 :       v += op1.col(i) * op2(i);
     100           0 :     return v;
     101           0 :   }
     102             : 
     103             :   /**
     104             :    * Perform the cholesky solves for edge, side, and interior projections
     105             :    */
     106             :   void choleskySolve(bool is_volume);
     107             : 
     108             :   /**
     109             :    * Assemble a small local system for cholesky solve
     110             :    */
     111             :   void choleskyAssembly(bool is_volume);
     112             : 
     113             : protected:
     114             :   FEProblemBase & _fe_problem;
     115             :   THREAD_ID _tid;
     116             : 
     117             :   /// Time
     118             :   Real & _t;
     119             : 
     120             :   /// The variable that this initial condition is acting upon.
     121             :   MooseVariableField<T> & _var;
     122             :   MooseVariableFE<T> * _fe_var;
     123             : 
     124             :   /// the finite element/volume assembly object
     125             :   Assembly & _assembly;
     126             : 
     127             :   /// The coordinate system type for this problem, references the value in Assembly
     128             :   const Moose::CoordinateSystemType & _coord_sys;
     129             : 
     130             :   /// The current element we are on will retrieving values at specific points in the domain. Note
     131             :   /// that this _IS_ valid even for nodes shared among several elements.
     132             :   const Elem * const & _current_elem;
     133             : 
     134             :   /// the volume of the current element
     135             :   const Real & _current_elem_volume;
     136             : 
     137             :   /// The current node if the point we are evaluating at also happens to be a node.
     138             :   /// Otherwise the pointer will be NULL.
     139             :   const Node * _current_node;
     140             : 
     141             :   /// The current quadrature point, contains the "nth" node number when visiting nodes.
     142             :   unsigned int _qp;
     143             : 
     144             :   /// Matrix storage member
     145             :   DenseMatrix<Real> _Ke;
     146             :   /// Linear b vector
     147             :   DenseVector<DataType> _Fe;
     148             :   /// Linear solution vector
     149             :   DenseVector<DataType> _Ue;
     150             : 
     151             :   /// The finite element type for the IC variable
     152             :   const libMesh::FEType & _fe_type;
     153             :   /// The type of continuity, e.g. C0, C1
     154             :   libMesh::FEContinuity _cont;
     155             : 
     156             :   /// The global DOF indices
     157             :   std::vector<dof_id_type> _dof_indices;
     158             :   /// Side/edge DOF indices
     159             :   std::vector<unsigned int> _side_dofs;
     160             : 
     161             :   /// The number of quadrature points for a given element
     162             :   unsigned int _n_qp;
     163             :   /// The number of nodes for a given element
     164             :   unsigned int _n_nodes;
     165             : 
     166             :   /// Whether the degree of freedom is fixed (true/false)
     167             :   std::vector<char> _dof_is_fixed;
     168             :   /// Stores the ids of the free dofs
     169             :   std::vector<int> _free_dof;
     170             : 
     171             :   /// The number of free dofs
     172             :   dof_id_type _free_dofs;
     173             :   /// The current dof being operated on
     174             :   dof_id_type _current_dof;
     175             :   /// number of dofs per node per variable
     176             :   dof_id_type _nc;
     177             : 
     178             :   /// pointers to shape functions
     179             :   const std::vector<std::vector<ValueType>> * _phi;
     180             :   /// pointers to shape function gradients
     181             :   const std::vector<std::vector<GradientShapeType>> * _dphi;
     182             :   /// pointers to the Jacobian * quadrature weights for current element
     183             :   const std::vector<Real> * _JxW;
     184             :   /// pointers to the xyz coordinates of the quadrature points for the current element
     185             :   const std::vector<Point> * _xyz_values;
     186             : 
     187             :   /// node counter
     188             :   unsigned int _n;
     189             :   /// the mesh dimension
     190             :   unsigned int _dim;
     191             : };
     192             : 
     193             : template <typename T>
     194             : InputParameters
     195      327040 : InitialConditionTempl<T>::validParams()
     196             : {
     197      327040 :   return InitialConditionBase::validParams();
     198             : }
     199             : 
     200             : typedef InitialConditionTempl<Real> InitialCondition;
     201             : typedef InitialConditionTempl<RealVectorValue> VectorInitialCondition;
     202             : typedef InitialConditionTempl<RealEigenVector> ArrayInitialCondition;
     203             : 
     204             : // Declare all the specializations, as the template specialization declaration below must know
     205             : template <>
     206             : void InitialConditionTempl<RealVectorValue>::setCZeroVertices();
     207             : template <>
     208             : RealVectorValue InitialConditionTempl<RealVectorValue>::gradientComponent(GradientType grad,
     209             :                                                                           unsigned int i);
     210             : template <>
     211             : RealEigenVector InitialConditionTempl<RealEigenVector>::gradientComponent(GradientType grad,
     212             :                                                                           unsigned int i);
     213             : template <>
     214             : void InitialConditionTempl<RealVectorValue>::setHermiteVertices();
     215             : template <>
     216             : void InitialConditionTempl<RealVectorValue>::setOtherCOneVertices();
     217             : template <>
     218             : void InitialConditionTempl<RealEigenVector>::choleskySolve(bool is_volume);
     219             : 
     220             : // Prevent implicit instantiation in other translation units where these classes are used
     221             : extern template class InitialConditionTempl<Real>;
     222             : extern template class InitialConditionTempl<RealVectorValue>;
     223             : extern template class InitialConditionTempl<RealEigenVector>;