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