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