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             : #include "ElemElemConstraint.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "ElementPairInfo.h"
      15             : #include "FEProblem.h"
      16             : #include "MooseMesh.h"
      17             : #include "MooseVariableFE.h"
      18             : #include "SystemBase.h"
      19             : 
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : InputParameters
      23           0 : ElemElemConstraint::validParams()
      24             : {
      25           0 :   InputParameters params = Constraint::validParams();
      26           0 :   params.addParam<unsigned int>("interface_id", 0, "The id of the interface.");
      27           0 :   return params;
      28           0 : }
      29             : 
      30           0 : ElemElemConstraint::ElemElemConstraint(const InputParameters & parameters)
      31             :   : Constraint(parameters),
      32             :     NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, false, false),
      33             :     NeighborMooseVariableInterface<Real>(
      34             :         this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
      35           0 :     _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      36           0 :     _dim(_mesh.dimension()),
      37           0 :     _interface_id(getParam<unsigned int>("interface_id")),
      38           0 :     _var(_sys.getFieldVariable<Real>(_tid, parameters.get<NonlinearVariableName>("variable"))),
      39             : 
      40           0 :     _current_elem(_assembly.elem()),
      41             : 
      42           0 :     _neighbor_elem(_assembly.neighbor()),
      43             : 
      44           0 :     _u(_var.sln()),
      45           0 :     _grad_u(_var.gradSln()),
      46             : 
      47           0 :     _phi(_assembly.phi(_var)),
      48           0 :     _grad_phi(_assembly.gradPhi(_var)),
      49             : 
      50           0 :     _test(_var.phi()),
      51           0 :     _grad_test(_var.gradPhi()),
      52             : 
      53           0 :     _phi_neighbor(_assembly.phiNeighbor(_var)),
      54           0 :     _grad_phi_neighbor(_assembly.gradPhiNeighbor(_var)),
      55             : 
      56           0 :     _test_neighbor(_var.phiNeighbor()),
      57           0 :     _grad_test_neighbor(_var.gradPhiNeighbor()),
      58             : 
      59           0 :     _u_neighbor(_var.slnNeighbor()),
      60           0 :     _grad_u_neighbor(_var.gradSlnNeighbor())
      61             : {
      62           0 :   addMooseVariableDependency(&_var);
      63           0 : }
      64             : 
      65             : void
      66           0 : ElemElemConstraint::reinit(const ElementPairInfo & element_pair_info)
      67             : {
      68           0 :   reinitConstraintQuadrature(element_pair_info);
      69           0 : }
      70             : 
      71             : void
      72           0 : ElemElemConstraint::reinitConstraintQuadrature(const ElementPairInfo & element_pair_info)
      73             : {
      74           0 :   _constraint_q_point.resize(element_pair_info._elem1_constraint_q_point.size());
      75           0 :   _constraint_weight.resize(element_pair_info._elem1_constraint_JxW.size());
      76           0 :   std::copy(element_pair_info._elem1_constraint_q_point.begin(),
      77             :             element_pair_info._elem1_constraint_q_point.end(),
      78             :             _constraint_q_point.begin());
      79           0 :   std::copy(element_pair_info._elem1_constraint_JxW.begin(),
      80             :             element_pair_info._elem1_constraint_JxW.end(),
      81             :             _constraint_weight.begin());
      82           0 : }
      83             : 
      84             : void
      85           0 : ElemElemConstraint::computeElemNeighResidual(Moose::DGResidualType type)
      86             : {
      87             :   bool is_elem;
      88           0 :   if (type == Moose::Element)
      89           0 :     is_elem = true;
      90             :   else
      91           0 :     is_elem = false;
      92             : 
      93           0 :   const VariableTestValue & test_space = is_elem ? _test : _test_neighbor;
      94           0 :   if (is_elem)
      95           0 :     prepareVectorTag(_assembly, _var.number());
      96             :   else
      97           0 :     prepareVectorTagNeighbor(_assembly, _var.number());
      98           0 :   for (_qp = 0; _qp < _constraint_q_point.size(); _qp++)
      99           0 :     for (_i = 0; _i < test_space.size(); _i++)
     100           0 :       _local_re(_i) += _constraint_weight[_qp] * computeQpResidual(type);
     101           0 :   accumulateTaggedLocalResidual();
     102           0 : }
     103             : 
     104             : void
     105           0 : ElemElemConstraint::computeResidual()
     106             : {
     107             :   // Compute the residual for this element
     108           0 :   computeElemNeighResidual(Moose::Element);
     109             : 
     110             :   // Compute the residual for the neighbor
     111           0 :   computeElemNeighResidual(Moose::Neighbor);
     112           0 : }
     113             : 
     114             : void
     115           0 : ElemElemConstraint::computeElemNeighJacobian(Moose::DGJacobianType type)
     116             : {
     117           0 :   const VariableTestValue & test_space =
     118           0 :       (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
     119           0 :   const VariableTestValue & loc_phi =
     120           0 :       (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
     121           0 :   prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), type);
     122             : 
     123           0 :   for (_qp = 0; _qp < _constraint_q_point.size(); _qp++)
     124           0 :     for (_i = 0; _i < test_space.size(); _i++)
     125           0 :       for (_j = 0; _j < loc_phi.size(); _j++)
     126           0 :         _local_ke(_i, _j) += _constraint_weight[_qp] * computeQpJacobian(type);
     127             : 
     128           0 :   accumulateTaggedLocalMatrix();
     129           0 : }
     130             : 
     131             : void
     132           0 : ElemElemConstraint::computeJacobian()
     133             : {
     134             :   // Compute element-element Jacobian
     135           0 :   computeElemNeighJacobian(Moose::ElementElement);
     136             : 
     137             :   // Compute element-neighbor Jacobian
     138           0 :   computeElemNeighJacobian(Moose::ElementNeighbor);
     139             : 
     140             :   // Compute neighbor-element Jacobian
     141           0 :   computeElemNeighJacobian(Moose::NeighborElement);
     142             : 
     143             :   // Compute neighbor-neighbor Jacobian
     144           0 :   computeElemNeighJacobian(Moose::NeighborNeighbor);
     145           0 : }