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 "XFEMSingleVariableConstraint.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "ElementPairInfo.h"
      15             : #include "FEProblem.h"
      16             : #include "GeometricCutUserObject.h"
      17             : #include "XFEM.h"
      18             : #include "Function.h"
      19             : 
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : registerMooseObject("XFEMApp", XFEMSingleVariableConstraint);
      23             : 
      24             : InputParameters
      25         384 : XFEMSingleVariableConstraint::validParams()
      26             : {
      27         384 :   InputParameters params = ElemElemConstraint::validParams();
      28         768 :   params.addParam<Real>("alpha",
      29         768 :                         100,
      30             :                         "Stabilization parameter in Nitsche's formulation and penalty factor "
      31             :                         "in the Penalty Method. In Nitsche's formulation this should be as "
      32             :                         "small as possible while the method is still stable; while in the "
      33             :                         "Penalty Method you want this to be quite large (e.g. 1e6).");
      34         768 :   params.addParam<FunctionName>("jump", 0, "Jump at the interface. Can be a Real or FunctionName.");
      35         768 :   params.addParam<FunctionName>(
      36         768 :       "jump_flux", 0, "Flux jump at the interface. Can be a Real or FunctionName.");
      37         768 :   params.addParam<UserObjectName>(
      38             :       "geometric_cut_userobject",
      39             :       "Name of GeometricCutUserObject associated with this constraint.");
      40         768 :   params.addParam<bool>(
      41             :       "use_penalty",
      42         768 :       false,
      43             :       "Use the Penalty instead of Nitsche (Nitsche only works for simple diffusion problems).");
      44         384 :   params.addClassDescription("Enforce constraints on the value or flux associated with a variable "
      45             :                              "at an XFEM interface.");
      46         384 :   return params;
      47           0 : }
      48             : 
      49         192 : XFEMSingleVariableConstraint::XFEMSingleVariableConstraint(const InputParameters & parameters)
      50             :   : ElemElemConstraint(parameters),
      51         384 :     _alpha(getParam<Real>("alpha")),
      52         192 :     _jump(getFunction("jump")),
      53         192 :     _jump_flux(getFunction("jump_flux")),
      54         576 :     _use_penalty(getParam<bool>("use_penalty"))
      55             : {
      56         576 :   _xfem = std::dynamic_pointer_cast<XFEM>(_fe_problem.getXFEM());
      57         192 :   if (_xfem == nullptr)
      58           0 :     mooseError("Problem casting to XFEM in XFEMSingleVariableConstraint");
      59             : 
      60             :   const UserObject * uo =
      61         384 :       &(_fe_problem.getUserObjectBase(getParam<UserObjectName>("geometric_cut_userobject")));
      62             : 
      63         192 :   if (dynamic_cast<const GeometricCutUserObject *>(uo) == nullptr)
      64           0 :     mooseError("UserObject casting to GeometricCutUserObject in XFEMSingleVariableConstraint");
      65             : 
      66         192 :   _interface_id = _xfem->getGeometricCutID(dynamic_cast<const GeometricCutUserObject *>(uo));
      67         192 : }
      68             : 
      69         576 : XFEMSingleVariableConstraint::~XFEMSingleVariableConstraint() {}
      70             : 
      71             : void
      72       62132 : XFEMSingleVariableConstraint::reinitConstraintQuadrature(const ElementPairInfo & element_pair_info)
      73             : {
      74       62132 :   _interface_normal = element_pair_info._elem1_normal;
      75       62132 :   ElemElemConstraint::reinitConstraintQuadrature(element_pair_info);
      76       62132 : }
      77             : 
      78             : Real
      79      964160 : XFEMSingleVariableConstraint::computeQpResidual(Moose::DGResidualType type)
      80             : {
      81             :   Real r = 0;
      82             : 
      83      964160 :   switch (type)
      84             :   {
      85      482080 :     case Moose::Element:
      86      482080 :       if (!_use_penalty)
      87             :       {
      88      418288 :         r -= (0.5 * _grad_u[_qp] * _interface_normal +
      89      418288 :               0.5 * _grad_u_neighbor[_qp] * _interface_normal) *
      90      418288 :              _test[_i][_qp];
      91      418288 :         r -= (_u[_qp] - _u_neighbor[_qp]) * 0.5 * _grad_test[_i][_qp] * _interface_normal;
      92      418288 :         r += 0.5 * _grad_test[_i][_qp] * _interface_normal * _jump.value(_t, _u[_qp]);
      93             :       }
      94      482080 :       r += 0.5 * _test[_i][_qp] * _jump_flux.value(_t, _u[_qp]);
      95      482080 :       r += _alpha * (_u[_qp] - _u_neighbor[_qp] - _jump.value(_t, _u[_qp])) * _test[_i][_qp];
      96      482080 :       break;
      97             : 
      98      482080 :     case Moose::Neighbor:
      99      482080 :       if (!_use_penalty)
     100             :       {
     101      418288 :         r += (0.5 * _grad_u[_qp] * _interface_normal +
     102      418288 :               0.5 * _grad_u_neighbor[_qp] * _interface_normal) *
     103      418288 :              _test_neighbor[_i][_qp];
     104      418288 :         r -= (_u[_qp] - _u_neighbor[_qp]) * 0.5 * _grad_test_neighbor[_i][_qp] * _interface_normal;
     105      418288 :         r += 0.5 * _grad_test_neighbor[_i][_qp] * _interface_normal *
     106      418288 :              _jump.value(_t, _u_neighbor[_qp]);
     107             :       }
     108      482080 :       r += 0.5 * _test_neighbor[_i][_qp] * _jump_flux.value(_t, _u_neighbor[_qp]);
     109      482080 :       r -= _alpha * (_u[_qp] - _u_neighbor[_qp] - _jump.value(_t, _u_neighbor[_qp])) *
     110      482080 :            _test_neighbor[_i][_qp];
     111      482080 :       break;
     112             :   }
     113      964160 :   return r;
     114             : }
     115             : 
     116             : Real
     117     1788416 : XFEMSingleVariableConstraint::computeQpJacobian(Moose::DGJacobianType type)
     118             : {
     119             :   Real r = 0;
     120             : 
     121     1788416 :   switch (type)
     122             :   {
     123      447104 :     case Moose::ElementElement:
     124      447104 :       if (!_use_penalty)
     125      395008 :         r += -0.5 * _grad_phi[_j][_qp] * _interface_normal * _test[_i][_qp] -
     126      395008 :              _phi[_j][_qp] * 0.5 * _grad_test[_i][_qp] * _interface_normal;
     127      447104 :       r += _alpha * _phi[_j][_qp] * _test[_i][_qp];
     128      447104 :       break;
     129             : 
     130      447104 :     case Moose::ElementNeighbor:
     131      447104 :       if (!_use_penalty)
     132      395008 :         r += -0.5 * _grad_phi_neighbor[_j][_qp] * _interface_normal * _test[_i][_qp] +
     133      395008 :              _phi_neighbor[_j][_qp] * 0.5 * _grad_test[_i][_qp] * _interface_normal;
     134      447104 :       r -= _alpha * _phi_neighbor[_j][_qp] * _test[_i][_qp];
     135      447104 :       break;
     136             : 
     137      447104 :     case Moose::NeighborElement:
     138      447104 :       if (!_use_penalty)
     139      395008 :         r += 0.5 * _grad_phi[_j][_qp] * _interface_normal * _test_neighbor[_i][_qp] -
     140      395008 :              _phi[_j][_qp] * 0.5 * _grad_test_neighbor[_i][_qp] * _interface_normal;
     141      447104 :       r -= _alpha * _phi[_j][_qp] * _test_neighbor[_i][_qp];
     142      447104 :       break;
     143             : 
     144      447104 :     case Moose::NeighborNeighbor:
     145      447104 :       if (!_use_penalty)
     146      395008 :         r += 0.5 * _grad_phi_neighbor[_j][_qp] * _interface_normal * _test_neighbor[_i][_qp] +
     147      395008 :              _phi_neighbor[_j][_qp] * 0.5 * _grad_test_neighbor[_i][_qp] * _interface_normal;
     148      447104 :       r += _alpha * _phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp];
     149      447104 :       break;
     150             :   }
     151             : 
     152     1788416 :   return r;
     153             : }