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      158819 : DiracKernelTempl<T>::validParams()
      25             : {
      26      158819 :   InputParameters params = DiracKernelBase::validParams();
      27             :   if (std::is_same<T, Real>::value)
      28      288834 :     params.registerBase("DiracKernel");
      29             :   else if (std::is_same<T, RealVectorValue>::value)
      30       28804 :     params.registerBase("VectorDiracKernel");
      31             :   else
      32             :     ::mooseError("unsupported DiracKernelTempl specialization");
      33      158819 :   return params;
      34           0 : }
      35             : 
      36             : template <typename T>
      37         999 : 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        1998 :     _var(this->mooseVariableField()),
      46         999 :     _phi(_assembly.phi(_var)),
      47         999 :     _grad_phi(_assembly.gradPhi(_var)),
      48         999 :     _test(_var.phi()),
      49         999 :     _grad_test(_var.gradPhi()),
      50         999 :     _u(_var.sln()),
      51        3996 :     _grad_u(_var.gradSln())
      52             : {
      53         999 :   addMooseVariableDependency(&this->mooseVariableField());
      54             : 
      55             :   // Stateful material properties are not allowed on DiracKernels
      56         999 :   statefulPropertiesAllowed(false);
      57         999 : }
      58             : 
      59             : template <typename T>
      60             : void
      61      303444 : DiracKernelTempl<T>::computeResidual()
      62             : {
      63      303444 :   prepareVectorTag(_assembly, _var.number());
      64             : 
      65      303444 :   const std::vector<unsigned int> * multiplicities =
      66      303444 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      67      303444 :   unsigned int local_qp = 0;
      68      303444 :   Real multiplicity = 1.0;
      69             : 
      70      615364 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      71             :   {
      72      311920 :     _current_point = _physical_point[_qp];
      73      311920 :     if (isActiveAtPoint(_current_elem, _current_point))
      74             :     {
      75      311920 :       if (!_drop_duplicate_points)
      76         237 :         multiplicity = (*multiplicities)[local_qp++];
      77             : 
      78     1697641 :       for (_i = 0; _i < _test.size(); _i++)
      79     1385721 :         _local_re(_i) += multiplicity * computeQpResidual();
      80             :     }
      81             :   }
      82             : 
      83      303444 :   accumulateTaggedLocalResidual();
      84      303444 : }
      85             : 
      86             : template <typename T>
      87             : void
      88       42235 : DiracKernelTempl<T>::computeJacobian()
      89             : {
      90       42235 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      91             : 
      92       42235 :   const std::vector<unsigned int> * multiplicities =
      93       42235 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      94       42235 :   unsigned int local_qp = 0;
      95       42235 :   Real multiplicity = 1.0;
      96             : 
      97       85774 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      98             :   {
      99       43539 :     _current_point = _physical_point[_qp];
     100       43539 :     if (isActiveAtPoint(_current_elem, _current_point))
     101             :     {
     102       43539 :       if (!_drop_duplicate_points)
     103          20 :         multiplicity = (*multiplicities)[local_qp++];
     104             : 
     105      224873 :       for (_i = 0; _i < _test.size(); _i++)
     106     1024624 :         for (_j = 0; _j < _phi.size(); _j++)
     107      843290 :           _local_ke(_i, _j) += multiplicity * computeQpJacobian();
     108             :     }
     109             :   }
     110             : 
     111       42235 :   accumulateTaggedLocalMatrix();
     112       42235 : }
     113             : 
     114             : template <typename T>
     115             : void
     116       42419 : DiracKernelTempl<T>::computeOffDiagJacobian(const unsigned int jvar_num)
     117             : {
     118       42419 :   if (jvar_num == _var.number())
     119             :   {
     120       42235 :     computeJacobian();
     121             :   }
     122             :   else
     123             :   {
     124         184 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     125             : 
     126         184 :     const std::vector<unsigned int> * multiplicities =
     127         184 :         _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
     128         184 :     unsigned int local_qp = 0;
     129         184 :     Real multiplicity = 1.0;
     130             : 
     131         368 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     132             :     {
     133         184 :       _current_point = _physical_point[_qp];
     134         184 :       if (isActiveAtPoint(_current_elem, _current_point))
     135             :       {
     136         184 :         if (!_drop_duplicate_points)
     137           0 :           multiplicity = (*multiplicities)[local_qp++];
     138             : 
     139         920 :         for (_i = 0; _i < _test.size(); _i++)
     140        3680 :           for (_j = 0; _j < _phi.size(); _j++)
     141        2944 :             _local_ke(_i, _j) += multiplicity * computeQpOffDiagJacobian(jvar_num);
     142             :       }
     143             :     }
     144             : 
     145         184 :     accumulateTaggedLocalMatrix();
     146             :   }
     147       42419 : }
     148             : 
     149             : template <typename T>
     150             : Real
     151      840346 : DiracKernelTempl<T>::computeQpJacobian()
     152             : {
     153      840346 :   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>;