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      158713 : DiracKernelTempl<T>::validParams()
      25             : {
      26      158713 :   InputParameters params = DiracKernelBase::validParams();
      27             :   if (std::is_same<T, Real>::value)
      28      144323 :     params.registerBase("DiracKernel");
      29             :   else if (std::is_same<T, RealVectorValue>::value)
      30       14390 :     params.registerBase("VectorDiracKernel");
      31             :   else
      32             :     ::mooseError("unsupported DiracKernelTempl specialization");
      33      158713 :   return params;
      34           0 : }
      35             : 
      36             : template <typename T>
      37         950 : 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        1900 :     _var(this->mooseVariableField()),
      46         950 :     _phi(_assembly.phi(_var)),
      47         950 :     _grad_phi(_assembly.gradPhi(_var)),
      48         950 :     _test(_var.phi()),
      49         950 :     _grad_test(_var.gradPhi()),
      50         950 :     _u(_var.sln()),
      51         950 :     _grad_u(_var.gradSln()),
      52         950 :     _drop_duplicate_points(parameters.get<bool>("drop_duplicate_points")),
      53         950 :     _point_not_found_behavior(
      54        1900 :         parameters.get<MooseEnum>("point_not_found_behavior").getEnum<PointNotFoundBehavior>())
      55             : {
      56         950 :   addMooseVariableDependency(&this->mooseVariableField());
      57             : 
      58             :   // Stateful material properties are not allowed on DiracKernels
      59         950 :   statefulPropertiesAllowed(false);
      60         950 : }
      61             : 
      62             : template <typename T>
      63             : void
      64      269579 : DiracKernelTempl<T>::computeResidual()
      65             : {
      66      269579 :   prepareVectorTag(_assembly, _var.number());
      67             : 
      68      269579 :   const std::vector<unsigned int> * multiplicities =
      69      269579 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      70      269579 :   unsigned int local_qp = 0;
      71      269579 :   Real multiplicity = 1.0;
      72             : 
      73      546646 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      74             :   {
      75      277067 :     _current_point = _physical_point[_qp];
      76      277067 :     if (isActiveAtPoint(_current_elem, _current_point))
      77             :     {
      78      277067 :       if (!_drop_duplicate_points)
      79         193 :         multiplicity = (*multiplicities)[local_qp++];
      80             : 
      81     1520570 :       for (_i = 0; _i < _test.size(); _i++)
      82     1243503 :         _local_re(_i) += multiplicity * computeQpResidual();
      83             :     }
      84             :   }
      85             : 
      86      269579 :   accumulateTaggedLocalResidual();
      87      269579 : }
      88             : 
      89             : template <typename T>
      90             : void
      91       37375 : DiracKernelTempl<T>::computeJacobian()
      92             : {
      93       37375 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      94             : 
      95       37375 :   const std::vector<unsigned int> * multiplicities =
      96       37375 :       _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
      97       37375 :   unsigned int local_qp = 0;
      98       37375 :   Real multiplicity = 1.0;
      99             : 
     100       75902 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     101             :   {
     102       38527 :     _current_point = _physical_point[_qp];
     103       38527 :     if (isActiveAtPoint(_current_elem, _current_point))
     104             :     {
     105       38527 :       if (!_drop_duplicate_points)
     106          16 :         multiplicity = (*multiplicities)[local_qp++];
     107             : 
     108      199599 :       for (_i = 0; _i < _test.size(); _i++)
     109      919576 :         for (_j = 0; _j < _phi.size(); _j++)
     110      758504 :           _local_ke(_i, _j) += multiplicity * computeQpJacobian();
     111             :     }
     112             :   }
     113             : 
     114       37375 :   accumulateTaggedLocalMatrix();
     115       37375 : }
     116             : 
     117             : template <typename T>
     118             : void
     119       37535 : DiracKernelTempl<T>::computeOffDiagJacobian(const unsigned int jvar_num)
     120             : {
     121       37535 :   if (jvar_num == _var.number())
     122             :   {
     123       37375 :     computeJacobian();
     124             :   }
     125             :   else
     126             :   {
     127         160 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     128             : 
     129         160 :     const std::vector<unsigned int> * multiplicities =
     130         160 :         _drop_duplicate_points ? NULL : &_local_dirac_kernel_info.getPoints()[_current_elem].second;
     131         160 :     unsigned int local_qp = 0;
     132         160 :     Real multiplicity = 1.0;
     133             : 
     134         320 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     135             :     {
     136         160 :       _current_point = _physical_point[_qp];
     137         160 :       if (isActiveAtPoint(_current_elem, _current_point))
     138             :       {
     139         160 :         if (!_drop_duplicate_points)
     140           0 :           multiplicity = (*multiplicities)[local_qp++];
     141             : 
     142         800 :         for (_i = 0; _i < _test.size(); _i++)
     143        3200 :           for (_j = 0; _j < _phi.size(); _j++)
     144        2560 :             _local_ke(_i, _j) += multiplicity * computeQpOffDiagJacobian(jvar_num);
     145             :       }
     146             :     }
     147             : 
     148         160 :     accumulateTaggedLocalMatrix();
     149             :   }
     150       37535 : }
     151             : 
     152             : template <typename T>
     153             : Real
     154      755944 : DiracKernelTempl<T>::computeQpJacobian()
     155             : {
     156      755944 :   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>;