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 "EigenADReal.h"
      12             : #include "EqualValueEmbeddedConstraint.h"
      13             : #include "FEProblem.h"
      14             : #include "DisplacedProblem.h"
      15             : #include "AuxiliarySystem.h"
      16             : #include "SystemBase.h"
      17             : #include "Assembly.h"
      18             : #include "MooseMesh.h"
      19             : #include "AddVariableAction.h"
      20             : #include "MooseEnum.h"
      21             : 
      22             : #include "libmesh/sparse_matrix.h"
      23             : 
      24             : registerMooseObject("MooseApp", EqualValueEmbeddedConstraint);
      25             : registerMooseObject("MooseApp", ADEqualValueEmbeddedConstraint);
      26             : 
      27             : template <bool is_ad>
      28             : InputParameters
      29       43197 : EqualValueEmbeddedConstraintTempl<is_ad>::validParams()
      30             : {
      31       43197 :   MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
      32       43197 :   InputParameters params = GenericNodeElemConstraint<is_ad>::validParams();
      33       43197 :   params.addClassDescription("This is a constraint enforcing overlapping portions of two blocks to "
      34             :                              "have the same variable value");
      35       43197 :   params.set<bool>("use_displaced_mesh") = false;
      36       43197 :   MooseEnum formulation("kinematic penalty", "kinematic");
      37       43197 :   params.addParam<MooseEnum>(
      38             :       "formulation", formulation, "Formulation used to enforce the constraint");
      39       43197 :   params.addRequiredParam<Real>(
      40             :       "penalty",
      41             :       "Penalty parameter used in constraint enforcement for kinematic and penalty formulations.");
      42             : 
      43       86394 :   return params;
      44       43197 : }
      45             : 
      46             : template <bool is_ad>
      47         204 : EqualValueEmbeddedConstraintTempl<is_ad>::EqualValueEmbeddedConstraintTempl(
      48             :     const InputParameters & parameters)
      49             :   : GenericNodeElemConstraint<is_ad>(parameters),
      50         204 :     _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
      51         204 :     _fe_problem(*parameters.get<FEProblem *>("_fe_problem")),
      52         204 :     _formulation(this->template getParam<MooseEnum>("formulation").template getEnum<Formulation>()),
      53         204 :     _penalty(this->template getParam<Real>("penalty")),
      54         408 :     _residual_copy(_sys.residualGhosted())
      55             : {
      56         204 :   _overwrite_secondary_residual = false;
      57         204 :   prepareSecondaryToPrimaryMap();
      58             :   if constexpr (is_ad)
      59             :   {
      60          84 :     if (_formulation == Formulation::KINEMATIC)
      61           4 :       this->paramError("formulation", "AD constraints cannot be used with KINEMATIC formulation.");
      62             :   }
      63         200 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67         204 : EqualValueEmbeddedConstraintTempl<is_ad>::prepareSecondaryToPrimaryMap()
      68             : {
      69             :   // get mesh pointLocator
      70         204 :   std::unique_ptr<libMesh::PointLocatorBase> pointLocator = _mesh.getPointLocator();
      71         204 :   pointLocator->enable_out_of_mesh_mode();
      72         204 :   const std::set<subdomain_id_type> allowed_subdomains{_primary};
      73             : 
      74             :   // secondary id and primary id
      75             :   dof_id_type sid, mid;
      76             : 
      77             :   // prepare _secondary_to_primary_map
      78         204 :   std::set<dof_id_type> unique_secondary_node_ids;
      79         204 :   const MeshBase & meshhelper = _mesh.getMesh();
      80       19896 :   for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(_secondary),
      81         204 :                                     meshhelper.active_subdomain_elements_end(_secondary)))
      82             :   {
      83       55692 :     for (auto & sn : elem->node_ref_range())
      84             :     {
      85       45948 :       sid = sn.id();
      86       45948 :       if (_secondary_to_primary_map.find(sid) == _secondary_to_primary_map.end())
      87             :       {
      88             :         // primary element
      89       31208 :         const Elem * me = pointLocator->operator()(sn, &allowed_subdomains);
      90       31208 :         if (me != NULL)
      91             :         {
      92        7140 :           mid = me->id();
      93        7140 :           _secondary_to_primary_map.insert(std::pair<dof_id_type, dof_id_type>(sid, mid));
      94        7140 :           _subproblem.addGhostedElem(mid);
      95             :         }
      96             :       }
      97             :     }
      98             :   }
      99         204 : }
     100             : 
     101             : template <bool is_ad>
     102             : bool
     103       43736 : EqualValueEmbeddedConstraintTempl<is_ad>::shouldApply()
     104             : {
     105             :   // primary element
     106       43736 :   auto it = _secondary_to_primary_map.find(_current_node->id());
     107             : 
     108       43736 :   if (it != _secondary_to_primary_map.end())
     109             :   {
     110       19880 :     const Elem * primary_elem = _mesh.elemPtr(it->second);
     111       19880 :     std::vector<Point> points = {*_current_node};
     112             : 
     113             :     // reinit variables on the primary element at the secondary point
     114       19880 :     _fe_problem.setNeighborSubdomainID(primary_elem, 0);
     115       19880 :     _fe_problem.reinitNeighborPhys(primary_elem, points, 0);
     116             : 
     117       19880 :     reinitConstraint();
     118             : 
     119       19880 :     return true;
     120       19880 :   }
     121       23856 :   return false;
     122             : }
     123             : 
     124             : template <bool is_ad>
     125             : void
     126       19880 : EqualValueEmbeddedConstraintTempl<is_ad>::reinitConstraint()
     127             : {
     128       19880 :   const Node * node = _current_node;
     129       19880 :   unsigned int sys_num = _sys.number();
     130       19880 :   dof_id_type dof_number = node->dof_number(sys_num, _var.number(), 0);
     131             : 
     132       19880 :   switch (_formulation)
     133             :   {
     134        5320 :     case Formulation::KINEMATIC:
     135        5320 :       _constraint_residual = -_residual_copy(dof_number);
     136        5320 :       break;
     137             : 
     138       14560 :     case Formulation::PENALTY:
     139       14560 :       _constraint_residual = _penalty * (_u_secondary[0] - _u_primary[0]);
     140       14560 :       break;
     141             : 
     142           0 :     default:
     143           0 :       mooseError("Invalid formulation");
     144             :       break;
     145             :   }
     146       19880 : }
     147             : 
     148             : template <bool is_ad>
     149             : Real
     150        4970 : EqualValueEmbeddedConstraintTempl<is_ad>::computeQpSecondaryValue()
     151             : {
     152        4970 :   return MetaPhysicL::raw_value(_u_secondary[_qp]);
     153             : }
     154             : 
     155             : template <bool is_ad>
     156             : GenericReal<is_ad>
     157       85680 : EqualValueEmbeddedConstraintTempl<is_ad>::computeQpResidual(Moose::ConstraintType type)
     158             : {
     159       85680 :   GenericReal<is_ad> resid = _constraint_residual;
     160             : 
     161       85680 :   switch (type)
     162             :   {
     163       12040 :     case Moose::Secondary:
     164             :     {
     165       12040 :       if (_formulation == Formulation::KINEMATIC)
     166             :       {
     167        2660 :         GenericReal<is_ad> pen_force = _penalty * (_u_secondary[_qp] - _u_primary[_qp]);
     168        2660 :         resid += pen_force;
     169           0 :       }
     170       12040 :       return _test_secondary[_i][_qp] * resid;
     171             :     }
     172             : 
     173       73640 :     case Moose::Primary:
     174       73640 :       return _test_primary[_i][_qp] * -resid;
     175             :   }
     176             : 
     177           0 :   return 0.0;
     178       42420 : }
     179             : 
     180             : template <bool is_ad>
     181             : Real
     182      402150 : EqualValueEmbeddedConstraintTempl<is_ad>::computeQpJacobian(Moose::ConstraintJacobianType type)
     183             : {
     184             :   mooseAssert(!is_ad,
     185             :               "In ADEqualValueEmbeddedConstraint, computeQpJacobian should not be called. "
     186             :               "Check computeJacobian implementation.");
     187             : 
     188      402150 :   unsigned int sys_num = _sys.number();
     189      402150 :   const Real penalty = MetaPhysicL::raw_value(_penalty);
     190             :   Real curr_jac, secondary_jac;
     191             : 
     192      402150 :   switch (type)
     193             :   {
     194       30870 :     case Moose::SecondarySecondary:
     195       30870 :       switch (_formulation)
     196             :       {
     197       14490 :         case Formulation::KINEMATIC:
     198       14490 :           curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
     199       14490 :                                   _connected_dof_indices[_j]);
     200       14490 :           return -curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
     201       16380 :         case Formulation::PENALTY:
     202       16380 :           return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
     203           0 :         default:
     204           0 :           mooseError("Invalid formulation");
     205             :       }
     206             : 
     207       18760 :     case Moose::SecondaryPrimary:
     208       18760 :       switch (_formulation)
     209             :       {
     210        8960 :         case Formulation::KINEMATIC:
     211        8960 :           return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
     212        9800 :         case Formulation::PENALTY:
     213        9800 :           return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
     214           0 :         default:
     215           0 :           mooseError("Invalid formulation");
     216             :       }
     217             : 
     218      219240 :     case Moose::PrimarySecondary:
     219      219240 :       switch (_formulation)
     220             :       {
     221      105840 :         case Formulation::KINEMATIC:
     222      105840 :           secondary_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
     223      105840 :                                        _connected_dof_indices[_j]);
     224      105840 :           return secondary_jac * _test_primary[_i][_qp];
     225      113400 :         case Formulation::PENALTY:
     226      113400 :           return -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp];
     227           0 :         default:
     228           0 :           mooseError("Invalid formulation");
     229             :       }
     230             : 
     231      133280 :     case Moose::PrimaryPrimary:
     232      133280 :       switch (_formulation)
     233             :       {
     234       64960 :         case Formulation::KINEMATIC:
     235       64960 :           return 0.0;
     236       68320 :         case Formulation::PENALTY:
     237       68320 :           return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
     238           0 :         default:
     239           0 :           mooseError("Invalid formulation");
     240             :       }
     241             : 
     242           0 :     default:
     243           0 :       mooseError("Unsupported type");
     244             :       break;
     245             :   }
     246             :   return 0.0;
     247             : }
     248             : 
     249             : template <bool is_ad>
     250             : Real
     251        5600 : EqualValueEmbeddedConstraintTempl<is_ad>::computeQpOffDiagJacobian(
     252             :     Moose::ConstraintJacobianType type, unsigned int /*jvar*/)
     253             : {
     254             :   mooseAssert(!is_ad,
     255             :               "In ADEqualValueEmbeddedConstraint, computeQpOffDiagJacobian should not be called. "
     256             :               "Check computeJacobian implementation.");
     257             : 
     258             :   Real curr_jac, secondary_jac;
     259        5600 :   unsigned int sys_num = _sys.number();
     260             : 
     261        5600 :   switch (type)
     262             :   {
     263           0 :     case Moose::SecondarySecondary:
     264           0 :       curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
     265           0 :                               _connected_dof_indices[_j]);
     266           0 :       return -curr_jac;
     267             : 
     268        1120 :     case Moose::SecondaryPrimary:
     269        1120 :       return 0.0;
     270             : 
     271           0 :     case Moose::PrimarySecondary:
     272           0 :       switch (_formulation)
     273             :       {
     274           0 :         case Formulation::KINEMATIC:
     275           0 :           secondary_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
     276           0 :                                        _connected_dof_indices[_j]);
     277           0 :           return secondary_jac * _test_primary[_i][_qp];
     278           0 :         case Formulation::PENALTY:
     279           0 :           return 0.0;
     280           0 :         default:
     281           0 :           mooseError("Invalid formulation");
     282             :       }
     283             : 
     284        4480 :     case Moose::PrimaryPrimary:
     285        4480 :       return 0.0;
     286             : 
     287           0 :     default:
     288           0 :       mooseError("Unsupported type");
     289             :       break;
     290             :   }
     291             : 
     292             :   return 0.0;
     293             : }
     294             : 
     295             : template class EqualValueEmbeddedConstraintTempl<false>;
     296             : template class EqualValueEmbeddedConstraintTempl<true>;