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             : // Moose includes
      11             : #include "DiracKernel.h"
      12             : #include "Assembly.h"
      13             : #include "DiracKernelBase.h"
      14             : #include "MooseError.h"
      15             : #include "SystemBase.h"
      16             : #include "Problem.h"
      17             : #include "MooseMesh.h"
      18             : 
      19             : #include "libmesh/libmesh_common.h"
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : template <typename T>
      23             : InputParameters
      24       35418 : DiracKernelTempl<T>::validParams()
      25             : {
      26       35418 :   InputParameters params = DiracKernelBase::validParams();
      27             :   if (std::is_same<T, Real>::value)
      28       64464 :     params.registerBase("DiracKernel");
      29             :   else if (std::is_same<T, RealVectorValue>::value)
      30        6372 :     params.registerBase("VectorDiracKernel");
      31             :   else
      32             :     ::mooseError("unsupported DiracKernelTempl specialization");
      33       35418 :   return params;
      34           0 : }
      35             : 
      36             : template <typename T>
      37         931 : DiracKernelTempl<T>::DiracKernelTempl(const InputParameters & parameters)
      38             :   : DiracKernelBase(parameters),
      39             :     MooseVariableInterface<T>(this,
      40             :                               false,
      41             :                               "variable",
      42             :                               Moose::VarKindType::VAR_SOLVER,
      43             :                               std::is_same<T, Real>::value ? Moose::VarFieldType::VAR_FIELD_STANDARD
      44             :                                                            : Moose::VarFieldType::VAR_FIELD_VECTOR),
      45        1862 :     _var(this->mooseVariableField()),
      46         931 :     _phi(_assembly.phi(_var)),
      47         931 :     _grad_phi(_assembly.gradPhi(_var)),
      48         931 :     _test(_var.phi()),
      49         931 :     _grad_test(_var.gradPhi()),
      50         931 :     _u(_var.sln()),
      51        3724 :     _grad_u(_var.gradSln())
      52             : {
      53         931 :   addMooseVariableDependency(&this->mooseVariableField());
      54             : 
      55             :   // Stateful material properties are not allowed on DiracKernels
      56         931 :   statefulPropertiesAllowed(false);
      57         931 : }
      58             : 
      59             : template <typename T>
      60             : void
      61      269223 : DiracKernelTempl<T>::computeResidual()
      62             : {
      63      269223 :   prepareVectorTag(_assembly, _var.number());
      64             : 
      65      269223 :   const std::vector<Real> * point_values =
      66      269223 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      67      269223 :   unsigned int local_qp = 0;
      68      269223 :   Real point_value = 1.0;
      69             : 
      70      545986 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      71             :   {
      72      276763 :     _current_point = _physical_point[_qp];
      73      276763 :     if (isActiveAtPoint(_current_elem, _current_point))
      74             :     {
      75      276763 :       if (!_drop_duplicate_points)
      76        8243 :         point_value = (*point_values)[local_qp++];
      77             : 
      78     1492914 :       for (_i = 0; _i < _test.size(); _i++)
      79     1216151 :         _local_re(_i) += point_value * computeQpResidual();
      80             :     }
      81             :   }
      82             : 
      83      269223 :   accumulateTaggedLocalResidual();
      84      269223 : }
      85             : 
      86             : template <typename T>
      87             : void
      88       37534 : DiracKernelTempl<T>::computeJacobian()
      89             : {
      90       37534 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      91             : 
      92       37534 :   const std::vector<Real> * point_values =
      93       37534 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      94       37534 :   unsigned int local_qp = 0;
      95       37534 :   Real point_value = 1.0;
      96             : 
      97       76228 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      98             :   {
      99       38694 :     _current_point = _physical_point[_qp];
     100       38694 :     if (isActiveAtPoint(_current_elem, _current_point))
     101             :     {
     102       38694 :       if (!_drop_duplicate_points)
     103         462 :         point_value = (*point_values)[local_qp++];
     104             : 
     105      199102 :       for (_i = 0; _i < _test.size(); _i++)
     106      894944 :         for (_j = 0; _j < _phi.size(); _j++)
     107      734536 :           _local_ke(_i, _j) += point_value * computeQpJacobian();
     108             :     }
     109             :   }
     110             : 
     111       37534 :   accumulateTaggedLocalMatrix();
     112       37534 : }
     113             : 
     114             : template <typename T>
     115             : void
     116       37698 : DiracKernelTempl<T>::computeOffDiagJacobian(const unsigned int jvar_num)
     117             : {
     118       37698 :   if (jvar_num == _var.number())
     119             :   {
     120       37534 :     computeJacobian();
     121             :   }
     122             :   else
     123             :   {
     124         164 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     125             : 
     126         164 :     const std::vector<Real> * point_values =
     127         164 :         _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
     128         164 :     unsigned int local_qp = 0;
     129         164 :     Real point_value = 1.0;
     130             : 
     131         328 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     132             :     {
     133         164 :       _current_point = _physical_point[_qp];
     134         164 :       if (isActiveAtPoint(_current_elem, _current_point))
     135             :       {
     136         164 :         if (!_drop_duplicate_points)
     137           0 :           point_value = (*point_values)[local_qp++];
     138             : 
     139         820 :         for (_i = 0; _i < _test.size(); _i++)
     140        3280 :           for (_j = 0; _j < _phi.size(); _j++)
     141        2624 :             _local_ke(_i, _j) += point_value * computeQpOffDiagJacobian(jvar_num);
     142             :       }
     143             :     }
     144             : 
     145         164 :     accumulateTaggedLocalMatrix();
     146             :   }
     147       37698 : }
     148             : 
     149             : template <typename T>
     150             : Real
     151      731912 : DiracKernelTempl<T>::computeQpJacobian()
     152             : {
     153      731912 :   return 0;
     154             : }
     155             : 
     156             : template <typename T>
     157             : Real
     158           0 : DiracKernelTempl<T>::computeQpOffDiagJacobian(unsigned int /*jvar*/)
     159             : {
     160           0 :   return 0;
     161             : }
     162             : 
     163             : // Explicitly instantiates the two versions of the DiracKernelTempl class
     164             : template class DiracKernelTempl<Real>;
     165             : template class DiracKernelTempl<RealVectorValue>;