KernelScalarBase

The KernelScalarBase 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 Kernel 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 interior of a domain. As usual, this piece of physics is referred to as the "residual" and is evaluated at integration quadrature points within that domain. To implement your own physics in MOOSE, you create your own kernel by subclassing the MOOSE KernelScalarBase class.

The Kernel scalar augmentation system supports the use of AD for residual calculations, as such there are two options for creating field-scalar coupling objects: KernelScalarBase and ADKernelScalarBase. 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 Kernel Scalar Coupling Classes

Each Kernel object has a focus field variable or spatial variable; its job is to contribute to the residual as well as the row of the Jacobian matrix. Herein, as in the source code of Kernels, this spatial variable will be called _var. In a coupled (multi-phyics) weak form, all domain integral terms containing the test function of _var are potential candidates for Kernel contributions.

The philosphy of the scalar augmentation class KernelScalarBase is to add a focus scalar variable referred to as _kappa to the Kernel object so that all terms in the coupled weak form that involve _var, _kappa, 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 KernelScalarBase 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() function must be overridden as usual for Kernels, 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 component d-_var-residual / d-_var
computeQpOffDiagJacobian(jvar_num): off-diagonal Jacobian component d-_var-residual / d-jvar
computeQpOffDiagJacobianScalar(svar_num): off-diagonal Jacobian component d-_var-residual / d-svar
computeScalarQpJacobian(): Jacobian component d-_kappa-residual / d-_kappa
computeScalarQpOffDiagJacobian(jvar_num): off-diagonal Jacobian component d-_kappa-residual / d-jvar
computeScalarQpOffDiagJacobianScalar(svar_num): off-diagonal Jacobian component d-_kappa-residual / d-svar

Examples of some of these methods are shown below in Examples from Source Code. 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.

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 Kernel documentation, you have access to several member variables inside your KernelScalarBase 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 Kernel 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

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

note:Parallelization of scalar contributions

While these quadrature loops are convenient for implementation in a single object, the speed of parallel execution may be slower due to the sequential assembly needed from each element assemblying to the same scalar variable _kappa. For greater speed, the developer may instead implement the terms for computeScalarQpResidual() and computeScalarQpJacobian() through a derived class of ElementIntegralUserObject as discussed at Coupling with Spatial Variables.

Examples from Source Code

As mentioned, the computeScalarQpResidual method should be overridden for both flavors of kernels, non-AD and AD. As an example, consider the scalar residual weak form term of the ScalarLMKernel class:

$(1)var element = document.getElementById("moose-equation-c2cba404-b635-4713-9c8a-b9edee34a4a3");katex.render(" F^{(\\lambda)} \\equiv \\int_{\\Omega} \\phi^h \\;\\text{d}\\Omega - V_0 = 0 ", element, {displayMode:true,throwOnError:false});$

The ScalarLMKernel class is implemented using the GenericKernelScalar template class to contain both the AD and non-AD versions within the same source files; the test sources files in the Tensor Mechanics module described at the bottom of this section appear more simply since they are non-AD only: Listing 5. The computeScalarQpResidual method for this class is provided in Listing 1, where _value/_pp_value is equal to $var element = document.getElementById("moose-equation-2c4bc6e4-670d-40ab-86cc-47f629df982a");katex.render("V_0", element, {displayMode:false,throwOnError:false});$ .

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

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

Meanwhile, the contribution to the 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-d186cf98-7e8a-4389-8db3-7d7df8270a43");katex.render("\\lambda^h", element, {displayMode:false,throwOnError:false});$ in this weak form).

$(2)var element = document.getElementById("moose-equation-8ecfc68c-ecb6-4570-8de9-b695c41eb2a4");katex.render(" F^{(\\phi)}_i \\equiv \\lambda^h \\int_{\\Omega} \\varphi_i \\;\\text{d}\\Omega ", element, {displayMode:true,throwOnError:false});$

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

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

This object also overrides the computeScalarQpOffDiagJacobian method to define the Jacobian term related to Eq. (1) as shown in Listing 3.

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

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

Notice that there is a conditional to confirm that the coupled jvar is the focus variable _var, otherwise it returns zero. Also, this method only returns a "Real" value since this method is only called by the non-AD version of the class during Jacobian computation; an assert is used to verify this intention.

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

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

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

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

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

The AD version of this object, ADScalarLMKernel, only requires the residual implementation. A solely AD source file would only need to override computeScalarQpResidual and computeQpResidual and leave all the Jacobian methods as base definitions, which return zero. See MortarScalarBase for examples of AD-only and non-AD separate classes.

As a more complicated example of the scalar augmentation system for kernels, the Tensor Mechanics test app contains headers, source, and test files for an alternative implementation of the "HomogenizedTotalLagrangianStressDivergence" system from the Tensor Mechanics module. This Kernel is designated with the suffix "S" to distinguish from the existing objects in the module. Also, there are other intermediate classes such as "TotalLagrangianStressDivergence" that are also designated with an "S" suffix. These other classes are needed since the lower class needs to also derive from KernelScalarBase. Meanwhile, they do not need the scalar_variable parameter and function identically to their original module source object; see the Parameters section for a comment about leaving this parameter blank.

The scalar augmentation system is designed such that multiple scalar variables can be coupled to an instance of the Kernel 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.

Examples of two schemes for decomposing the coupling terms and having multiple scalar variables are contained in the source files of the Tensor Mechanics module test directory as well as input files 2drow.i and 2dsole.i, with listings below. The comments within these header and source files serve as documentation and should be consulted to visualize how the rows and columns of the relevant residual and Jacobian contributions are handled. The suffix "R" refers to assembling the entire row of the Jacobian in one object, and the suffix "A" refers to assembling symmetric pairs of the residual and Jacobian; see the header file for clarification.

Listing 5: Organization of spatial and scalar variable contributions by row.

/// Total Lagrangian formulation with most homogenization terms (one disp_xyz field and one scalar)
/// The macro_gradient variable is split into two scalars: the first component called '_hvar'
/// herein and all other components called '_avar' herein. For parameter _beta = 0, the primary
/// scalar (_kappa) is _hvar and the coupled scalar is _avar. For parameter _beta = 1, the primary
/// scalar (_kappa) is _avar and the coupled scalar is _hvar. Just like the primary field variable
/// (_var) is either disp_x or disp_y or disp_z depending on _alpha.
///
/// Thus, each instance of HomogenizedTotalLagrangianStressDivergenceR acts on one field variable
/// (_disp_alpha) and one scalar variable (_hvar_beta). The job of the kernel is to assemble the
/// residual of all dofs of _disp_alpha and of all dofs of _hvar_beta (namely, selected rows).
/// Also, it assembles the ENTIRE row for _disp_alpha and _hvar_beta (namely the columns
/// from all dofs of all _disp field variables and all dofs of all scalar variables _hvar and
/// _avar). The rows for the other field/scalar variables are handled by other instances of the
/// kernel, according to the flags compute_scalar_residuals and compute_field_residuals.
/// When compute_field_residuals is given, only component=_alpha matters and beta = {0,1}

[Kernels]
  [sdx]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = true
    compute_scalar_residuals = false
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_y
    component = 1
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = true
    compute_scalar_residuals = false
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sd0]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = false
    compute_scalar_residuals = true
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sd1]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_y
    component = 1
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    compute_field_residuals = false
    compute_scalar_residuals = true
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

Listing 6: Organization of spatial and scalar variable contributions by symmetric pairs.

/// Total Lagrangian formulation with most homogenization terms (one disp_xyz field and one scalar)
/// The macro_gradient variable is split into two scalars: the first component called '_hvar'
/// herein and all other components called '_avar' herein. For parameter _beta = 0, the primary
/// scalar (_kappa) is _hvar and the coupled scalar is _avar. For parameter _beta = 1, the primary
/// scalar (_kappa) is _avar and the coupled scalar is _hvar. Just like the primary field variable
/// (_var) is either disp_x or disp_y or disp_z depending on _alpha.
///
/// Thus, each instance of HomogenizedTotalLagrangianStressDivergenceA acts on one field variable
/// (_disp_alpha) and one scalar variable (_hvar_beta). The job of the kernel is to assemble the
/// residual of all dofs of _disp_alpha and of all dofs of _hvar_beta (namely, selected entries).
/// It assembles a symmetric portion of the Jacobian for _disp_alpha and _hvar_beta, with some
/// logical checks to access only the particular desired terms (often for _alpha or _beta = 0).
/// The entries for the other field/scalar variables are handled by other instances of the
/// kernel, which have other values of _alpha AND _beta. The logical checks ensure the proper
/// decomposition of the jobs.
///
/// In summary, for x=disp_x etc. and h=_hvar and a=_avar, then the contributions of the instances are
/// _alpha=0, _beta=0
/// R = [Rx,  00,  00,  Rh,  00 ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, 000
///      Jhx, 000, 000, Jhh, Jha]
/// _alpha=1, _beta=0
/// R = [00,  Ry,  00,  00,  00 ]^T
/// J = [Jyx, Jyy, Jyz, Jyh, 000
///      000, Jhy, 000, 000, 000]
/// _alpha=2, _beta=0
/// R = [00,  00,  Rz,  00,  00 ]^T
/// J = [Jzx, Jzy, Jzz, Jzh, 000
///      000, 000, Jhz, 000, 000]
/// _alpha=0, _beta=1
/// R = [00,  00,  00,  00,  Ra ]^T
/// J = [000, 000, 000, 000, Jxa
///      Jax, 000, 000, Jah, Jaa]
/// _alpha=1, _beta=1
/// R = [00,  00,  00,  00,  00 ]^T
/// J = [000, 000, 000, 000, Jya
///      000, Jay, 000, 000, 000]
/// _alpha=2, _beta=1
/// R = [00,  00,  00,  00,  00 ]^T
/// J = [000, 000, 000, 000, Jza
///      000, 000, Jaz, 000, 000]
///
/// In this manner, the full R and J are obtained with NO duplication of jobs:
/// R = [Rx,  Ry,  Rz,  Rh,  Ra ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, Jxa
///      Jyx, Jyy, Jyz, Jyh, Jya
///      Jzx, Jzy, Jzz, Jzh, Jza
///      Jhx, Jhy, Jhz, Jhh, Jha
///      Jax, Jay, Jaz, Jah, Jaa]
///
class HomogenizedTotalLagrangianStressDivergenceA : public TotalLagrangianStressDivergenceS
{
public:
  static InputParameters validParams();
  HomogenizedTotalLagrangianStressDivergenceA(const InputParameters & parameters);

protected:
  // Add overrides to base class contributions to only happen for _beta==0, to happen only once
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobianDisplacement(unsigned int alpha, unsigned int beta) override;

  /**
   * Method for computing the scalar part of residual for _kappa
   */
  virtual void computeScalarResidual() override;

  /**
   * Method for computing the scalar variable part of Jacobian for d-_kappa-residual / d-_kappa
   */
  virtual void computeScalarJacobian() override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-jvar
   * jvar is looped over all field variables, which herein is just disp_x and disp_y
   */
  virtual void computeScalarOffDiagJacobian(const unsigned int jvar_num) override;

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

  /**
   * Method for computing an off-diagonal jacobian component d-_var-residual / d-svar.
   * svar is looped over all scalar variables, which herein is just _kappa and _kappa_other
   */
  virtual void computeOffDiagJacobianScalarLocal(const unsigned int svar_num) override;

  /**
   * Method for computing d-_var-residual / d-_svar at quadrature points.
   */
  virtual Real computeQpOffDiagJacobianScalar(const unsigned int svar_num) override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-svar
   * svar is looped over other scalar variables, which herein is just _kappa_other
   */
  virtual void computeScalarOffDiagJacobianScalar(const unsigned int svar_num) override;

protected:
  /// Which component of the scalar vector residual this constraint is responsible for
  const unsigned int _beta;

  /// (Pointer to) Scalar variable this kernel 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;

  /// Type of each constraint (stress or strain) for each component
  HomogenizationA::ConstraintMap _cmap;

  /// The constraint type; initialize with 'none'
  HomogenizationA::ConstraintType _ctype = HomogenizationA::ConstraintType::None;

  /// Used internally to iterate over each scalar component
  unsigned int _m;
  unsigned int _n;
  unsigned int _a;
  unsigned int _b;
}

[Kernels]
  [sdx0]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy0]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_y
    component = 1
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdx1]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_x
    component = 0
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy1]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_y
    component = 1
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

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 kernel derived from KernelScalarBase:

scalar_variable: the focus scalar variable of the kernel, 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 Kernel, 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 ScalarLMKernel object with the scalar_variable parameter renamed to kappa:

[Kernels]
  [diff]
    type = Diffusion
    variable = u
  []

  [ffnk]
    type = BodyForce
    variable = u
    function = ffn
  []

  [sk_lm]
    type = ScalarLMKernel
    variable = u
    kappa = lambda
    pp_name = pp
    value = 2.666666666666666
  []
[]

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 are also some optional parameters that can be supplied to KernelScalarBase classes. They are:

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.
compute_field_residuals: Whether to compute residuals for the primal field variable. If several KernelScalarBase objects are used in the input file to compute different rows (i.e. different variables) of the global residual, then some objects can be targeted to field variable rows and others to scalar variable rows.

(moose/framework/src/kernels/ScalarLMKernel.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 "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(moose/framework/src/kernels/ScalarLMKernel.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 "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(moose/framework/src/kernels/ScalarLMKernel.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 "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(moose/framework/src/kernels/ScalarLMKernel.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 "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(moose/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceR.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 "TotalLagrangianStressDivergenceS.h"

// Helpers common to the whole homogenization system
namespace HomogenizationR
{
/// Moose constraint type, for input
const MultiMooseEnum constraintType("strain stress none");
/// Constraint type: stress/PK stress or strain/deformation gradient
enum class ConstraintType
{
  Strain,
  Stress,
  None
};
typedef std::map<std::pair<unsigned int, unsigned int>, std::pair<ConstraintType, const Function *>>
    ConstraintMap;
}

/// Total Lagrangian formulation with most homogenization terms (one disp_xyz field and one scalar)
/// The macro_gradient variable is split into two scalars: the first component called '_hvar'
/// herein and all other components called '_avar' herein. For parameter _beta = 0, the primary
/// scalar (_kappa) is _hvar and the coupled scalar is _avar. For parameter _beta = 1, the primary
/// scalar (_kappa) is _avar and the coupled scalar is _hvar. Just like the primary field variable
/// (_var) is either disp_x or disp_y or disp_z depending on _alpha.
///
/// Thus, each instance of HomogenizedTotalLagrangianStressDivergenceR acts on one field variable
/// (_disp_alpha) and one scalar variable (_hvar_beta). The job of the kernel is to assemble the
/// residual of all dofs of _disp_alpha and of all dofs of _hvar_beta (namely, selected rows).
/// Also, it assembles the ENTIRE row for _disp_alpha and _hvar_beta (namely the columns
/// from all dofs of all _disp field variables and all dofs of all scalar variables _hvar and
/// _avar). The rows for the other field/scalar variables are handled by other instances of the
/// kernel, according to the flags compute_scalar_residuals and compute_field_residuals.
/// When compute_field_residuals is given, only component=_alpha matters and beta = {0,1} is looped.
/// When compute_scalar_residuals is given, only prime_scalar=_beta matters and alpha = {0,1,2} is looped.
///
/// In summary, for x=disp_x etc. and h=_hvar and a=_avar, then the contributions of the instances are
/// _alpha=0
/// R = [Rx,  00,  00,  00,  00 ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, Jxa]
/// _alpha=1
/// R = [00,  Ry,  00,  00,  00 ]^T
/// J = [Jyx, Jyy, Jyz, Jyh, Jya]
/// _alpha=2
/// R = [00,  00,  Rz,  00,  00 ]^T
/// J = [Jzx, Jzy, Jzz, Jzh, Jza]
/// _beta=0
/// R = [00,  00,  00,  Rh,  00 ]^T
/// J = [Jhx, Jhy, Jhz, Jhh, Jha]
/// _beta=1
/// R = [00,  00,  00,  00,  Ra ]^T
/// J = [Jax, Jay, Jaz, Jah, Jaa]
///
/// In this manner, the full R and J are obtained with NO duplication of jobs:
/// R = [Rx,  Ry,  Rz,  Rh,  Ra ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, Jxa
///      Jyx, Jyy, Jyz, Jyh, Jya
///      Jzx, Jzy, Jzz, Jzh, Jza
///      Jhx, Jhy, Jhz, Jhh, Jha
///      Jax, Jay, Jaz, Jah, Jaa]
///
class HomogenizedTotalLagrangianStressDivergenceR : public TotalLagrangianStressDivergenceS
{
public:
  static InputParameters validParams();
  HomogenizedTotalLagrangianStressDivergenceR(const InputParameters & parameters);

protected:
  // Add overrides to base class contributions to only happen for _beta==0, to happen only once
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobianDisplacement(unsigned int alpha, unsigned int beta) override;

  /**
   * Method for computing the scalar part of residual for _kappa
   */
  virtual void computeScalarResidual() override;

  /**
   * Method for computing the scalar variable part of Jacobian for d-_kappa-residual / d-_kappa
   */
  virtual void computeScalarJacobian() override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-jvar
   * jvar is looped over all field variables, which herein is just disp_x and disp_y
   */
  virtual void computeScalarOffDiagJacobian(const unsigned int jvar_num) override;

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

  /**
   * Method for computing an off-diagonal jacobian component d-_var-residual / d-svar.
   * svar is looped over all scalar variables, which herein is just _kappa and _kappa_other
   */
  virtual void computeOffDiagJacobianScalarLocal(const unsigned int svar_num) override;

  /**
   * Method for computing d-_var-residual / d-_svar at quadrature points.
   */
  virtual Real computeQpOffDiagJacobianScalar(const unsigned int /*svar_num*/) override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-svar
   * svar is looped over other scalar variables, which herein is just _kappa_other
   */
  virtual void computeScalarOffDiagJacobianScalar(const unsigned int svar_num) override;

protected:
  /// Which component of the scalar vector residual this constraint is responsible for
  const unsigned int _beta;

  /// (Pointer to) Scalar variable this kernel 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;

  /// Type of each constraint (stress or strain) for each component
  HomogenizationR::ConstraintMap _cmap;

  /// The constraint type; initialize with 'none'
  HomogenizationR::ConstraintType _ctype = HomogenizationR::ConstraintType::None;

  /// Used internally to iterate over each scalar component
  unsigned int _m;
  unsigned int _n;
  unsigned int _a;
  unsigned int _b;
};

(moose/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/2drow.i)

# 2D with mixed conditions on stress/strain

[GlobalParams]
  displacements = 'disp_x disp_y'
  large_kinematics = false
[]

[Mesh]
  [base]
    type = FileMeshGenerator
    file = '2d.exo'
  []

  [sidesets]
    type = SideSetsFromNormalsGenerator
    input = base
    normals = '-1 0 0
                1 0 0
                0 -1 0
                0 1 0'
    fixed_normal = true
    new_boundary = 'left right bottom top'
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [hvar]
    family = SCALAR
    order = FIRST
  []
  [hvarA]
    family = SCALAR
    order = SECOND
  []
[]

[AuxVariables]
  [sxx]
    family = MONOMIAL
    order = CONSTANT
  []
  [syy]
    family = MONOMIAL
    order = CONSTANT
  []
  [sxy]
    family = MONOMIAL
    order = CONSTANT
  []
  [exx]
    family = MONOMIAL
    order = CONSTANT
  []
  [eyy]
    family = MONOMIAL
    order = CONSTANT
  []
  [exy]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Kernels]
  [sdx]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = true
    compute_scalar_residuals = false
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_y
    component = 1
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = true
    compute_scalar_residuals = false
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sd0]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    compute_field_residuals = false
    compute_scalar_residuals = true
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sd1]
    type = HomogenizedTotalLagrangianStressDivergenceR
    variable = disp_y
    component = 1
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    compute_field_residuals = false
    compute_scalar_residuals = true
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[AuxKernels]
  [sxx]
    type = RankTwoAux
    variable = sxx
    rank_two_tensor = pk1_stress
    index_i = 0
    index_j = 0
  []
  [syy]
    type = RankTwoAux
    variable = syy
    rank_two_tensor = pk1_stress
    index_i = 1
    index_j = 1
  []
  [sxy]
    type = RankTwoAux
    variable = sxy
    rank_two_tensor = pk1_stress
    index_i = 0
    index_j = 1
  []
  [exx]
    type = RankTwoAux
    variable = exx
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 0
  []
  [eyy]
    type = RankTwoAux
    variable = eyy
    rank_two_tensor = mechanical_strain
    index_i = 1
    index_j = 1
  []
  [exy]
    type = RankTwoAux
    variable = exy
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 1
  []
[]

[Functions]
  [strain11]
    type = ParsedFunction
    value = '4.0e-2*t'
  []
  [strain22]
    type = ParsedFunction
    value = '-2.0e-2*t'
  []
  [strain12]
    type = ParsedFunction
    value = '1.0e-2*t'
  []
  [stress11]
    type = ParsedFunction
    value = '400*t'
  []
  [stress22]
    type = ParsedFunction
    value = '-200*t'
  []
  [stress12]
    type = ParsedFunction
    value = '100*t'
  []
[]

[BCs]
  [Periodic]
    [x]
      variable = disp_x
      auto_direction = 'x y'
    []
    [y]
      variable = disp_y
      auto_direction = 'x y'
    []
  []

  [fix1_x]
    type = DirichletBC
    boundary = "fix1"
    variable = disp_x
    value = 0
  []
  [fix1_y]
    type = DirichletBC
    boundary = "fix1"
    variable = disp_y
    value = 0
  []

  [fix2_y]
    type = DirichletBC
    boundary = "fix2"
    variable = disp_y
    value = 0
  []
[]

[Materials]
  [elastic_tensor_1]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 100000.0
    poissons_ratio = 0.3
    block = '1'
  []
  [elastic_tensor_2]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 120000.0
    poissons_ratio = 0.21
    block = '2'
  []
  [elastic_tensor_3]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 80000.0
    poissons_ratio = 0.4
    block = '3'
  []
  [compute_stress]
    type = ComputeLagrangianLinearElasticStress
  []
  [compute_strain]
    type = ComputeLagrangianStrain
    homogenization_gradient_names = 'homogenization_gradient'
  []
  [compute_homogenization_gradient]
    type = ComputeHomogenizedLagrangianStrainA
    macro_gradientA = hvar
    macro_gradient = hvarA
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

[Postprocessors]
  [sxx]
    type = ElementAverageValue
    variable = sxx
    execute_on = 'initial timestep_end'
  []
  [syy]
    type = ElementAverageValue
    variable = syy
    execute_on = 'initial timestep_end'
  []
  [sxy]
    type = ElementAverageValue
    variable = sxy
    execute_on = 'initial timestep_end'
  []
  [exx]
    type = ElementAverageValue
    variable = exx
    execute_on = 'initial timestep_end'
  []
  [eyy]
    type = ElementAverageValue
    variable = eyy
    execute_on = 'initial timestep_end'
  []
  [exy]
    type = ElementAverageValue
    variable = exy
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient

  solve_type = 'newton'
#  solve_type = 'PJFNK'
  line_search = none

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  l_max_its = 2
  l_tol = 1e-14
  nl_max_its = 30
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10

  start_time = 0.0
  dt = 0.2
  dtmin = 0.2
  end_time = 1.0
[]

[Outputs]
  csv = true
[]

(moose/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceA.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 "TotalLagrangianStressDivergenceS.h"

// Helpers common to the whole homogenization system
namespace HomogenizationA
{
/// Moose constraint type, for input
const MultiMooseEnum constraintType("strain stress none");
/// Constraint type: stress/PK stress or strain/deformation gradient
enum class ConstraintType
{
  Strain,
  Stress,
  None
};
typedef std::map<std::pair<unsigned int, unsigned int>, std::pair<ConstraintType, const Function *>>
    ConstraintMap;
}

/// Total Lagrangian formulation with most homogenization terms (one disp_xyz field and one scalar)
/// The macro_gradient variable is split into two scalars: the first component called '_hvar'
/// herein and all other components called '_avar' herein. For parameter _beta = 0, the primary
/// scalar (_kappa) is _hvar and the coupled scalar is _avar. For parameter _beta = 1, the primary
/// scalar (_kappa) is _avar and the coupled scalar is _hvar. Just like the primary field variable
/// (_var) is either disp_x or disp_y or disp_z depending on _alpha.
///
/// Thus, each instance of HomogenizedTotalLagrangianStressDivergenceA acts on one field variable
/// (_disp_alpha) and one scalar variable (_hvar_beta). The job of the kernel is to assemble the
/// residual of all dofs of _disp_alpha and of all dofs of _hvar_beta (namely, selected entries).
/// It assembles a symmetric portion of the Jacobian for _disp_alpha and _hvar_beta, with some
/// logical checks to access only the particular desired terms (often for _alpha or _beta = 0).
/// The entries for the other field/scalar variables are handled by other instances of the
/// kernel, which have other values of _alpha AND _beta. The logical checks ensure the proper
/// decomposition of the jobs.
///
/// In summary, for x=disp_x etc. and h=_hvar and a=_avar, then the contributions of the instances are
/// _alpha=0, _beta=0
/// R = [Rx,  00,  00,  Rh,  00 ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, 000
///      Jhx, 000, 000, Jhh, Jha]
/// _alpha=1, _beta=0
/// R = [00,  Ry,  00,  00,  00 ]^T
/// J = [Jyx, Jyy, Jyz, Jyh, 000
///      000, Jhy, 000, 000, 000]
/// _alpha=2, _beta=0
/// R = [00,  00,  Rz,  00,  00 ]^T
/// J = [Jzx, Jzy, Jzz, Jzh, 000
///      000, 000, Jhz, 000, 000]
/// _alpha=0, _beta=1
/// R = [00,  00,  00,  00,  Ra ]^T
/// J = [000, 000, 000, 000, Jxa
///      Jax, 000, 000, Jah, Jaa]
/// _alpha=1, _beta=1
/// R = [00,  00,  00,  00,  00 ]^T
/// J = [000, 000, 000, 000, Jya
///      000, Jay, 000, 000, 000]
/// _alpha=2, _beta=1
/// R = [00,  00,  00,  00,  00 ]^T
/// J = [000, 000, 000, 000, Jza
///      000, 000, Jaz, 000, 000]
///
/// In this manner, the full R and J are obtained with NO duplication of jobs:
/// R = [Rx,  Ry,  Rz,  Rh,  Ra ]^T
/// J = [Jxx, Jxy, Jxz, Jxh, Jxa
///      Jyx, Jyy, Jyz, Jyh, Jya
///      Jzx, Jzy, Jzz, Jzh, Jza
///      Jhx, Jhy, Jhz, Jhh, Jha
///      Jax, Jay, Jaz, Jah, Jaa]
///
class HomogenizedTotalLagrangianStressDivergenceA : public TotalLagrangianStressDivergenceS
{
public:
  static InputParameters validParams();
  HomogenizedTotalLagrangianStressDivergenceA(const InputParameters & parameters);

protected:
  // Add overrides to base class contributions to only happen for _beta==0, to happen only once
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobianDisplacement(unsigned int alpha, unsigned int beta) override;

  /**
   * Method for computing the scalar part of residual for _kappa
   */
  virtual void computeScalarResidual() override;

  /**
   * Method for computing the scalar variable part of Jacobian for d-_kappa-residual / d-_kappa
   */
  virtual void computeScalarJacobian() override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-jvar
   * jvar is looped over all field variables, which herein is just disp_x and disp_y
   */
  virtual void computeScalarOffDiagJacobian(const unsigned int jvar_num) override;

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

  /**
   * Method for computing an off-diagonal jacobian component d-_var-residual / d-svar.
   * svar is looped over all scalar variables, which herein is just _kappa and _kappa_other
   */
  virtual void computeOffDiagJacobianScalarLocal(const unsigned int svar_num) override;

  /**
   * Method for computing d-_var-residual / d-_svar at quadrature points.
   */
  virtual Real computeQpOffDiagJacobianScalar(const unsigned int svar_num) override;

  /**
   * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-svar
   * svar is looped over other scalar variables, which herein is just _kappa_other
   */
  virtual void computeScalarOffDiagJacobianScalar(const unsigned int svar_num) override;

protected:
  /// Which component of the scalar vector residual this constraint is responsible for
  const unsigned int _beta;

  /// (Pointer to) Scalar variable this kernel 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;

  /// Type of each constraint (stress or strain) for each component
  HomogenizationA::ConstraintMap _cmap;

  /// The constraint type; initialize with 'none'
  HomogenizationA::ConstraintType _ctype = HomogenizationA::ConstraintType::None;

  /// Used internally to iterate over each scalar component
  unsigned int _m;
  unsigned int _n;
  unsigned int _a;
  unsigned int _b;
};

(moose/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/2dsole.i)

# 2D with mixed conditions on stress/strain

[GlobalParams]
  displacements = 'disp_x disp_y'
  large_kinematics = false
[]

[Mesh]
  [base]
    type = FileMeshGenerator
    file = '2d.exo'
  []

  [sidesets]
    type = SideSetsFromNormalsGenerator
    input = base
    normals = '-1 0 0
                1 0 0
                0 -1 0
                0 1 0'
    fixed_normal = true
    new_boundary = 'left right bottom top'
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [hvar]
    family = SCALAR
    order = FIRST
  []
  [hvarA]
    family = SCALAR
    order = SECOND
  []
[]

[AuxVariables]
  [sxx]
    family = MONOMIAL
    order = CONSTANT
  []
  [syy]
    family = MONOMIAL
    order = CONSTANT
  []
  [sxy]
    family = MONOMIAL
    order = CONSTANT
  []
  [exx]
    family = MONOMIAL
    order = CONSTANT
  []
  [eyy]
    family = MONOMIAL
    order = CONSTANT
  []
  [exy]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Kernels]
  [sdx0]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_x
    component = 0
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy0]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_y
    component = 1
    macro_var = hvar
    macro_other = hvarA
    prime_scalar = 0
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdx1]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_x
    component = 0
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
  [sdy1]
    type = HomogenizedTotalLagrangianStressDivergenceA
    variable = disp_y
    component = 1
    macro_var = hvarA
    macro_other = hvar
    prime_scalar = 1
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[AuxKernels]
  [sxx]
    type = RankTwoAux
    variable = sxx
    rank_two_tensor = pk1_stress
    index_i = 0
    index_j = 0
  []
  [syy]
    type = RankTwoAux
    variable = syy
    rank_two_tensor = pk1_stress
    index_i = 1
    index_j = 1
  []
  [sxy]
    type = RankTwoAux
    variable = sxy
    rank_two_tensor = pk1_stress
    index_i = 0
    index_j = 1
  []
  [exx]
    type = RankTwoAux
    variable = exx
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 0
  []
  [eyy]
    type = RankTwoAux
    variable = eyy
    rank_two_tensor = mechanical_strain
    index_i = 1
    index_j = 1
  []
  [exy]
    type = RankTwoAux
    variable = exy
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 1
  []
[]

[Functions]
  [strain11]
    type = ParsedFunction
    value = '4.0e-2*t'
  []
  [strain22]
    type = ParsedFunction
    value = '-2.0e-2*t'
  []
  [strain12]
    type = ParsedFunction
    value = '1.0e-2*t'
  []
  [stress11]
    type = ParsedFunction
    value = '400*t'
  []
  [stress22]
    type = ParsedFunction
    value = '-200*t'
  []
  [stress12]
    type = ParsedFunction
    value = '100*t'
  []
[]

[BCs]
  [Periodic]
    [x]
      variable = disp_x
      auto_direction = 'x y'
    []
    [y]
      variable = disp_y
      auto_direction = 'x y'
    []
  []

  [fix1_x]
    type = DirichletBC
    boundary = "fix1"
    variable = disp_x
    value = 0
  []
  [fix1_y]
    type = DirichletBC
    boundary = "fix1"
    variable = disp_y
    value = 0
  []

  [fix2_y]
    type = DirichletBC
    boundary = "fix2"
    variable = disp_y
    value = 0
  []
[]

[Materials]
  [elastic_tensor_1]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 100000.0
    poissons_ratio = 0.3
    block = '1'
  []
  [elastic_tensor_2]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 120000.0
    poissons_ratio = 0.21
    block = '2'
  []
  [elastic_tensor_3]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 80000.0
    poissons_ratio = 0.4
    block = '3'
  []
  [compute_stress]
    type = ComputeLagrangianLinearElasticStress
  []
  [compute_strain]
    type = ComputeLagrangianStrain
    homogenization_gradient_names = 'homogenization_gradient'
  []
  [compute_homogenization_gradient]
    type = ComputeHomogenizedLagrangianStrainA
    macro_gradientA = hvar
    macro_gradient = hvarA
    constraint_types = ${constraint_types}
    targets = ${targets}
  []
[]

[Postprocessors]
  [sxx]
    type = ElementAverageValue
    variable = sxx
    execute_on = 'initial timestep_end'
  []
  [syy]
    type = ElementAverageValue
    variable = syy
    execute_on = 'initial timestep_end'
  []
  [sxy]
    type = ElementAverageValue
    variable = sxy
    execute_on = 'initial timestep_end'
  []
  [exx]
    type = ElementAverageValue
    variable = exx
    execute_on = 'initial timestep_end'
  []
  [eyy]
    type = ElementAverageValue
    variable = eyy
    execute_on = 'initial timestep_end'
  []
  [exy]
    type = ElementAverageValue
    variable = exy
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient

  solve_type = 'newton'
#  solve_type = 'PJFNK'
  line_search = none

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  l_max_its = 2
  l_tol = 1e-14
  nl_max_its = 30
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10

  start_time = 0.0
  dt = 0.2
  dtmin = 0.2
  end_time = 1.0
[]

[Outputs]
  csv = true
[]

(moose/test/tests/kernels/scalar_kernel_constraint/scalar_constraint_kernel.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  nx = 2
  ny = 2
  elem_type = QUAD9
[]

[Functions]
  [exact_fn]
    type = ParsedFunction
    value = 'x*x+y*y'
  []

  [ffn]
    type = ParsedFunction
    value = -4
  []

  [bottom_bc_fn]
    type = ParsedFunction
    value = -2*y
  []

  [right_bc_fn]
    type = ParsedFunction
    value =  2*x
  []

  [top_bc_fn]
    type = ParsedFunction
    value =  2*y
  []

  [left_bc_fn]
    type = ParsedFunction
    value = -2*x
  []
[]

[Variables]
  [u]
    family = LAGRANGE
    order = SECOND
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
[]

[Kernels]
  [diff]
    type = Diffusion
    variable = u
  []

  [ffnk]
    type = BodyForce
    variable = u
    function = ffn
  []

  [sk_lm]
    type = ScalarLMKernel
    variable = u
    kappa = lambda
    pp_name = pp
    value = 2.666666666666666
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [bottom]
    type = FunctionNeumannBC
    variable = u
    boundary = 'bottom'
    function = bottom_bc_fn
  []
  [right]
    type = FunctionNeumannBC
    variable = u
    boundary = 'right'
    function = right_bc_fn
  []
  [top]
    type = FunctionNeumannBC
    variable = u
    boundary = 'top'
    function = top_bc_fn
  []
  [left]
    type = FunctionNeumannBC
    variable = u
    boundary = 'left'
    function = left_bc_fn
  []
[]

[Postprocessors]
  # integrate the volume of domain since original objects set
  # int(phi)=V0, rather than int(phi-V0)=0
  [pp]
    type = FunctionElementIntegral
    function = 1
    execute_on = initial
  []
  [l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
    execute_on = 'initial timestep_end'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
    solve_type = 'NEWTON'
  []
[]

[Executioner]
  type = Steady
  nl_rel_tol = 1e-9
  l_tol = 1.e-10
  nl_max_its = 10
  # This example builds an indefinite matrix, so "-pc_type hypre -pc_hypre_type boomeramg" cannot
  # be used reliably on this problem
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  # This is a linear problem, so we don't need to recompute the
  # Jacobian. This isn't a big deal for a Steady problems, however, as
  # there is only one solve.
  solve_type = 'LINEAR'
[]

[Outputs]
  exodus = true
  hide = lambda
[]

(moose/test/tests/kernels/scalar_kernel_constraint/scalar_constraint_kernel.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  nx = 2
  ny = 2
  elem_type = QUAD9
[]

[Functions]
  [exact_fn]
    type = ParsedFunction
    value = 'x*x+y*y'
  []

  [ffn]
    type = ParsedFunction
    value = -4
  []

  [bottom_bc_fn]
    type = ParsedFunction
    value = -2*y
  []

  [right_bc_fn]
    type = ParsedFunction
    value =  2*x
  []

  [top_bc_fn]
    type = ParsedFunction
    value =  2*y
  []

  [left_bc_fn]
    type = ParsedFunction
    value = -2*x
  []
[]

[Variables]
  [u]
    family = LAGRANGE
    order = SECOND
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
[]

[Kernels]
  [diff]
    type = Diffusion
    variable = u
  []

  [ffnk]
    type = BodyForce
    variable = u
    function = ffn
  []

  [sk_lm]
    type = ScalarLMKernel
    variable = u
    kappa = lambda
    pp_name = pp
    value = 2.666666666666666
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [bottom]
    type = FunctionNeumannBC
    variable = u
    boundary = 'bottom'
    function = bottom_bc_fn
  []
  [right]
    type = FunctionNeumannBC
    variable = u
    boundary = 'right'
    function = right_bc_fn
  []
  [top]
    type = FunctionNeumannBC
    variable = u
    boundary = 'top'
    function = top_bc_fn
  []
  [left]
    type = FunctionNeumannBC
    variable = u
    boundary = 'left'
    function = left_bc_fn
  []
[]

[Postprocessors]
  # integrate the volume of domain since original objects set
  # int(phi)=V0, rather than int(phi-V0)=0
  [pp]
    type = FunctionElementIntegral
    function = 1
    execute_on = initial
  []
  [l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
    execute_on = 'initial timestep_end'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
    solve_type = 'NEWTON'
  []
[]

[Executioner]
  type = Steady
  nl_rel_tol = 1e-9
  l_tol = 1.e-10
  nl_max_its = 10
  # This example builds an indefinite matrix, so "-pc_type hypre -pc_hypre_type boomeramg" cannot
  # be used reliably on this problem
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  # This is a linear problem, so we don't need to recompute the
  # Jacobian. This isn't a big deal for a Steady problems, however, as
  # there is only one solve.
  solve_type = 'LINEAR'
[]

[Outputs]
  exodus = true
  hide = lambda
[]