LCOV - code coverage report
Current view: top level - src/constraints - MortarConstraintBase.C (source / functions) Hit Total Coverage
Test: idaholab/moose framework: 2bf808 Lines: 108 126 85.7 %
Date: 2025-07-17 01:28:37 Functions: 6 6 100.0 %
Legend: Lines: hit not hit

          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 "MortarConstraintBase.h"
      11             : #include "FEProblemBase.h"
      12             : #include "Assembly.h"
      13             : #include "MooseVariableFE.h"
      14             : 
      15             : #include "libmesh/string_to_enum.h"
      16             : 
      17             : InputParameters
      18      216431 : MortarConstraintBase::validParams()
      19             : {
      20      216431 :   InputParameters params = Constraint::validParams();
      21      216431 :   params += MortarConsumerInterface::validParams();
      22      216431 :   params += TwoMaterialPropertyInterface::validParams();
      23             : 
      24             :   // Whether on a displaced or undisplaced mesh, coupling ghosting will only happen for
      25             :   // cross-interface elements
      26      216431 :   params.addRelationshipManager("AugmentSparsityOnInterface",
      27             :                                 Moose::RelationshipManagerType::COUPLING,
      28           0 :                                 [](const InputParameters & obj_params, InputParameters & rm_params)
      29             :                                 {
      30        3688 :                                   rm_params.set<bool>("use_displaced_mesh") =
      31        3688 :                                       obj_params.get<bool>("use_displaced_mesh");
      32        7376 :                                   rm_params.set<BoundaryName>("secondary_boundary") =
      33        7376 :                                       obj_params.get<BoundaryName>("secondary_boundary");
      34        7376 :                                   rm_params.set<BoundaryName>("primary_boundary") =
      35        7376 :                                       obj_params.get<BoundaryName>("primary_boundary");
      36        7376 :                                   rm_params.set<SubdomainName>("secondary_subdomain") =
      37        7376 :                                       obj_params.get<SubdomainName>("secondary_subdomain");
      38        7376 :                                   rm_params.set<SubdomainName>("primary_subdomain") =
      39        7376 :                                       obj_params.get<SubdomainName>("primary_subdomain");
      40        3688 :                                   rm_params.set<bool>("ghost_higher_d_neighbors") =
      41        3688 :                                       obj_params.get<bool>("ghost_higher_d_neighbors");
      42        3688 :                                 });
      43             : 
      44             :   // If the LM is ever a Lagrange variable, we will attempt to obtain its dof indices from the
      45             :   // process that owns the node. However, in order for the dof indices to get set for the Lagrange
      46             :   // variable, the process that owns the node needs to have local copies of any lower-d elements
      47             :   // that have the connected node. Note that the geometric ghosting done here is different than that
      48             :   // done by the AugmentSparsityOnInterface RM, even when ghost_point_neighbors is true. The latter
      49             :   // ghosts equal-manifold lower-dimensional secondary element point neighbors and their interface
      50             :   // couplings. This ghosts lower-dimensional point neighbors of higher-dimensional elements.
      51             :   // Neither is guaranteed to be a superset of the other. For instance ghosting of lower-d point
      52             :   // neighbors (AugmentSparsityOnInterface with ghost_point_neighbors = true) is only guaranteed to
      53             :   // ghost those lower-d point neighbors on *processes that own lower-d elements*. And you may have
      54             :   // a process that only owns higher-dimensional elements
      55             :   //
      56             :   // Note that in my experience it is only important for the higher-d lower-d point neighbors to be
      57             :   // ghosted when forming sparsity patterns and so I'm putting this here instead of at the
      58             :   // MortarConsumerInterface level
      59      216431 :   params.addRelationshipManager("GhostHigherDLowerDPointNeighbors",
      60             :                                 Moose::RelationshipManagerType::GEOMETRIC);
      61             : 
      62      216431 :   params.addParam<VariableName>("secondary_variable", "Primal variable on secondary surface.");
      63      216431 :   params.addParam<VariableName>(
      64             :       "primary_variable",
      65             :       "Primal variable on primary surface. If this parameter is not provided then the primary "
      66             :       "variable will be initialized to the secondary variable");
      67      216431 :   params.makeParamNotRequired<NonlinearVariableName>("variable");
      68      216431 :   params.setDocString(
      69             :       "variable",
      70             :       "The name of the lagrange multiplier variable that this constraint is applied to. This "
      71             :       "parameter may not be supplied in the case of using penalty methods for example");
      72      649293 :   params.addParam<bool>(
      73      432862 :       "compute_primal_residuals", true, "Whether to compute residuals for the primal variable.");
      74      649293 :   params.addParam<bool>(
      75      432862 :       "compute_lm_residuals", true, "Whether to compute Lagrange Multiplier residuals");
      76      649293 :   params.addParam<MooseEnum>(
      77             :       "quadrature",
      78      432862 :       MooseEnum("DEFAULT FIRST SECOND THIRD FOURTH", "DEFAULT"),
      79             :       "Quadrature rule to use on mortar segments. For 2D mortar DEFAULT is recommended. "
      80             :       "For 3D mortar, QUAD meshes are integrated using triangle mortar segments. "
      81             :       "While DEFAULT quadrature order is typically sufficiently accurate, exact integration of "
      82             :       "QUAD mortar faces requires SECOND order quadrature for FIRST variables and FOURTH order "
      83             :       "quadrature for SECOND order variables.");
      84      649293 :   params.addParam<bool>(
      85             :       "use_petrov_galerkin",
      86      432862 :       false,
      87             :       "Whether to use the Petrov-Galerkin approach for the mortar-based constraints. If set to "
      88             :       "true, we use the standard basis as the test function and dual basis as "
      89             :       "the shape function for the interpolation of the Lagrange multiplier variable.");
      90      216431 :   params.addCoupledVar("aux_lm",
      91             :                        "Auxiliary Lagrange multiplier variable that is utilized together with the "
      92             :                        "Petrov-Galerkin approach.");
      93      216431 :   return params;
      94           0 : }
      95             : 
      96        1232 : MortarConstraintBase::MortarConstraintBase(const InputParameters & parameters)
      97             :   : Constraint(parameters),
      98             :     NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, false, false),
      99             :     MortarConsumerInterface(this),
     100             :     TwoMaterialPropertyInterface(this, Moose::EMPTY_BLOCK_IDS, getBoundaryIDs()),
     101             :     MooseVariableInterface<Real>(this,
     102             :                                  true,
     103        2464 :                                  isParamValid("variable") ? "variable" : "secondary_variable",
     104             :                                  Moose::VarKindType::VAR_SOLVER,
     105             :                                  Moose::VarFieldType::VAR_FIELD_STANDARD),
     106        1232 :     _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
     107        2464 :     _var(isParamValid("variable")
     108        1232 :              ? &_subproblem.getStandardVariable(_tid, parameters.getMooseType("variable"))
     109             :              : nullptr),
     110        1232 :     _secondary_var(
     111        2464 :         isParamValid("secondary_variable")
     112        2464 :             ? _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("secondary_variable"))
     113        1232 :             : _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("primary_variable"))),
     114        1232 :     _primary_var(
     115        2464 :         isParamValid("primary_variable")
     116        1232 :             ? _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("primary_variable"))
     117             :             : _secondary_var),
     118             : 
     119        1232 :     _compute_primal_residuals(getParam<bool>("compute_primal_residuals")),
     120        1232 :     _compute_lm_residuals(!_var ? false : getParam<bool>("compute_lm_residuals")),
     121        1232 :     _test_dummy(),
     122        1232 :     _use_dual(_var ? _var->useDual() : false),
     123        1232 :     _tangents(_assembly.tangents()),
     124        1232 :     _coord(_assembly.mortarCoordTransformation()),
     125        1232 :     _q_point(_assembly.qPointsMortar()),
     126        1232 :     _use_petrov_galerkin(getParam<bool>("use_petrov_galerkin")),
     127        1232 :     _aux_lm_var(isCoupled("aux_lm") ? getVar("aux_lm", 0) : nullptr),
     128        2464 :     _test(_var
     129        1232 :               ? ((_use_petrov_galerkin && _aux_lm_var) ? _aux_lm_var->phiLower() : _var->phiLower())
     130             :               : _test_dummy),
     131        1232 :     _test_secondary(_secondary_var.phiFace()),
     132        1232 :     _test_primary(_primary_var.phiFaceNeighbor()),
     133        1232 :     _grad_test_secondary(_secondary_var.gradPhiFace()),
     134        1232 :     _grad_test_primary(_primary_var.gradPhiFaceNeighbor()),
     135        1232 :     _interior_secondary_elem(_assembly.elem()),
     136        1232 :     _interior_primary_elem(_assembly.neighbor()),
     137        4928 :     _displaced(getParam<bool>("use_displaced_mesh"))
     138             : {
     139        1232 :   if (_use_dual)
     140         124 :     _assembly.activateDual();
     141             : 
     142        1232 :   if (_use_petrov_galerkin && (!_use_dual))
     143           0 :     paramError("use_petrov_galerkin",
     144             :                "We need to set `use_dual = true` while using the Petrov-Galerkin approach");
     145             : 
     146        1232 :   if (_use_petrov_galerkin && ((!isParamValid("aux_lm")) || _aux_lm_var == nullptr))
     147           0 :     paramError("use_petrov_galerkin",
     148             :                "We need to specify an auxiliary variable `aux_lm` while using the Petrov-Galerkin "
     149             :                "approach");
     150             : 
     151        1232 :   if (_use_petrov_galerkin && _aux_lm_var->useDual())
     152           0 :     paramError("aux_lm",
     153             :                "Auxiliary LM variable needs to use standard shape function, i.e., set `use_dual = "
     154             :                "false`.");
     155             : 
     156             :   // Note parameter is discretization order, we then convert to quadrature order
     157        1232 :   const MooseEnum p_order = getParam<MooseEnum>("quadrature");
     158             :   // If quadrature not DEFAULT, set mortar qrule
     159        1232 :   if (p_order != "DEFAULT")
     160             :   {
     161           0 :     Order q_order = static_cast<Order>(2 * Utility::string_to_enum<Order>(p_order) + 1);
     162           0 :     _assembly.setMortarQRule(q_order);
     163             :   }
     164             : 
     165        1232 :   if (_var)
     166         668 :     addMooseVariableDependency(_var);
     167        1232 :   addMooseVariableDependency(&_secondary_var);
     168        1232 :   addMooseVariableDependency(&_primary_var);
     169        1232 : }
     170             : 
     171             : void
     172      188598 : MortarConstraintBase::computeResidual()
     173             : {
     174      188598 :   precalculateResidual();
     175             : 
     176      188598 :   if (_compute_primal_residuals)
     177             :   {
     178             :     // Compute the residual for the secondary interior primal dofs
     179      188598 :     computeResidual(Moose::MortarType::Secondary);
     180             : 
     181             :     // Compute the residual for the primary interior primal dofs.
     182      188598 :     computeResidual(Moose::MortarType::Primary);
     183             :   }
     184             : 
     185      188598 :   if (_compute_lm_residuals)
     186             :     // Compute the residual for the lower dimensional LM dofs (if we even have an LM variable)
     187       89334 :     computeResidual(Moose::MortarType::Lower);
     188      188598 : }
     189             : 
     190             : void
     191       89628 : MortarConstraintBase::computeJacobian()
     192             : {
     193       89628 :   precalculateResidual();
     194             : 
     195       89628 :   if (_compute_primal_residuals)
     196             :   {
     197             :     // Compute the jacobian for the secondary interior primal dofs
     198       89628 :     computeJacobian(Moose::MortarType::Secondary);
     199             : 
     200             :     // Compute the jacobian for the primary interior primal dofs.
     201       89628 :     computeJacobian(Moose::MortarType::Primary);
     202             :   }
     203             : 
     204       89628 :   if (_compute_lm_residuals)
     205             :     // Compute the jacobian for the lower dimensional LM dofs (if we even have an LM variable)
     206       48972 :     computeJacobian(Moose::MortarType::Lower);
     207       89628 : }
     208             : 
     209             : void
     210       15741 : MortarConstraintBase::zeroInactiveLMDofs(const std::unordered_set<const Node *> & inactive_lm_nodes,
     211             :                                          const std::unordered_set<const Elem *> & inactive_lm_elems)
     212             : {
     213             :   // If no LM variable has been defined, skip
     214       15741 :   if (!_var)
     215        6503 :     return;
     216             : 
     217        9238 :   const auto sn = _sys.number();
     218        9238 :   const auto vn = _var->number();
     219             : 
     220             :   // If variable is nodal, zero DoFs based on inactive LM nodes
     221        9238 :   if (_var->isNodal())
     222             :   {
     223        6058 :     for (const auto node : inactive_lm_nodes)
     224             :     {
     225             :       // Allow mixed Lagrange orders between primal and LM
     226           0 :       if (!node->n_comp(sn, vn))
     227           0 :         continue;
     228             : 
     229           0 :       const auto dof_index = node->dof_number(sn, vn, 0);
     230             :       // No scaling; this is not physics
     231           0 :       if (_assembly.computingJacobian())
     232           0 :         addJacobianElement(
     233             :             _assembly, /*element_value=*/1, dof_index, dof_index, /*scaling_factor=*/1);
     234           0 :       if (_assembly.computingResidual())
     235             :       {
     236           0 :         const Real lm_value = _var->getNodalValue(*node);
     237           0 :         addResiduals(_assembly,
     238           0 :                      std::array<Real, 1>{{lm_value}},
     239           0 :                      std::array<dof_id_type, 1>{{dof_index}},
     240             :                      /*scaling_factor=*/1);
     241             :       }
     242             :     }
     243             :   }
     244             :   // If variable is elemental, zero based on inactive LM elems
     245             :   else
     246             :   {
     247        3540 :     for (const auto el : inactive_lm_elems)
     248             :     {
     249         360 :       const auto n_comp = el->n_comp(sn, vn);
     250             : 
     251         720 :       for (const auto comp : make_range(n_comp))
     252             :       {
     253         360 :         const auto dof_index = el->dof_number(sn, vn, comp);
     254             :         // No scaling; this is not physics
     255         360 :         if (_assembly.computingJacobian())
     256         120 :           addJacobianElement(
     257             :               _assembly, /*element_value=*/1, dof_index, dof_index, /*scaling_factor=*/1);
     258         360 :         if (_assembly.computingResidual())
     259             :         {
     260         240 :           const Real lm_value = _var->getElementalValue(el, comp);
     261         240 :           addResiduals(_assembly,
     262           0 :                        std::array<Real, 1>{{lm_value}},
     263         240 :                        std::array<dof_id_type, 1>{{dof_index}},
     264             :                        /*scaling_factor=*/1);
     265             :         }
     266             :       }
     267             :     }
     268             :   }
     269             : }

Generated by: LCOV version 1.14