MortarScalarBase

The MortarScalarBase is part of the scalar augmentation system to complement the ScalarKernel class. Its principal purpose is to add standard quadrature loops and assembly routines to handle the contributions from a single added scalar variable to a Mortar Constraint class, including the entire row of the Jacobian. These routines will assist with representing weak form terms in a partial differential equation involving scalar variables integrated over the surface of a non-conforming interface (mortar surface). As usual, this piece of physics is referred to as the "residual" and is evaluated at integration quadrature points along that interface. To implement your own physics in MOOSE, you create your own mortar constraint by subclassing the MOOSE MortarScalarBase class.

The mortar scalar augmentation system supports the use of AD for residual calculations, as such there are two options for creating field-scalar coupling objects: MortarScalarBase and ADMortarScalarBase. To further understand automatic differentiation, please refer to the Automatic Differentiation page for more information.

Developers should read all sections; users can find Parameters described at the bottom.

Creation of Mortar Constraint Scalar Coupling Classes

MortarConstraint objects are designed to handle weak-form contributions from both surface-based lower-dimensional (Lagrange multiplier) variables and domain-based (primary and secondary) spatial variables. For these three focus spatial variables, its job is to contribute to the associated rows of the residual and the Jacobian matrix. Herein, as in the source code of MortarConstraints, these spatial variables will be called (see Mortar Parameters for the analogous input parameter names):

_var: Corresponds to a Lagrange Multiplier variable that lives on the lower dimensional block on the secondary face
_secondary_var: Primal variable on the secondary side of the mortar interface (lives on the interior elements)
_primary_var: Primal variable on the primary side of the mortar interface (lives on the interior elements).

The three focus spatial variables _var, _secondary_var, and _primary_var are each indicated by the name var below, and are differentiated by the mortar_type flag as one of "Primary". "Secondary", or "Lower". In a coupled (multi-phyics) weak form, all interface integral terms containing the test function of one of these three variables are potential candidates for MortarConstraint contributions.

The philosphy of the scalar augmentation class MortarScalarBase is to add a single focus scalar variable referred to as _kappa to the MortarConstraint object so that all terms in the coupled weak form that involve the spatial variables, scalar variable, and/or their test functions can be assembled in one or multiple class instances. This philosophy is similar to how the lower dimensional variable _lambda is added to the element faces of DGKernel and IntegratedBC objects associated with the hybrid finite element method (HFEM). Documentation for that approach can be found HFEMDiffusion and HFEMDirichletBC along with the base classes DGLowerDKernel and LowerDIntegratedBC.

In a MortarScalarBase subclass, a naming scheme is established for the quadrature point methods of the two variable types: methods contributing to the test function of _kappa have "Scalar" near the front and methods contributing to the trial function of scalar variables in the Jacobian have "Scalar" at the end. The computeScalarQpResidual() function should be overridden (see Parameters for cases where the scalar should be suppressed). The computeQpResidual(mortar_type) function must be overridden as usual for MortarConstraint, although it may return zero.

For non-AD objects, several contributions to the Jacobian matrix can be optionally overridden for use in Newton-type nonlinear solvers. As mentioned later, the developer should choose and document which terms (rows) of the residual and terms (rows and columns) of the Jacobian will be attributed to an instance of the developed class. These choices can be motivated by whether some terms in the weak form can be or have already been implemented within other MOOSE classes.

computeQpJacobian(jacobian_type, jvar_num): Jacobian component d-var-residual / d-jvar
computeQpOffDiagJacobianScalar(mortar_type, svar_num): Jacobian component d-var-residual / d-svar
computeScalarQpJacobian(): Jacobian component d-_kappa-residual / d-_kappa
computeScalarQpOffDiagJacobian(mortar_type, jvar_num): Jacobian component d-_kappa-residual / d-jvar
computeScalarQpOffDiagJacobianScalar(svar_num): Jacobian component d-_kappa-residual / d-svar

Similar to the mortar_type flag mentioned above, the jacobian_type flag is used to distinguish the couplings between the test and trial functions of the three focus spatial variables. For example, jacobian_type = LowerSecondary indicates that the linearized weak form term from the _var test function and _secondary_var trial function should be evaluated. All nine combinations from the three focus variables are visited during the loops that call computeQpJacobian(jacobian_type, jvar_num). Also, loops over the coupled variables wrap around these quadrature loops: the integer for the spatial variable is jvar_num and the integer for the scalar variable is svar_num. Note that the jvar_num integer is different for _var and _secondary_var and for any other spatial variables coupled by the developer through derived classes of MortarScalarBase, so consult the examples of these methods below in Examples from Source Code for how to query the current target using logical tests.

Also, there are some pre-calculation routines that are called within the quadrature loop once before the loop over spatial variable test and shape functions as well as before the loop over scalar components. These methods are useful for material or stabilization calculations.

initScalarQpResidual(): evaluations depending on qp but independent of test functions
initScalarQpJacobian(jvar_num): evaluations depending on qp but independent of test and shape functions
initScalarQpOffDiagJacobian(jsvar): evaluations depending on qp but independent of test and shape functions

In addition to those mentioned in the MortarConstraints documentation, you have access to several member variables inside your MortarScalarBase class for computing the residual and Jacobian values in the above mentioned functions:

_h, _l: indices for the current test and trial scalar component respectively.
_kappa: value of the scalar variable this MortarConstraint operates on; indexed by _h (i.e. _kappa[_h]).
_kappa_var: ID number of this scalar variable; useful to differentiate from others.
_k_order: order (number of components) of the scalar variable.

Since the test and trial "shape" functions of a scalar are "1", variables are not required for that value. Examples of the source codes below demonstrate this fact.

warning:AD global indexing required

ADMortarScalarBase only works with MOOSE configured with global AD indexing (the default).

Examples from Source Code

As mentioned, the computeScalarQpResidual method should be overridden for both flavors of mortar constraints, non-AD and AD. As an example, consider the scalar residual weak form terms of the PeriodicSegmentalConstraint class:

$(1)var element = document.getElementById("moose-equation-399765d1-d08f-4203-b02f-2c468d4bac49");katex.render(" \\langle \\vec{\\kappa} \\cdot \\llbracket \\vec{x} \\rrbracket, \\lambda \\rangle _{\\Gamma^+} - \\langle \\vec{\\kappa} \\cdot \\llbracket \\vec{x} \\rrbracket, \\vec{\\sigma} \\cdot \\hat{n} \\rangle _{\\Gamma^+} ", element, {displayMode:true,throwOnError:false});$

The computeScalarQpResidual method for the non-AD version of this class is provided in Listing 1, where _kappa_aux is equal to $var element = document.getElementById("moose-equation-1a7f6059-4f47-445d-83c2-74a5479b5259");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ .

Listing 1: The C++ weak-form residual statement of Eq. (1).

Real
PeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }

Meanwhile, the contribution to the lower spatial variable residual of this object is associated with Eq. (2) and implemented in Listing 2 (note that the scalar variable _kappa is termed as $var element = document.getElementById("moose-equation-5ab7b25b-eafd-4f81-bb00-8d1d1afdeac9");katex.render("\\vec{\\epsilon}", element, {displayMode:false,throwOnError:false});$ in this weak form).

$(2)var element = document.getElementById("moose-equation-59032b99-23cc-4ea9-b762-1e163fb554a6");katex.render(" \\langle \\mu , \\vec{\\epsilon} \\cdot \\llbracket \\vec{x} \\rrbracket \\rangle _{\\Gamma^+} ", element, {displayMode:true,throwOnError:false});$

Listing 2: The C++ weak-form residual statement of Eq. (2).

Real
PeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }

The other contributions to the lower spatial variable residual of this object is associated with Eq. (3) and implemented in Listing 3 within the EqualValueConstraint class.

$(3)var element = document.getElementById("moose-equation-d51774ae-4d2a-42d4-a5bf-fa26b8b2a630");katex.render(" -\\langle \\mu , \\llbracket u \\rrbracket \\rangle _{\\Gamma^+} ", element, {displayMode:true,throwOnError:false});$

Listing 3: The C++ weak-form residual statement of Eq. (3).

case Moose::MortarType::Lower:
    {
      auto residual = (_u_primary[_qp] - _u_secondary[_qp]) * _test[_i][_qp];

      if (_stabilize)
      {
        // secondary
        residual -= _delta * _lower_secondary_volume * _test[_i][_qp] *
                    (_lambda[_qp] - _diff_secondary[_qp] * _grad_u_secondary[_qp] * _normals[_qp]);

        // primary
        residual -=
            _delta * _lower_primary_volume * _test[_i][_qp] *
            (_lambda[_qp] + _diff_primary[_qp] * _grad_u_primary[_qp] * _normals_primary[_qp]);
      }

For an example of the contributions to _primary_var and _secondary_var residuals, consider the penalty version of the periodic constraint associated with Eq. (4) and implemented in Listing 4 within the PenaltyPeriodicSegmentalConstraint class.

$(4)var element = document.getElementById("moose-equation-b766f396-3c4f-46df-b493-c6f3ac1b4756");katex.render(" - \\langle \\llbracket w \\rrbracket,\\tau \\vec{\\epsilon} \\cdot \\llbracket \\vec{x} \\rrbracket \\rangle _{\\Gamma^+} ", element, {displayMode:true,throwOnError:false});$

Listing 4: The C++ weak-form residual statement of Eq. (4).

Real
PenaltyPeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  // Compute penalty parameter times x-jump times average heat flux
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }

The PeriodicSegmentalConstraint class also overrides the computeScalarQpOffDiagJacobian method to define the Jacobian term related to Eq. (1) as shown in Listing 5.

Listing 5: The C++ weak-form Jacobian for d-_kappa-residual / d-jvar.

Real
PeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(const Moose::MortarType mortar_type,
                                                            const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      if (!_var || _var->number() != jvar_num)
        return 0;
      break;
    default:
      return 0;
  }

Notice that there is a conditional to confirm that the coupled jvar is the focus lower variable _var, otherwise it returns zero.

Similarly, it also overrides the computeQpOffDiagJacobianScalar method to define the Jacobian term related to Eq. (2) as shown in Listing 6.

Listing 6: The C++ weak-form Jacobian for d-_var-residual / d-svar.

Real
PeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(const Moose::MortarType mortar_type,
                                                            const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test[_i][_qp];
      break;
    default:
      return 0;
  }

Also notice the conditional that confirms the coupled svar is the focus scalar _kappa, otherwise it returns zero.

The AD version of this object, ADPeriodicSegmentalConstraint, only requires the residual implementation; as such it overrides computeScalarQpResidual and computeQpResidual as follows.

Listing 7: The C++ AD weak-form residual statement of Eq. (1).

ADPeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADReal r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }

Listing 8: The C++ AD weak-form residual statement of Eq. (2).

ADPeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADRealVectorValue kappa_vec(_kappa[0], 0, 0);
  Moose::derivInsert(kappa_vec(0).derivatives(), _kappa_var_ptr->dofIndices()[0], 1);
  if (_k_order == 2)
  {
    kappa_vec(1) = _kappa[1];
    Moose::derivInsert(kappa_vec(1).derivatives(), _kappa_var_ptr->dofIndices()[1], 1);
  }

Depending upon the weak form and its coupling terms between spatial and scalar variables, not all of the methods listed in Creation of Mortar Constraint Scalar Coupling Classes need to be overridden.

The scalar augmentation system is designed such that multiple scalar variables can be coupled to an instance of the MortarConstraint class, each focusing on one scalar from the list. This approach is similar to how Tensor Mechanics module classes operator on one component variable of the displacement vector field and are coupled to the other components. The developer can decide how to organize the coupling and off-diagonal Jacobian terms in a logical way and document this for the user.

An example for decomposing the coupling terms and having multiple scalar variables are contained in the source files of the framework test directory as well as input file testperiodicsole.i, with listings below. The comments within these header and source files should be consulted to visualize how the rows and columns of the relevant residual and Jacobian contributions are handled.

Listing 9: Decomposition of spatial and scalar variable contributions by x and y components.

/// Test object to illustrate coupling between mortar spatial variables
/// and two separate scalar variables, kappa and kappa_other. The solution
/// is equivalent to PenaltyPeriodicSegmentalConstraint.
///
/// The kappa variable is split into two scalars: the first component called '_kappa'
/// herein and all other components called '_kappa_other' herein.
/// This decomposition is ONLY valid for 2D problems with a macro scalar vector [kappa_x kappa_y]
/// applied as a periodic boundary condition. For parameter _alpha = 0, the primary
/// scalar (_kappa) is kappa_x and the coupled scalar is kappa_y. For parameter _alpha = 1, the primary
/// scalar (_kappa) is kappa_y and the coupled scalar is kappa_x.
/// That mapping MUST be handled in the input file.
///
/// Thus, each instance of TestPeriodicSole acts on both primary/secondary field variables
/// and one scalar variable (_svar_alpha). The job of the constraint is to assemble the
/// residual and Jacobian terms arising from _svar_alpha (namely, selected entries).
/// It assembles one row of the Jacobian for the scalar variable as well as the couplings of that
/// scalar to the primary and secondary spatial variables, using logical checks with _alpha.
/// The entries for the other field/scalar variables are handled by the other instance of the
/// constraint, which have the opposite value of _alpha. The logical checks ensure the proper
/// decomposition of the jobs.
///
/// In summary, for x=disp_x etc. and k=kappa_x and a=kappa_y, then the contributions of the instances are
/// PenaltyEqualValueConstraint
/// R = [Rss+p, Rps+p, 0, 0]^T
/// J = [Jss, Jsp, 000, 000
///      Jps, Jpp, 000, 000
///      000, 000, 000, 000
///      000, 000, 000, 000]
/// TestPeriodicSole
/// _alpha=0
/// R = [Rsk, Rpk,  Rk,  00]^T
/// J = [000, 000, Jsk, 000
///      000, 000, Jpk, 000
///      Jks, Jkp, Jkk, Jka
///      000, 000, 000, 000]
/// _alpha=1
/// R = [Rsa, Rpa,  00,  Ra]^T
/// J = [000, 000, 000, Jsa
///      000, 000, 000, Jpa
///      000, 000, 000, 000
///      Jas, Jap, Jak, Jaa]
///
/// In this manner, the full R and J are obtained with NO duplication of jobs:
/// R = [Rs,  Rp,  Rk,  Ra ]^T
/// J = [Jss, Jsp, Jsk, Jsa
///      Jps, Jpp, Jpk, Jpa
///      Jks, Jkp, Jkk, Jka
///      Jas, Jap, Jak, Jaa]
///

class TestPeriodicSole : public DerivativeMaterialInterface<MortarScalarBase>
{
public:
  static InputParameters validParams();
  TestPeriodicSole(const InputParameters & parameters);

protected:
  virtual void precalculateResidual() override;
  virtual void precalculateJacobian() override;
  virtual void initScalarQpResidual() override;

  /**
   * Method for computing the residual at quadrature points
   */
  virtual Real computeQpResidual(Moose::MortarType mortar_type) override;
  Real computeQpJacobian(Moose::ConstraintJacobianType /*jacobian_type*/,
                         unsigned int /*jvar*/) override
  {
    return 0;
  }

[Constraints]
  [mortarlr]
    type = PenaltyEqualValueConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodiclrx]
    type = TestPeriodicSole
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    kappa = kappa_x
    kappa_aux = kappa_aux
    component = 0
    kappa_other = kappa_y
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodiclry]
    type = TestPeriodicSole
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    kappa = kappa_y
    kappa_aux = kappa_aux
    component = 1
    kappa_other = kappa_x
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [mortarbt]
    type = PenaltyEqualValueConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodicbtx]
    type = TestPeriodicSole
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    kappa = kappa_x
    kappa_aux = kappa_aux
    component = 0
    kappa_other = kappa_y
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodicbty]
    type = TestPeriodicSole
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    kappa = kappa_y
    kappa_aux = kappa_aux
    component = 1
    kappa_other = kappa_x
    correct_edge_dropping = true
    compute_scalar_residuals = true
    penalty_value = 1.e3
  []
[]

note:Displaced mesh features untested

The displaced mesh features are not yet tested for the scalar augmentation system.

Parameters

There is one required parameters the user must supply for a mortar constraint derived from MortarScalarBase:

scalar_variable: the focus scalar variable of the mortar constraint, for which assembly of the residual and Jacobian contributions will occur. It must be a MooseScalarVariable. This parameter may be renamed in a derived class to be more physically meaningful.

If the scalar_variable parameter is not specified, then the derived class will behave identically to a regular MortarConstraint, namely without any scalar functionality. This feature is useful if the scalar augmentation in inserted into a class structure with several levels and not all derived members use scalar variables.

As an example, the parameter listing is shown below for the PeriodicSegmentalConstraint object with the scalar_variable parameter renamed to epsilon:

[Constraints]
  [mortarlr]
    type = EqualValueConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    variable = lm1
    correct_edge_dropping = true
  []
  [periodiclr]
    type = PeriodicSegmentalConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm1
    correct_edge_dropping = true
  []
  [mortarbt]
    type = EqualValueConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    variable = lm2
    correct_edge_dropping = true
  []
  [periodicbt]
    type = PeriodicSegmentalConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm2
    correct_edge_dropping = true
  []
[]

Note: to avoid an error message "Variable 'kappa' does not exist in this system", the following block should be added to the input file:

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

There is one optional parameters that can be supplied to`MortarScalarBase` classes:

compute_scalar_residuals: Whether to compute scalar residuals. This will automatically be set to false if a scalar_variable parameter is not supplied. Other cases where the user may want to set this to false is during testing when the scalar variable is an AuxVariable and not a solution variable in the system. This parameter has a similar usage as the compute_lm_residuals and compute_primal_residuals for all Mortar objects.

(moose/framework/src/constraints/PeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", PeriodicSegmentalConstraint);

InputParameters
PeriodicSegmentalConstraint::validParams()
{
  InputParameters params = MortarScalarBase::validParams();
  params.addClassDescription(
      "PeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

PeriodicSegmentalConstraint::PeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<MortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

Real
PeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }
  Real r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

Real
PeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

Real
PeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(const Moose::MortarType mortar_type,
                                                            const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

Real
PeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(const Moose::MortarType mortar_type,
                                                            const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      if (!_var || _var->number() != jvar_num)
        return 0;
      break;
    default:
      return 0;
  }

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      jac *= (*_phi)[_j][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

(moose/framework/src/constraints/PeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", PeriodicSegmentalConstraint);

InputParameters
PeriodicSegmentalConstraint::validParams()
{
  InputParameters params = MortarScalarBase::validParams();
  params.addClassDescription(
      "PeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

PeriodicSegmentalConstraint::PeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<MortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

Real
PeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }
  Real r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

Real
PeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

Real
PeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(const Moose::MortarType mortar_type,
                                                            const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

Real
PeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(const Moose::MortarType mortar_type,
                                                            const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      if (!_var || _var->number() != jvar_num)
        return 0;
      break;
    default:
      return 0;
  }

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      jac *= (*_phi)[_j][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

(moose/framework/src/constraints/EqualValueConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "EqualValueConstraint.h"
#include "SubProblem.h"
#include "FEProblem.h"
#include "Assembly.h"

registerMooseObject("MooseApp", EqualValueConstraint);

InputParameters
EqualValueConstraint::validParams()
{
  InputParameters params = ADMortarConstraint::validParams();
  params.addClassDescription(
      "EqualValueConstraint enforces solution continuity between secondary and "
      "primary sides of a mortar interface using lagrange multipliers");
  params.addRangeCheckedParam<Real>(
      "delta", 0, "0<=delta<=1", "The coefficient for stabilizing terms");
  params.addParam<MaterialPropertyName>(
      "diff_secondary", 1, "The diffusivity on the secondary side");
  params.addParam<MaterialPropertyName>("diff_primary", 1, "The diffusivity on the primary side");
  return params;
}

EqualValueConstraint::EqualValueConstraint(const InputParameters & parameters)
  : ADMortarConstraint(parameters),
    _lower_secondary_volume(_assembly.lowerDElemVolume()),
    _lower_primary_volume(_assembly.neighborLowerDElemVolume()),
    _delta(getParam<Real>("delta")),
    _diff_secondary(getADMaterialProperty<Real>("diff_secondary")),
    _diff_primary(getADMaterialProperty<Real>("diff_primary")),
    _stabilize(_delta > TOLERANCE * TOLERANCE)
{
}

ADReal
EqualValueConstraint::computeQpResidual(Moose::MortarType mortar_type)
{
  switch (mortar_type)
  {
    case Moose::MortarType::Secondary:
    {
      // The sign choice here makes it so that for the true solution: lambda = normals_secondary *
      // diff_secondary * grad_u_secondary
      auto residual = -_lambda[_qp] * _test_secondary[_i][_qp];

      if (_stabilize)
        residual += _delta * _lower_secondary_volume *
                    (_diff_secondary[_qp] * _grad_test_secondary[_i][_qp] * _normals[_qp]) *
                    (_lambda[_qp] - _diff_secondary[_qp] * _grad_u_secondary[_qp] * _normals[_qp]);

      return residual;
    }

    case Moose::MortarType::Primary:
    {
      // The sign choice here makes it so that for the true solution: lambda = -normals_primary *
      // diff_primary * grad_u_primary
      auto residual = _lambda[_qp] * _test_primary[_i][_qp];

      if (_stabilize)
        residual +=
            _delta * _lower_primary_volume *
            (_diff_primary[_qp] * _grad_test_primary[_i][_qp] * _normals_primary[_qp]) *
            (-_lambda[_qp] - _diff_primary[_qp] * _grad_u_primary[_qp] * _normals_primary[_qp]);

      return residual;
    }

    case Moose::MortarType::Lower:
    {
      auto residual = (_u_primary[_qp] - _u_secondary[_qp]) * _test[_i][_qp];

      if (_stabilize)
      {
        // secondary
        residual -= _delta * _lower_secondary_volume * _test[_i][_qp] *
                    (_lambda[_qp] - _diff_secondary[_qp] * _grad_u_secondary[_qp] * _normals[_qp]);

        // primary
        residual -=
            _delta * _lower_primary_volume * _test[_i][_qp] *
            (_lambda[_qp] + _diff_primary[_qp] * _grad_u_primary[_qp] * _normals_primary[_qp]);
      }

      return residual;
    }

    default:
      return 0;
  }
}

(moose/framework/src/constraints/PenaltyPeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PenaltyPeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", PenaltyPeriodicSegmentalConstraint);

InputParameters
PenaltyPeriodicSegmentalConstraint::validParams()
{
  InputParameters params = MortarScalarBase::validParams();
  params.addClassDescription(
      "PenaltyPeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using a penalty approach "
      "(no Lagrange multipliers needed). Must be used alongside PenaltyEqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");
  params.addParam<Real>(
      "penalty_value",
      1.0,
      "Penalty value used to impose a generalized force capturing the mortar constraint equation");

  return params;
}

PenaltyPeriodicSegmentalConstraint::PenaltyPeriodicSegmentalConstraint(
    const InputParameters & parameters)
  : DerivativeMaterialInterface<MortarScalarBase>(parameters),
    _temp_jump_global(),
    _tau_s(),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma")),
    _pen_scale(getParam<Real>("penalty_value"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

// Compute the stability parameters to use for all quadrature points
void
PenaltyPeriodicSegmentalConstraint::precalculateResidual()
{
  precalculateStability();
}
void
PenaltyPeriodicSegmentalConstraint::precalculateJacobian()
{
  precalculateStability();
}

// Compute the temperature jump for current quadrature point
void
PenaltyPeriodicSegmentalConstraint::initScalarQpResidual()
{
  precalculateMaterial();
}

Real
PenaltyPeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  // Compute penalty parameter times x-jump times average heat flux
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }
  Real r = _tau_s * (kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Secondary:
      r *= _test_secondary[_i][_qp];
      break;
    case Moose::MortarType::Primary:
      r *= -_test_primary[_i][_qp];
      break;
    case Moose::MortarType::Lower:
      return 0;
    default:
      return 0;
  }
  return r;
}

Real
PenaltyPeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  Real r = _tau_s * _temp_jump_global;
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);

  r *= -dx(_h);

  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_vec(1) = _kappa[1];
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r += dx(_h) * _tau_s * (kappa_vec * dx);
  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

Real
PenaltyPeriodicSegmentalConstraint::computeScalarQpJacobian()
{
  // Stability/penalty term for Jacobian of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);

  Real jac = dx(_h) * _tau_s * dx(_l);

  return jac;
}

Real
PenaltyPeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(
    const Moose::MortarType mortar_type, const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = _tau_s;

  switch (mortar_type)
  {

    case Moose::MortarType::Secondary: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test_secondary[_i][_qp] * dx(_h);
      break;
    case Moose::MortarType::Primary: // Residual_sign -1  ddeltaU_ddisp sign -1;
      jac *= -_test_primary[_i][_qp] * dx(_h);
      break;

    default:
      return 0;
  }
  return jac;
}

Real
PenaltyPeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(
    const Moose::MortarType mortar_type, const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Secondary:
      if (jvar_num != _secondary_var.number())
        return 0;
      break;
    case Moose::MortarType::Primary:
      if (jvar_num != _primary_var.number())
        return 0;
      break;
    default:
      return 0;
  }

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = _tau_s;

  switch (mortar_type)
  {
    case Moose::MortarType::Secondary:
      jac *= (*_phi)[_j][_qp] * dx(_h);
      break;
    case Moose::MortarType::Primary:
      jac *= -(*_phi)[_j][_qp] * dx(_h);
      break;

    default:
      return 0;
  }
  return jac;
}

void
PenaltyPeriodicSegmentalConstraint::precalculateStability()
{
  // Example showing how the penalty could be loaded from some function
  _tau_s = _pen_scale;
}

// Compute temperature jump and flux average/jump
void
PenaltyPeriodicSegmentalConstraint::precalculateMaterial()
{
  _temp_jump_global = (_u_primary[_qp] - _u_secondary[_qp]);
}

(moose/framework/src/constraints/PeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", PeriodicSegmentalConstraint);

InputParameters
PeriodicSegmentalConstraint::validParams()
{
  InputParameters params = MortarScalarBase::validParams();
  params.addClassDescription(
      "PeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

PeriodicSegmentalConstraint::PeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<MortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

Real
PeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }
  Real r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

Real
PeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

Real
PeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(const Moose::MortarType mortar_type,
                                                            const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

Real
PeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(const Moose::MortarType mortar_type,
                                                            const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      if (!_var || _var->number() != jvar_num)
        return 0;
      break;
    default:
      return 0;
  }

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      jac *= (*_phi)[_j][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

(moose/framework/src/constraints/PeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", PeriodicSegmentalConstraint);

InputParameters
PeriodicSegmentalConstraint::validParams()
{
  InputParameters params = MortarScalarBase::validParams();
  params.addClassDescription(
      "PeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

PeriodicSegmentalConstraint::PeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<MortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

Real
PeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  RealVectorValue kappa_vec(_kappa[0], 0, 0);
  if (_k_order == 2)
    kappa_vec(1) = _kappa[1];
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
  }
  Real r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

Real
PeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

Real
PeriodicSegmentalConstraint::computeQpOffDiagJacobianScalar(const Moose::MortarType mortar_type,
                                                            const unsigned int svar_num)
{
  if (svar_num != _kappa_var)
    return 0;

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower: // Residual_sign -1  ddeltaU_ddisp sign 1;
      jac *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

Real
PeriodicSegmentalConstraint::computeScalarQpOffDiagJacobian(const Moose::MortarType mortar_type,
                                                            const unsigned int jvar_num)
{
  // Test if jvar is the ID of the primary variables and not some other random variable
  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      if (!_var || _var->number() != jvar_num)
        return 0;
      break;
    default:
      return 0;
  }

  // Stability/penalty term for Jacobian
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  Real jac = -dx(_h);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      jac *= (*_phi)[_j][_qp];
      break;
    default:
      return 0;
  }
  return jac;
}

(moose/framework/src/constraints/ADPeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ADPeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", ADPeriodicSegmentalConstraint);

InputParameters
ADPeriodicSegmentalConstraint::validParams()
{
  InputParameters params = ADMortarScalarBase::validParams();
  params.addClassDescription(
      "ADPeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

ADPeriodicSegmentalConstraint::ADPeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<ADMortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

ADReal
ADPeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADRealVectorValue kappa_vec(_kappa[0], 0, 0);
  Moose::derivInsert(kappa_vec(0).derivatives(), _kappa_var_ptr->dofIndices()[0], 1);
  if (_k_order == 2)
  {
    kappa_vec(1) = _kappa[1];
    Moose::derivInsert(kappa_vec(1).derivatives(), _kappa_var_ptr->dofIndices()[1], 1);
  }
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
    Moose::derivInsert(kappa_vec(1).derivatives(), _kappa_var_ptr->dofIndices()[1], 1);
    Moose::derivInsert(kappa_vec(2).derivatives(), _kappa_var_ptr->dofIndices()[2], 1);
  }
  ADReal r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

ADReal
ADPeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADReal r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

(moose/framework/src/constraints/ADPeriodicSegmentalConstraint.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ADPeriodicSegmentalConstraint.h"

registerMooseObject("MooseApp", ADPeriodicSegmentalConstraint);

InputParameters
ADPeriodicSegmentalConstraint::validParams()
{
  InputParameters params = ADMortarScalarBase::validParams();
  params.addClassDescription(
      "ADPeriodicSegmentalConstraint enforces macro-micro periodic conditions between "
      "secondary and primary sides of a mortar interface using Lagrange multipliers."
      "Must be used alongside EqualValueConstraint.");
  params.renameCoupledVar("scalar_variable", "epsilon", "Primary coupled scalar variable");
  params.addRequiredCoupledVar("sigma", "Controlled scalar averaging variable");

  return params;
}

ADPeriodicSegmentalConstraint::ADPeriodicSegmentalConstraint(const InputParameters & parameters)
  : DerivativeMaterialInterface<ADMortarScalarBase>(parameters),
    _kappa_aux_ptr(getScalarVar("sigma", 0)),
    _ka_order(_kappa_aux_ptr->order()),
    _kappa_aux(coupledScalarValue("sigma"))
{
  if (_kappa_aux_ptr->kind() != Moose::VarKindType::VAR_AUXILIARY)
    paramError("sigma",
               "Must assign auxiliary scalar variable to sigma, rather than nonlinear variable");
}

ADReal
ADPeriodicSegmentalConstraint::computeQpResidual(const Moose::MortarType mortar_type)
{
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADRealVectorValue kappa_vec(_kappa[0], 0, 0);
  Moose::derivInsert(kappa_vec(0).derivatives(), _kappa_var_ptr->dofIndices()[0], 1);
  if (_k_order == 2)
  {
    kappa_vec(1) = _kappa[1];
    Moose::derivInsert(kappa_vec(1).derivatives(), _kappa_var_ptr->dofIndices()[1], 1);
  }
  else if (_k_order == 3)
  {
    kappa_vec(1) = _kappa[1];
    kappa_vec(2) = _kappa[2];
    Moose::derivInsert(kappa_vec(1).derivatives(), _kappa_var_ptr->dofIndices()[1], 1);
    Moose::derivInsert(kappa_vec(2).derivatives(), _kappa_var_ptr->dofIndices()[2], 1);
  }
  ADReal r = -(kappa_vec * dx);

  switch (mortar_type)
  {
    case Moose::MortarType::Lower:
      r *= _test[_i][_qp];
      break;
    default:
      return 0;
  }
  return r;
}

ADReal
ADPeriodicSegmentalConstraint::computeScalarQpResidual()
{
  // Stability/penalty term for residual of scalar variable
  RealVectorValue dx(_phys_points_primary[_qp] - _phys_points_secondary[_qp]);
  ADReal r = -dx(_h) * _lambda[_qp];

  RealVectorValue kappa_aux_vec(_kappa_aux[0], 0, 0);
  if (_k_order == 2)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
  }
  else if (_k_order == 3)
  {
    kappa_aux_vec(1) = _kappa_aux[1];
    kappa_aux_vec(2) = _kappa_aux[2];
  }

  r -= dx(_h) * (kappa_aux_vec * _normals[_qp]);

  return r;
}

(moose/test/include/constraints/TestPeriodicSole.h)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "MortarScalarBase.h"
#include "DerivativeMaterialInterface.h"
#include "MooseVariableScalar.h"
#include "Assembly.h"

/// Test object to illustrate coupling between mortar spatial variables
/// and two separate scalar variables, kappa and kappa_other. The solution
/// is equivalent to PenaltyPeriodicSegmentalConstraint.
///
/// The kappa variable is split into two scalars: the first component called '_kappa'
/// herein and all other components called '_kappa_other' herein.
/// This decomposition is ONLY valid for 2D problems with a macro scalar vector [kappa_x kappa_y]
/// applied as a periodic boundary condition. For parameter _alpha = 0, the primary
/// scalar (_kappa) is kappa_x and the coupled scalar is kappa_y. For parameter _alpha = 1, the primary
/// scalar (_kappa) is kappa_y and the coupled scalar is kappa_x.
/// That mapping MUST be handled in the input file.
///
/// Thus, each instance of TestPeriodicSole acts on both primary/secondary field variables
/// and one scalar variable (_svar_alpha). The job of the constraint is to assemble the
/// residual and Jacobian terms arising from _svar_alpha (namely, selected entries).
/// It assembles one row of the Jacobian for the scalar variable as well as the couplings of that
/// scalar to the primary and secondary spatial variables, using logical checks with _alpha.
/// The entries for the other field/scalar variables are handled by the other instance of the
/// constraint, which have the opposite value of _alpha. The logical checks ensure the proper
/// decomposition of the jobs.
///
/// In summary, for x=disp_x etc. and k=kappa_x and a=kappa_y, then the contributions of the instances are
/// PenaltyEqualValueConstraint
/// R = [Rss+p, Rps+p, 0, 0]^T
/// J = [Jss, Jsp, 000, 000
///      Jps, Jpp, 000, 000
///      000, 000, 000, 000
///      000, 000, 000, 000]
/// TestPeriodicSole
/// _alpha=0
/// R = [Rsk, Rpk,  Rk,  00]^T
/// J = [000, 000, Jsk, 000
///      000, 000, Jpk, 000
///      Jks, Jkp, Jkk, Jka
///      000, 000, 000, 000]
/// _alpha=1
/// R = [Rsa, Rpa,  00,  Ra]^T
/// J = [000, 000, 000, Jsa
///      000, 000, 000, Jpa
///      000, 000, 000, 000
///      Jas, Jap, Jak, Jaa]
///
/// In this manner, the full R and J are obtained with NO duplication of jobs:
/// R = [Rs,  Rp,  Rk,  Ra ]^T
/// J = [Jss, Jsp, Jsk, Jsa
///      Jps, Jpp, Jpk, Jpa
///      Jks, Jkp, Jkk, Jka
///      Jas, Jap, Jak, Jaa]
///

class TestPeriodicSole : public DerivativeMaterialInterface<MortarScalarBase>
{
public:
  static InputParameters validParams();
  TestPeriodicSole(const InputParameters & parameters);

protected:
  virtual void precalculateResidual() override;
  virtual void precalculateJacobian() override;
  virtual void initScalarQpResidual() override;

  /**
   * Method for computing the residual at quadrature points
   */
  virtual Real computeQpResidual(Moose::MortarType mortar_type) override;
  Real computeQpJacobian(Moose::ConstraintJacobianType /*jacobian_type*/,
                         unsigned int /*jvar*/) override
  {
    return 0;
  };

  /**
   * Method for computing the scalar part of residual at quadrature points
   */
  virtual Real computeScalarQpResidual() override;

  /**
   * Method for computing the scalar variable part of Jacobian at
   * quadrature points
   */
  virtual Real computeScalarQpJacobian() override;

  // using MortarScalarBase::computeOffDiagJacobianScalar;

  /**
   * Method for computing d-_var-residual / d-_kappa at quadrature points.
   */
  virtual Real computeQpOffDiagJacobianScalar(Moose::MortarType mortar_type,
                                              const unsigned int svar_num) override;

  /**
   * Method for computing d-_kappa-residual / d-_var at quadrature points.
   */
  virtual Real computeScalarQpOffDiagJacobian(Moose::MortarType mortar_type,
                                              const unsigned int jvar_num) override;

  /**
   * Method for computing an off-diagonal jacobian component at quadrature points.
   */
  virtual Real computeScalarQpOffDiagJacobianScalar(const unsigned int svar_num) override;

  // Compute T jump and heat flux average/jump
  void precalculateMaterial();
  // Compute four stability tensors
  void precalculateStability();

protected:
  /// the temperature jump in global and interface coordiantes;
  /// TM-analogy: _displacement_jump_global, _interface_displacement_jump
  ///@{
  Real _temp_jump_global;
  ///@}

  /// The four stability parameters from the VMDG method
  ///@{
  Real _tau_s;
  ///@}

  /// The controlled scalar variable ID
  const unsigned int _kappa_aux_var;

  /// Order of the homogenization variable, used in several places
  const unsigned int _ka_order;

  /// The controlled scalar variable
  const VariableValue & _kappa_aux;

  /// Which component of the vector residual this constraint is responsible for
  const unsigned int _alpha;

  /// (Pointer to) Scalar variable this constraint operates on
  const MooseVariableScalar * const _kappao_var_ptr;

  /// The unknown scalar variable ID
  const unsigned int _kappao_var;

  /// Order of the scalar variable, used in several places
  const unsigned int _ko_order;

  /// Reference to the current solution at the current quadrature point
  const VariableValue & _kappa_other;

  const Real & _current_elem_volume;

  const Real & _current_side_volume;

  /// Input property to allow user modifying penalty parameter
  const Real _pen_scale;

private:
  /// hard code the penalty for now
  const Real pencoef = 1.0;
};

(moose/test/tests/mortar/periodic_segmental_constraint/testperiodicsole.i)

[Mesh]
  [left_block]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -1.0
    xmax = 1.0
    ymin = -1.0
    ymax = 1.0
    nx = 2
    ny = 2
    elem_type = QUAD4
  []
  [left_block_sidesets]
    type = RenameBoundaryGenerator
    input = left_block
    old_boundary = '0 1 2 3'
    new_boundary = '10 11 12 13'
  []
  [left_block_id]
    type = SubdomainIDGenerator
    input = left_block_sidesets
    subdomain_id = 1
  []
  [left]
    type = LowerDBlockFromSidesetGenerator
    input = left_block_id
    sidesets = '13'
    new_block_id = '10003'
    new_block_name = 'secondary_left'
  []
  [right]
    type = LowerDBlockFromSidesetGenerator
    input = left
    sidesets = '11'
    new_block_id = '10001'
    new_block_name = 'primary_right'
  []
  [bottom]
    type = LowerDBlockFromSidesetGenerator
    input = right
    sidesets = '10'
    new_block_id = '10000'
    new_block_name = 'secondary_bottom'
  []
  [top]
    type = LowerDBlockFromSidesetGenerator
    input = bottom
    sidesets = '12'
    new_block_id = '10002'
    new_block_name = 'primary_top'
  []

  [corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = top
  []
[]

[Variables]
  [u]
    order = FIRST
    family = LAGRANGE
  []
  [kappa_x]
    order = FIRST
    family = SCALAR
  []
  [kappa_y]
    order = FIRST
    family = SCALAR
  []
[]

[AuxVariables]
  [kappa_aux]
    order = SECOND
    family = SCALAR
  []
  [./flux_x]
      order = FIRST
      family = MONOMIAL
  [../]
  [./flux_y]
      order = FIRST
      family = MONOMIAL
  [../]
[]

[AuxScalarKernels]
  [kappa]
    type = FunctionScalarAux
    variable = kappa_aux
    function = '1 3'
    execute_on = initial #timestep_end
  []
[]

[AuxKernels]
  [./flux_x]
    type = DiffusionFluxAux
    diffusivity = 'conductivity'
    variable = flux_x
    diffusion_variable = u
    component = x
    block = 1
  [../]
  [./flux_y]
    type = DiffusionFluxAux
    diffusivity = 'conductivity'
    variable = flux_y
    diffusion_variable = u
    component = y
    block = 1
  [../]
[]

[Kernels]
  [diff1]
    type = Diffusion
    variable = u
    block = 1
  []
[]

[Materials]
  [k1]
    type = GenericConstantMaterial
    prop_names = 'conductivity'
    prop_values = 1.0
    block = 1
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [fix_right]
    type = DirichletBC
    variable = u
    boundary = pinned_node
    value = 0
  []
[]

[Constraints]
  [mortarlr]
    type = PenaltyEqualValueConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodiclrx]
    type = TestPeriodicSole
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    kappa = kappa_x
    kappa_aux = kappa_aux
    component = 0
    kappa_other = kappa_y
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodiclry]
    type = TestPeriodicSole
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    kappa = kappa_y
    kappa_aux = kappa_aux
    component = 1
    kappa_other = kappa_x
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [mortarbt]
    type = PenaltyEqualValueConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodicbtx]
    type = TestPeriodicSole
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    kappa = kappa_x
    kappa_aux = kappa_aux
    component = 0
    kappa_other = kappa_y
    correct_edge_dropping = true
    penalty_value = 1.e3
  []
  [periodicbty]
    type = TestPeriodicSole
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    kappa = kappa_y
    kappa_aux = kappa_aux
    component = 1
    kappa_other = kappa_x
    correct_edge_dropping = true
    compute_scalar_residuals = true
    penalty_value = 1.e3
  []
[]

[Preconditioning]
  [smp]
    full = true
    type = SMP
  []
[]

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  solve_type = NEWTON
[]

[Postprocessors]
  [max]
    type = ElementExtremeValue
    variable = 'flux_x'
  []
[]

[Outputs]
  csv = true
[]

(moose/test/tests/mortar/periodic_segmental_constraint/periodic_simple2d.i)

[Mesh]
  [left_block]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -1.0
    xmax = 1.0
    ymin = -1.0
    ymax = 1.0
    nx = 2
    ny = 2
    elem_type = QUAD9
  []
  [left_block_sidesets]
    type = RenameBoundaryGenerator
    input = left_block
    old_boundary = '0 1 2 3'
    new_boundary = '10 11 12 13'
  []
  [left_block_id]
    type = SubdomainIDGenerator
    input = left_block_sidesets
    subdomain_id = 1
  []
  [left]
    type = LowerDBlockFromSidesetGenerator
    input = left_block_id
    sidesets = '13'
    new_block_id = '10003'
    new_block_name = 'secondary_left'
  []
  [right]
    type = LowerDBlockFromSidesetGenerator
    input = left
    sidesets = '11'
    new_block_id = '10001'
    new_block_name = 'primary_right'
  []
  [bottom]
    type = LowerDBlockFromSidesetGenerator
    input = right
    sidesets = '10'
    new_block_id = '10000'
    new_block_name = 'secondary_bottom'
  []
  [top]
    type = LowerDBlockFromSidesetGenerator
    input = bottom
    sidesets = '12'
    new_block_id = '10002'
    new_block_name = 'primary_top'
  []

  [corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = top
  []
[]

[Variables]
  [u]
    order = SECOND
    family = LAGRANGE
  []
  [epsilon]
    order = SECOND
    family = SCALAR
  []
  [./lm1]
    order = FIRST
    family = LAGRANGE
    block = secondary_left
  [../]
  [./lm2]
    order = FIRST
    family = LAGRANGE
    block = secondary_bottom
  [../]
[]

[AuxVariables]
  [sigma]
    order = SECOND
    family = SCALAR
  []
[]

[AuxScalarKernels]
  [sigma]
    type = FunctionScalarAux
    variable = sigma
    function = '1 2'
    execute_on = initial #timestep_end
  []
[]

[Kernels]
  [diff1]
    type = Diffusion
    variable = u
    block = 1
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [fix_right]
    type = DirichletBC
    variable = u
    boundary = pinned_node
    value = 0
  []
[]

[Constraints]
  [mortarlr]
    type = EqualValueConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    variable = lm1
    correct_edge_dropping = true
  []
  [periodiclr]
    type = PeriodicSegmentalConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm1
    correct_edge_dropping = true
  []
  [mortarbt]
    type = EqualValueConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    variable = lm2
    correct_edge_dropping = true
  []
  [periodicbt]
    type = PeriodicSegmentalConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm2
    correct_edge_dropping = true
  []
[]

[Preconditioning]
  [smp]
    full = true
    type = SMP
  []
[]

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15'
  solve_type = NEWTON
[]

[Outputs]
  exodus = true
  csv = true
[]

(moose/test/tests/mortar/periodic_segmental_constraint/periodic_simple2d.i)

[Mesh]
  [left_block]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -1.0
    xmax = 1.0
    ymin = -1.0
    ymax = 1.0
    nx = 2
    ny = 2
    elem_type = QUAD9
  []
  [left_block_sidesets]
    type = RenameBoundaryGenerator
    input = left_block
    old_boundary = '0 1 2 3'
    new_boundary = '10 11 12 13'
  []
  [left_block_id]
    type = SubdomainIDGenerator
    input = left_block_sidesets
    subdomain_id = 1
  []
  [left]
    type = LowerDBlockFromSidesetGenerator
    input = left_block_id
    sidesets = '13'
    new_block_id = '10003'
    new_block_name = 'secondary_left'
  []
  [right]
    type = LowerDBlockFromSidesetGenerator
    input = left
    sidesets = '11'
    new_block_id = '10001'
    new_block_name = 'primary_right'
  []
  [bottom]
    type = LowerDBlockFromSidesetGenerator
    input = right
    sidesets = '10'
    new_block_id = '10000'
    new_block_name = 'secondary_bottom'
  []
  [top]
    type = LowerDBlockFromSidesetGenerator
    input = bottom
    sidesets = '12'
    new_block_id = '10002'
    new_block_name = 'primary_top'
  []

  [corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = top
  []
[]

[Variables]
  [u]
    order = SECOND
    family = LAGRANGE
  []
  [epsilon]
    order = SECOND
    family = SCALAR
  []
  [./lm1]
    order = FIRST
    family = LAGRANGE
    block = secondary_left
  [../]
  [./lm2]
    order = FIRST
    family = LAGRANGE
    block = secondary_bottom
  [../]
[]

[AuxVariables]
  [sigma]
    order = SECOND
    family = SCALAR
  []
[]

[AuxScalarKernels]
  [sigma]
    type = FunctionScalarAux
    variable = sigma
    function = '1 2'
    execute_on = initial #timestep_end
  []
[]

[Kernels]
  [diff1]
    type = Diffusion
    variable = u
    block = 1
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [fix_right]
    type = DirichletBC
    variable = u
    boundary = pinned_node
    value = 0
  []
[]

[Constraints]
  [mortarlr]
    type = EqualValueConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    variable = lm1
    correct_edge_dropping = true
  []
  [periodiclr]
    type = PeriodicSegmentalConstraint
    primary_boundary = '11'
    secondary_boundary = '13'
    primary_subdomain = 'primary_right'
    secondary_subdomain = 'secondary_left'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm1
    correct_edge_dropping = true
  []
  [mortarbt]
    type = EqualValueConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    variable = lm2
    correct_edge_dropping = true
  []
  [periodicbt]
    type = PeriodicSegmentalConstraint
    primary_boundary = '12'
    secondary_boundary = '10'
    primary_subdomain = 'primary_top'
    secondary_subdomain = 'secondary_bottom'
    secondary_variable = u
    epsilon = epsilon
    sigma = sigma
    variable = lm2
    correct_edge_dropping = true
  []
[]

[Preconditioning]
  [smp]
    full = true
    type = SMP
  []
[]

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       NONZERO               1e-15'
  solve_type = NEWTON
[]

[Outputs]
  exodus = true
  csv = true
[]