Conditions

Sources. Thesis §4.3.3.4 (discretization and mortar operators, pp. 100–113), §4.3.4.3 (frictional algebraic form), §4.3.5 (mesh tying); code: custom_conditions/paired_condition.{h,cpp}, custom_conditions/mortar_contact_condition.{h,cpp}, custom_conditions/ALM_frictionless_mortar_contact_condition.h, custom_conditions/ALM_frictionless_components_mortar_contact_condition.h, custom_conditions/ALM_frictional_mortar_contact_condition.h, custom_conditions/penalty_frictionless_mortar_contact_condition.h, custom_conditions/penalty_frictional_mortar_contact_condition.h, the four *_axisym_condition.{h,cpp}, custom_conditions/mesh_tying_mortar_condition.{h,cpp}, custom_conditions/mpc_mortar_contact_condition.{h,cpp}, contact_structural_mechanics_application.{h,cpp} (prototypes and Register()), automatic_differentiation/*/*_template.cpp, kratos/includes/mortar_classes.h, custom_utilities/mortar_explicit_contribution_utilities.{h,cpp}.

A condition in Kratos is the object that contributes a local left-hand side (LHS) and right-hand side (RHS) to the global system for a piece of boundary. In this application every contact pair (one slave face and one master face that the search has matched) is a condition: it owns both geometries, integrates the mortar operators \(\mathbf{D}\) and \(\mathbf{M}\) over the exact intersection of the two faces, and evaluates the weak contact constraints and their consistent linearization for its slave nodes. This page describes how these classes are organized, how a local system is computed, what differs between the formulation families, how the registered names are built, and which values the conditions read from the ProcessInfo and the Properties. The mathematics is documented in the Theory section: Frictionless contact, Frictional contact, Mesh tying, Mortar integration and dual Lagrange multipliers and Automatic differentiation.

Class hierarchy

Figure: the condition classes of `custom_conditions/`, their template parameters and the registered-name decoder.

Condition                                          (Kratos core, includes/condition.h)
└── PairedCondition                                paired_condition.{h,cpp}
    ├── MortarContactCondition<TDim, TNumNodes, TFrictional, TNormalVariation, TNumNodesMaster>
    │   │                                          mortar_contact_condition.{h,cpp}
    │   ├── AugmentedLagrangianMethodFrictionlessMortarContactCondition            FRICTIONLESS
    │   │   └── AugmentedLagrangianMethodFrictionlessMortarContactAxisymCondition
    │   ├── AugmentedLagrangianMethodFrictionlessComponentsMortarContactCondition  FRICTIONLESS_COMPONENTS
    │   ├── AugmentedLagrangianMethodFrictionalMortarContactCondition              FRICTIONAL
    │   │   └── AugmentedLagrangianMethodFrictionalMortarContactAxisymCondition
    │   ├── PenaltyMethodFrictionlessMortarContactCondition                        FRICTIONLESS_PENALTY
    │   │   └── PenaltyMethodFrictionlessMortarContactAxisymCondition
    │   └── PenaltyMethodFrictionalMortarContactCondition                          FRICTIONAL_PENALTY
    │       └── PenaltyMethodFrictionalMortarContactAxisymCondition
    ├── MeshTyingMortarCondition<TDim, TNumNodes, TNumNodesMaster>
    └── MPCMortarContactCondition<TDim, TNumNodes, TNumNodesMaster>

Three layers can be distinguished:

1. PairedCondition (non-templated) only adds to Condition the notion of a second geometry and a cached normal for it.
2. MortarContactCondition contains everything that is common to the contact formulations: the segmentation of the pair, the integration of the mortar operators and of their derivatives, the isolation check, the integration-method selection, and the dispatch to the family-specific CalculateLocalLHS / CalculateLocalRHS. It never assembles anything itself: its own CalculateLocalLHS / CalculateLocalRHS, EquationIdVector, GetDofList and GetActiveInactiveValue raise KRATOS_ERROR (mortar_contact_condition.cpp:537-590).
3. The families are specializations of the FrictionalCase template parameter. Their headers are hand-written (DoF lists, MatrixSize, state encoding, extra members), while the bodies of CalculateLocalLHS / CalculateLocalRHS in the .cpp files are generated by the sympy scripts of automatic_differentiation/ and are, therefore, very large (from 27 480 lines for the penalty frictionless condition to 170 787 lines for the ALM frictional one). The axisymmetric variants derive from the 2D instantiation of a family and override only the integration weight.

MeshTyingMortarCondition and MPCMortarContactCondition do not derive from MortarContactCondition: the first one because its operators are constant and computed once (there is no active set, no derivative of the operators), the second one because it does not contribute to the system at all but feeds a multipoint constraint.

Template parameters

The pairs <TDim, TNumNodes, TNumNodesMaster> that exist are <2,2,2>, <3,3,3>, <3,4,4>, <3,3,4> and <3,4,3>; combined with the five FrictionalCase values and the two values of TNormalVariation, mortar_contact_condition.cpp:751-808 instantiates the 50 base classes explicitly. The BelongType used by the segmentation to remember which slave/master edge produced each cell vertex is also chosen from these parameters (PointBelongsLine2D2N, PointBelongsTriangle3D3N, PointBelongsQuadrilateral3D4N, PointBelongsTriangle3D3NQuadrilateral3D4N, PointBelongsQuadrilateral3D4NTriangle3D3N, all in kratos/includes/mortar_classes.h).

The FrictionalCase enum

Defined in contact_structural_mechanics_application_variables.h:

enum class FrictionalCase {
    FRICTIONLESS            = 0,   // ALM, scalar normal Lagrange multiplier
    FRICTIONLESS_COMPONENTS = 1,   // ALM, vector Lagrange multiplier, tangential part penalized to zero
    FRICTIONAL              = 2,   // ALM, vector Lagrange multiplier, Coulomb friction
    FRICTIONLESS_PENALTY    = 3,   // pure penalty, displacement DoFs only
    FRICTIONAL_PENALTY      = 4    // pure penalty with friction, displacement DoFs only
};

The enum drives compile-time decisions in mortar_contact_condition.h:140-159:

$$ \texttt{MatrixSize} = \begin{cases} d\,(n_s + n_m) + n_s & \texttt{FRICTIONLESS} \\ d\,(2 n_s + n_m) & \texttt{FRICTIONLESS\_COMPONENTS},\ \texttt{FRICTIONAL} \\ d\,(n_s + n_m) & \texttt{FRICTIONLESS\_PENALTY},\ \texttt{FRICTIONAL\_PENALTY} \end{cases} $$

with \(d\) = TDim, \(n_s\) = TNumNodes and \(n_m\) = TNumNodesMaster, and IsFrictional = (TFrictional == FRICTIONAL || TFrictional == FRICTIONAL_PENALTY), which selects the MortarOperatorWithDerivatives<TDim, TNumNodes, IsFrictional, TNumNodesMaster> and DualLagrangeMultiplierOperatorsWithDerivatives<...> types (the frictional versions also store the derivatives needed by the slip). A second enum in the same header, NormalDerivativesComputation, is read at run time from CONSIDER_NORMAL_VARIATION (see PairedCondition below and the Variables and flags reference).

PairedCondition

PairedCondition (paired_condition.h) is “basically equal to the base condition, with a pointer to the paired geometry”. The two faces are stored in a single Kratos-core CouplingGeometry<Node> (using CouplingGeometryType = CouplingGeometry<Node>, paired_condition.h:71), created in the constructors as Kratos::make_shared<CouplingGeometryType>(pGeometry, pPairedGeometry) (paired_condition.h:104-126). The accessors are:

Note (naming inversion). CouplingGeometry calls its first geometry Master and its second one Slave, in the sense of “the geometry that owns the integration”. Contact mechanics (and this application) use the opposite convention, so CouplingGeometryType::Master is the slave face and CouplingGeometryType::Slave is the master face. Always use GetParentGeometry() (slave) and GetPairedGeometry() (master) and never the CouplingGeometry enumerators directly. The slave normal, by contrast, is not cached in the condition: it is read from the condition’s own NORMAL value (this->GetValue(NORMAL), mortar_contact_condition.cpp:302), which the search process sets.

The three life-cycle hooks of PairedCondition (paired_condition.cpp:68-120) only manage mPairedNormal:

• Initialize and InitializeSolutionStep recompute it from GetPairedGeometry().UnitNormal(center);
• InitializeNonLinearIteration recomputes it only when rCurrentProcessInfo[CONSIDER_NORMAL_VARIATION] (cast to NormalDerivativesComputation) is different from NO_DERIVATIVES_COMPUTATION, so that with ELEMENTAL_DERIVATIVES, NODAL_ELEMENTAL_DERIVATIVES or NO_DERIVATIVES_COMPUTATION_WITH_NORMAL_UPDATE the master normal follows the deformation inside the Newton loop.

mPairedNormal is serialized (save/load, paired_condition.cpp:125-138). PairedCondition has three Create overloads; the one taking (NewId, pGeometry, pProperties, pPairedGeom) is what the search uses to build a pair from the two skin conditions (mpReferenceCondition->Create(rConditionId, pObjectSlave->pGetGeometry(), pProperties, pObjectMaster->pGetGeometry()), base_contact_search_process.cpp:811).

Life-cycle hooks

Kratos calls five hooks on every condition during a time step (Architecture). The table summarizes what each class does in them; a dash means “only forwards to the base class”.

The assembly methods (CalculateLocalSystem, CalculateLeftHandSide, CalculateRightHandSide) are called by the builder-and-solver inside each Newton iteration, and AddExplicitContribution(const ProcessInfo&) is called by ContactUtilities::ComputeExplicitContributionConditions before the active-set check; both are described next. Serialization (save/load) stores mPairedNormal in the base, and mPreviousMortarOperators plus mPreviousMortarOperatorsInitialized in the frictional and MPC classes (ALM_frictional_mortar_contact_condition.h:525-533, mpc_mortar_contact_condition.h:592-599), so that a restart resumes with a consistent slip.

MortarContactCondition: how a local system is computed

CalculateLocalSystem, CalculateLeftHandSide and CalculateRightHandSide (mortar_contact_condition.cpp:157-208) just resize the local containers to MatrixSize (ResizeLHS / ResizeRHS) and call the master driver CalculateConditionSystem(rLHS, rRHS, rProcessInfo, ComputeLHS, ComputeRHS) (mortar_contact_condition.cpp:290-434). CalculateMassMatrix and CalculateDampingMatrix resize to 0x0: a contact condition has no inertia. The driver executes the following steps.

CalculateConditionSystem(LHS, RHS, ProcessInfo, ComputeLHS, ComputeRHS)        mortar_contact_condition.cpp:290
  1. slave = GetParentGeometry(); n_slave = GetValue(NORMAL)
     derivative_data.Initialize(slave, ProcessInfo)        <- u1, X1, NormalSlave, PenaltyParameter[i] = node.GetValue(INITIAL_PENALTY),
                                                              ScaleFactor = ProcessInfo[SCALE_FACTOR] (+ TangentFactor, u1old if frictional)
     if TNormalVariation: DerivativesUtilities::CalculateDeltaNormalSlave(...)
  2. integration_utility = ExactMortarIntegrationUtility(INTEGRATION_ORDER_CONTACT, DISTANCE_THRESHOLD,
                                                         0, ZERO_TOLERANCE_FACTOR, CONSIDER_TESSELLATION)
     master = GetPairedGeometry(); n_master = GetPairedNormal()
  3. is_inside = CheckIsolatedElement(DELTA_TIME) ? false
                 : integration_utility.GetExactIntegration(slave, n_slave, master, n_master, cells)   l. 334
     integration_area = GetTotalArea(slave, cells)
  4. if is_inside and integration_area / slave.Area() > 1e-5:                                          l. 340
       derivative_data.UpdateMasterPair(master, ProcessInfo)      <- u2, X2 (+ u2old)
       if TNormalVariation: CalculateDeltaNormalMaster(...)
       dual_LM = DerivativesUtilities::CalculateAeAndDeltaAe(slave, n_slave, master, derivative_data,
                    general_variables, consider_normal_variation, cells, method, axisym_coefficient)  l. 354
       for each cell in cells:                                                                         l. 358
         decomp_geom = Line2D2 / Triangle3D3 built from the cell vertices (global coordinates)
         skip if bad shape (LengthCheck in 2D with CheckThresholdCoefficient = 1e-12, HeronCheck in 3D)
         for each Gauss point of decomp_geom (GetIntegrationMethod()):                                 l. 380
           derivative_data.ResetDerivatives()
           map GP -> slave local coordinates (local_point_parent)
           MortarExplicitContributionUtilities::CalculateKinematics(...)        <- N_slave, Phi (dual), N_master, DetJ  l. 390
           w = GP.Weight() * GetAxisymmetricCoefficient(general_variables)
           if ComputeLHS:
             (3D) CalculateDeltaCellVertex(...); CalculateDeltaDetjSlave(...); CalculateDeltaN(...)
             mortar_operators.CalculateDeltaMortarOperators(general_variables, derivative_data, w)      l. 403
           else:
             mortar_operators.CalculateMortarOperators(general_variables, w)                            l. 405
       active_inactive = GetActiveInactiveValue(slave)                                                 l. 411
       if ComputeLHS: CalculateLocalLHS(LHS, mortar_operators, derivative_data, active_inactive, ProcessInfo)
       if ComputeRHS: CalculateLocalRHS(RHS, mortar_operators, derivative_data, active_inactive, ProcessInfo)
     else:                                                                                             l. 421
       Set(ISOLATED, true); ZeroLHS(LHS); ZeroRHS(RHS)

Some remarks on each stage:

Segmentation (step 3). IntegrationUtility is ExactMortarIntegrationUtility<TDim, TNumNodes, true, TNumNodesMaster> (mortar_contact_condition.h:156); the true asks the utility to return, for every vertex of every cell, the belonging information (which slave/master node or edge intersection generated it), needed later by CalculateDeltaCellVertex to differentiate the cell geometry with respect to the displacements. GetExactIntegration projects the master onto the slave plane along the slave normal, clips the two polygons and triangulates the intersection (see Mortar integration and dual Lagrange multipliers). The result is a list of TDim -vertex cells expressed in the slave local coordinates. The distance_threshold (default 1.0e24) and the zero_tolerance_factor (default 1.0) come from the ProcessInfo; the INTEGRATION_ORDER_CONTACT (default 2) and CONSIDER_TESSELLATION (default false) from the pair Properties.

Isolation. CheckIsolatedElement (mortar_contact_condition.cpp:440-532) currently only returns this->Is(ISOLATED) (a long heuristic based on the relative motion of the two faces is kept commented out). Hence, once a pair is flagged ISOLATED (because its intersection is empty or its area is below \(10^{-5}\) of the slave area) it stays zeroed until Initialize resets the flag (mortar_contact_condition.cpp:95), which happens when the search rebuilds the pairs. Isolated slave nodes are handled by the block builder-and-solver (Builder and solvers).

Dual shape functions (step 4, Ae). CalculateAeAndDeltaAe (custom_utilities/derivatives_utilities.h:345) integrates \(\mathbf{D}_e\) and \(\mathbf{M}_e\) over the cells (Popp’s eq. 3.65) and returns the coefficient matrix \(\mathbf{A}_e = \mathbf{D}_e \mathbf{M}_e^{-1}\) that turns the standard shape functions \(N_k\) into the dual ones \(\Phi_j = \sum_k A_{e,jk} N_k\), together with its directional derivatives \(\Delta \mathbf{A}_e\) (Popp’s eq. 4.58); the boolean dual_LM it returns is false when the matrix could not be inverted safely, in which case the standard shape functions are used for the multipliers. The axisymmetric coefficient is passed so that \(\mathbf{A}_e\) is consistent with the weighted integration.

Cell and Gauss loops. Each cell is a Line2D2<Point> or Triangle3D3<Point> in global coordinates; degenerate cells are skipped (MortarUtilities::LengthCheck / HeronCheck, mortar_contact_condition.cpp:369-373). For every Gauss point of the cell (GetIntegrationMethod() below), CalculateKinematics evaluates the slave shape functions \(N^{(1)}\) and the dual ones \(\Phi\) at the parent local point, projects the point onto the master to get \(N^{(2)}\), and computes the Jacobian determinant of the cell (MortarKinematicVariablesWithDerivatives, mortar_classes.h:444). When the LHS is requested the derivatives of the cell vertices (3D only), of the Jacobian and of all the shape functions are updated first and CalculateDeltaMortarOperators accumulates both the operators and their derivatives (MortarOperatorWithDerivatives, mortar_classes.h:1195-1297, Popp’s eqs. 4.32–4.38):

$$ D_{jk} \mathrel{+}= w \, \Phi_j N^{(1)}_k \, \det J, \qquad M_{jl} \mathrel{+}= w \, \Phi_j N^{(2)}_l \, \det J, \qquad \Delta D_{jk}, \Delta M_{jl} \text{ for each of the } d\,(n_s + n_m) \text{ displacement DoFs}. $$

For a pure RHS evaluation only CalculateMortarOperators is called: the residual does not need the derivatives, which makes CalculateRightHandSide (used, for example, by the line-search strategy and by the convergence criteria) much cheaper than CalculateLocalSystem.

Active/inactive encoding and dispatch (steps GetActiveInactiveValue and CalculateLocalLHS/RHS). The base class declares GetActiveInactiveValue as virtual and the families implement it (see the table in the next section). The integer is passed as rActiveInactive to the family-specific CalculateLocalLHS(rLocalLHS, rMortarConditionMatrices, rDerivativeData, rActiveInactive, rCurrentProcessInfo) and CalculateLocalRHS(...) (mortar_contact_condition.h:531-552), whose bodies are the sympy-generated code. The interface has kept this integer since the first generators, which emitted if (rActiveInactive == N) {...} else if ... chains (generate_frictionless_mortar_condition.py:187-200); the code that is currently committed comes from the _non_zero generators and from the frictional generator, which branch per node directly on the flags (if (r_geometry[i].IsNot(ACTIVE)) { ... } else if (r_geometry[i].Is(SLIP)) {...}), so the integer is computed but the generated bodies do not test it. See Automatic differentiation for the generation pipeline and the BEGIN/END AD REPLACEMENT markers.

GetIntegrationMethod and INTEGRATION_ORDER_CONTACT

GetIntegrationMethod() (mortar_contact_condition.h:583-601, identical in mesh_tying_mortar_condition.h:638-650) maps the INTEGRATION_ORDER_CONTACT property of the pair to the Gauss rule used on every integration cell:

The same integer is handed to ExactMortarIntegrationUtility, which uses the same switch (kratos/utilities/exact_mortar_segmentation_utility.cpp:1324-1346) for the auxiliary integration it performs during the clipping. Because the integrand on a cell is a product of polynomials of low degree (the shape functions of both sides are linear/bilinear and the cell Jacobian is constant on a triangle), the default order 2 is exact for the frictionless operators on linear geometries; higher orders only matter for quadrilaterals and for the axisymmetric weight. The property is written by SearchBaseProcess from the integration_order key of the contact process settings (Contact process settings reference).

The formulation families

Figure: nodal DoFs, historical variables, non-historical values and flags used by each `mortar_type`; the conditions below only see the slave/master displacements and the slave multipliers listed in their DoF list.

Every family orders its local system as [master displacements, slave displacements, slave multipliers] (/* ORDER - [ MASTER, SLAVE, LAMBDA ] */ in every EquationIdVector, for example ALM_frictional_mortar_contact_condition_template.cpp:174-215), with the components of each node consecutive (X, Y, and Z in 3D).

Additional facts that the table cannot hold:

• Where the ALM parameters come from. The generated bodies read ScaleFactor = rDerivativeData.ScaleFactor, PenaltyParameter = rDerivativeData.PenaltyParameter (one value per slave node), TangentFactor = rDerivativeData.TangentFactor (frictional only), DynamicFactor = GetVariableVector(slave, DYNAMIC_FACTOR) and LM/LMNormal from the nodal multiplier variables; DerivativeData::Initialize fills them from node.GetValue(INITIAL_PENALTY), rCurrentProcessInfo[SCALE_FACTOR] and rCurrentProcessInfo[TANGENT_FACTOR] (mortar_classes.h:709-710, 960). The generated frictionless code therefore evaluates, for an active node \(j\), the augmented pressure \(\bar{\lambda}_{n,j} = k \lambda_{n,j} + \varepsilon_j \tilde{g}_{n,j}\) with \(\tilde{g}_{n,j} = \mathbf{n}_j \cdot (\mathbf{D} \mathbf{x}^{(1)} - \mathbf{M} \mathbf{x}^{(2)})_j\) and, for an inactive node, the term \(-\frac{k^2}{\varepsilon_j} \lambda_{n,j}\) that returns the multiplier to zero (Frictionless contact, thesis §4.3.3.4.3).
• Consistent slip (frictional families). The slip increment needs the operators of the previous converged step: InitializeSolutionStep calls ComputePreviousMortarOperators once per step (guarded by mPreviousMortarOperatorsInitialized, which FinalizeSolutionStep resets; ALM_frictional_mortar_contact_condition_template.cpp:89-117). The generated code then compares the Frobenius norms of \(\mathbf{D}-\mathbf{D}_{old}\) and \(\mathbf{M}-\mathbf{M}_{old}\) with rCurrentProcessInfo[OPERATOR_THRESHOLD]: if both exceed it the objective slip \((\mathbf{D}-\mathbf{D}_{old})\mathbf{x}^{(1)} - (\mathbf{M}-\mathbf{M}_{old})\mathbf{x}^{(2)}\) is used, otherwise the non-objective \(\mathbf{D}(\mathbf{x}^{(1)}-\mathbf{x}^{(1)}_{old}) - \mathbf{M}(\mathbf{x}^{(2)}-\mathbf{x}^{(2)}_{old})\), and the condition flag MODIFIED records the choice (this->Set(MODIFIED, !is_objetive), ALM_frictional_mortar_contact_condition.cpp:201-205). The friction coefficient is read per slave node with GetFrictionCoefficient() (FRICTION_COEFFICIENT non-historical value, ALM_frictional_mortar_contact_condition.h:456-467), still marked // TODO: Define the "CL" or friction law; the FrictionalLaw classes are described in Frictional laws and MPC constraint. See Frictional contact for the stick/slip branches (thesis §4.3.4.3).
• Penalty families. Without multipliers, the augmented pressure reduces to \(\varepsilon_j \tilde{g}_{n,j}\) and the stick force to TangentFactor\(\cdot \varepsilon_j\) times the slip; the penalty frictionless AddExplicitContribution additionally accumulates NODAL_AREA, which the penalty process uses to scale the nodal penalty. These conditions are the ones used by the explicit dynamic solver.
• Explicit contributions. AddExplicitContribution(const ProcessInfo&) is the only method allowed to write on the nodes; ContactUtilities::ComputeExplicitContributionConditions calls it on the ComputingContact model part at the beginning of every step and iteration to refresh the weighted gap (and slip) that the active-set check and the Predict of the strategy need (Architecture). The overloads taking a destination variable (mortar_contact_condition.cpp:264-285) compute the nodal area into a double variable via MortarExplicitContributionUtilities::ComputeNodalArea and raise an error for array_1d destinations.

Note. PenaltyMethodFrictionalMortarContactCondition::AddExplicitContribution passes false, false after mPreviousMortarOperators (penalty_frictional_mortar_contact_condition.cpp:129), so the IntegrationOrder argument receives 0 and ExactMortarIntegrationUtility falls back to GI_GAUSS_2: the explicit gap/slip of this family ignores INTEGRATION_ORDER_CONTACT (the implicit assembly does honor it). Harmless with the default order.

Local DoF ordering

EquationIdVector and GetDofList fill the local vector in the order master displacements, slave displacements, slave multipliers, node by node and component by component. For the smallest ALM frictionless condition, ALMFrictionlessMortarContactCondition2D2N (MatrixSize = 2·(2+2)+2 = 10), the local indices are

index : 0   1   2   3   4   5   6   7   8        9
DoF   : u2x u2y u2x u2y u1x u1y u1x u1y lambda_n lambda_n
node  : m0  m0  m1  m1  s0  s0  s1  s1  s0       s1

where m denotes master nodes, s slave nodes, u2 the master displacement, u1 the slave displacement and lambda_n the LAGRANGE_MULTIPLIER_CONTACT_PRESSURE DoF (ALM_frictionless_mortar_contact_condition_template.cpp:96-133). For the vector-multiplier families the last block has TDim entries per slave node (VECTOR_LAGRANGE_MULTIPLIER_X, _Y, _Z); for the penalty families it is absent. The generated code addresses rLocalLHS(i, j) with these absolute indices, which is why every specialization is a separate function with hard-coded sizes and why a change in this ordering requires regenerating all the .cpp files. The mesh-tying condition follows the same order with dof_size components per node, and the MPC condition only has the two displacement blocks.

The resulting sizes for every registered geometry pair are:

Active/inactive and stick/slip encodings

The two GetActiveInactiveValue encodings turn the nodal flags of the slave face into one integer:

• Frictionless families (ALM_frictionless_mortar_contact_condition.h:385-393, identical in the components and penalty frictionless headers): value += 1 << i_node for every node that Is(ACTIVE). For a triangle with nodes 0 and 2 active the value is \(2^0 + 2^2 = 5\); the range is \(0 \ldots 2^{n_s}-1\) (0–15 for a quadrilateral).
• Frictional families (ALM_frictional_mortar_contact_condition.h:437-450, penalty_frictional_mortar_contact_condition.h:471-484): value += pow(3, i_node) for an active node that Is(SLIP) and value += 2 * pow(3, i_node) for an active node that is not SLIP (stick); inactive nodes contribute 0. For a triangle with node 0 sticking, node 1 slipping and node 2 inactive the value is \(2 \cdot 3^0 + 1 \cdot 3^1 = 5\); the range is \(0 \ldots 3^{n_s}-1\) (0–80 for a quadrilateral).

The committed generated code branches per node instead of on the integer. The frictional structure, repeated for every slave node i in both CalculateLocalLHS and CalculateLocalRHS, is

// NODE i
if (r_geometry[i].IsNot(ACTIVE)) {            // INACTIVE: -k^2/eps * lambda_i  (returns the multiplier to zero)
    ...
} else if (r_geometry[i].Is(SLIP)) {          // ACTIVE-SLIP
    if (is_objetive) { ... }                  //   objective slip     (D - D_old) x1 - (M - M_old) x2
    else             { ... }                  //   non-objective slip D (x1 - x1_old) - M (x2 - x2_old)
} else {                                      // ACTIVE-STICK
    if (is_objetive) { ... } else { ... }
}

(ALM_frictional_mortar_contact_condition.cpp:207-1327 for node 0 of the <2,2,false,2> specialization; the // ACTIVE-SLIP, // OBJECTIVE-SLIP and // ACTIVE-STICK comments mark the branches). Five branches per node times up to four nodes, each containing the fully expanded linearization, is what makes this file 170 787 lines long. The frictionless families have two branches per node (inactive / active). The is_objetive boolean and the MODIFIED flag are computed once per call, before the node loop (l. 201-205).

Axisymmetric variants

The four *_axisym_condition.{h,cpp} classes are templated on <TNumNodes, TNormalVariation> (plus TNumNodesMaster for the penalty frictional one) and derive from the <2, TNumNodes, TNormalVariation> instantiation of their family; only <2, false> and <2, true> are instantiated. They override three methods (ALM_frictionless_mortar_contact_axisym_condition.cpp:74-113, same lines in the other three files):

• IsAxisymmetric() returns true;
• GetAxisymmetricCoefficient(rVariables) returns \(2 \pi r / t\) with t = GetProperties()[THICKNESS];
• CalculateRadius(rVariables) returns \(r = \sum_i N^{(1)}_i x_i\), the current \(x\) coordinate of the Gauss point on the slave (the reference-configuration alternative is kept commented out).

Because the base class multiplies every Gauss weight by GetAxisymmetricCoefficient (mortar_contact_condition.cpp:393) and passes it to CalculateAeAndDeltaAe, the operators \(\mathbf{D}\), \(\mathbf{M}\) and \(\mathbf{A}_e\) become integrals over the revolved surface, consistent with the axisymmetric solid elements that divide by THICKNESS. The base implementation returns 1.0 (mortar_contact_condition.cpp:703-706). The Axisym infix in the registered names selects these classes; the Python processes pick them when axisymmetric is true in the contact process settings (alm_contact_process.py:325-346, penalty_contact_process.py:197-213).

Note (MatrixSize shadowing). Each axisymmetric header re-declares static constexpr IndexType MatrixSize instead of inheriting it. Two of them disagree with their base: ALM_frictional_mortar_contact_axisym_condition.h:119 declares 2 * (TNumNodes + TNumNodes) + TNumNodes (= 10 for two nodes) while the frictional base is TDim * (TNumNodes + TNumNodes + TNumNodesMaster) (= 12), and penalty_frictionless_mortar_contact_axisym_condition.h:120 declares the same 2 * (TNumNodes + TNumNodes) + TNumNodes (= 10) while the penalty base is TDim * (TNumNodes + TNumNodesMaster) (= 8). The mismatch is harmless because the derived constants are never used: ResizeLHS, ResizeRHS, ZeroLHS and ZeroRHS are implemented in the base class with the base constant (mortar_contact_condition.cpp:711-740), and the generated CalculateLocalLHS/RHS are specializations of the base family with hard-coded sizes.

Registered names

KratosContactStructuralMechanicsApplication::Register() (contact_structural_mechanics_application.cpp:123-261) registers 68 conditions and one constraint (ContactMasterSlaveConstraint). The prototypes are members of the application class (contact_structural_mechanics_application.h) built in its constructor with LineType = Line2D2<Node>, TriangleType = Triangle3D3<Node> and QuadrilateralType = Quadrilateral3D4<Node> geometries for the slave and the master (contact_structural_mechanics_application.cpp:25-120).

Decoder

A contact condition name is assembled as

$$ \underbrace{\texttt{ALM} \;\vert\; \texttt{Penalty}}_{\text{method}} \; \underbrace{[\texttt{NV}]}_{\text{normal variation}} \; \underbrace{\texttt{Frictionless} \;\vert\; \texttt{FrictionlessComponents} \;\vert\; \texttt{Frictional}}_{\text{family}} \; \underbrace{[\texttt{Axisym}]}_{\text{2D only}} \; \texttt{MortarContactCondition} \; \underbrace{\texttt{2D2N} \;\vert\; \texttt{3D3N} \;\vert\; \texttt{3D4N} \;\vert\; \texttt{3D3N4N} \;\vert\; \texttt{3D4N3N}}_{\text{geometry pair}} $$

and the two non-contact families are MeshTyingMortarCondition<pair> and MPCMortarContactCondition<pair>. The geometry suffix reads <TDim>D<TNumNodes>N[<TNumNodesMaster>N], the master count being appended only when it differs from the slave count. BaseContactSearchProcess builds exactly this string at run time: condition_name << "Condition" << TDim << "D" << TNumNodes << "N" << final_string (base_contact_search_process.cpp:99), where the Python process supplies the prefix (_get_condition_name, for example "ALMNVFrictionalAxisymMortarContact") and final_string = "<n_master>N" when the two sides have a different number of nodes (search_base_process.py:307-320).

Full list

There is no axisymmetric FrictionlessComponents condition, no axisymmetric 3D condition (axisymmetry is a 2D concept) and no NV variant of the mesh-tying or MPC conditions (they do not linearize the normals). Users never write these names in ProjectParameters.json: the contact processes derive them from contact_type, normal_variation and axisymmetric (Contact process settings reference).

Note (pairs cannot come from an .mdpa). A Begin Conditions <name> block of an .mdpa file creates conditions through Create(NewId, rThisNodes, pProperties), which for a PairedCondition wraps the nodes in a CouplingGeometry whose second geometry is nullptr (paired_condition.h:104-114). Such a condition fails Check() with YOU HAVE NOT INITIALIZED THE PAIR GEOMETRY IN THE MortarContactCondition. Pairs must therefore always be created by the search (or, in C++ tests, with r_prototype.Create(id, p_slave_geometry, p_properties, p_master_geometry) as in tests/cpp_tests/conditions/test_mesh_tying_condition.cpp:280). The .mdpa only provides the skin conditions (any surface/line condition of the StructuralMechanicsApplication) from which the search builds the pairs.

Note (mesh-tying prototypes). The prototypes mMeshTyingMortarCondition3D3N4N and mMeshTyingMortarCondition3D4N3N are constructed with matching geometries (TriangleType/TriangleType and QuadrilateralType/QuadrilateralType, contact_structural_mechanics_application.cpp:34-35) while every contact prototype of the same suffix uses the true mixed pair (for example mALMFrictionlessMortarContactCondition3D3N4N at line 49 uses TriangleType/QuadrilateralType). The prototype geometry is only a placeholder that Create(NewId, pGeom, pProperties, pMasterGeom) replaces with the real faces, so the registered condition works; the inconsistency only affects code that inspects the prototype itself.

MeshTyingMortarCondition

MeshTyingMortarCondition<TDim, TNumNodes, TNumNodesMaster> (mesh_tying_mortar_condition.{h,cpp}) glues two non-matching meshes with dual Lagrange multipliers (Mesh tying). Its differences with respect to the contact conditions are:

• Generic in the tied variable. Initialize (mesh_tying_mortar_condition.cpp:72-100) reads the TYING_VARIABLE property (default "DISPLACEMENT"). If the name is a Variable<double>, mpDoFVariables = {var} and mpLMVariables = {SCALAR_LAGRANGE_MULTIPLIER}; if it is a Variable<array_1d<double,3>>, the _X, _Y (and _Z in 3D) components are paired with VECTOR_LAGRANGE_MULTIPLIER_X/Y/Z; anything else raises Compatible variables are: double or array_1d<double, 3>. The local size is therefore dynamic, dof_size * (2 * TNumNodes + TNumNodesMaster) with dof_size = mpDoFVariables.size() (mesh_tying_mortar_condition.cpp:260, 652), and there is no MatrixSize constant.
• Operators computed once and cached. The rest of Initialize (l. 101-194) performs the same segmentation, Ae and Gauss loops as CalculateConditionSystem but without derivatives (the operators of a tied interface do not change with the deformation in the small-deformation setting the condition targets) and stores the result in mMortarConditionMatrices. A pair with an empty intersection is switched off with this->Set(ACTIVE, false) (l. 189).
• Hand-written LHS/RHS (no AD). CalculateConditionSystem (l. 362-395) only fills a DofData struct with the current u1, u2 and LagrangeMultipliers (InitializeDofData / UpdateMasterPair) and calls CalculateLocalLHS (l. 528-580) and CalculateLocalRHS (l. 585-640). With \(k\) the scale factor (SCALE_FACTOR if present in the ProcessInfo, else BUILD_SCALE_FACTOR, else 1), the block structure in the order [master, slave, multipliers] is

$$ \mathbf{K}_{tying} = \begin{bmatrix} \mathbf{0} & \mathbf{0} & -k\,\mathbf{M}^T \\ \mathbf{0} & \mathbf{0} & k\,\mathbf{D}^T \\ -k\,\mathbf{M} & k\,\mathbf{D} & \mathbf{0} \end{bmatrix}, \qquad \mathbf{r}_{tying} = \begin{bmatrix} k\,\mathbf{M}^T \boldsymbol{\lambda} \\ -k\,\mathbf{D}^T \boldsymbol{\lambda} \\ -k\,(\mathbf{D}\mathbf{u}^{(1)} - \mathbf{M}\mathbf{u}^{(2)}) \end{bmatrix}, $$

each scalar entry being repeated for the dof_size components of the node (rLocalLHS(row + j*dof_size + i, col + k*dof_size + i), mesh_tying_mortar_condition.cpp:556-575). This is the discrete form of the tying constraint \(\mathbf{D}\mathbf{u}^{(1)} = \mathbf{M}\mathbf{u}^{(2)}\) of the LaTeX derivation in automatic_differentiation/mesh_tying_mortar_condition/mesh_tying_mortar_condition.tex; the automatic_differentiation/mesh_tying_mortar_condition/ generator is legacy (its README states that AD was removed from mesh tying).

• Registered for five pairs (mesh_tying_mortar_condition.cpp:863-867) and driven by MeshTyingProcess (python_scripts/mesh_tying_process.py), which sets TYING_VARIABLE on the pair properties and flags the model part RIGID. The only condition-level C++ test of the application, tests/cpp_tests/conditions/test_mesh_tying_condition.cpp, exercises this class.

MPCMortarContactCondition

MPCMortarContactCondition<TDim, TNumNodes, TNumNodesMaster> (mpc_mortar_contact_condition.{h,cpp}) implements the alternative multipoint-constraint route: the contact constraint is imposed by a ContactMasterSlaveConstraint that ties the slave displacements to the master ones, and the condition exists only to compute and refresh that constraint (Frictional laws and MPC constraint). Consequently:

• MatrixSize = TDim * (TNumNodes + TNumNodesMaster) (mpc_mortar_contact_condition.h:142) and there are no multiplier DoFs. CalculateLeftHandSide returns a zero matrix (mpc_mortar_contact_condition.cpp:351-362). CalculateRightHandSide returns a zero vector except when the condition Is(SLIP): then, for every slave node that Is(SLIP), it adds the explicit Coulomb force \(-\mu\, (\mathbf{R}\cdot\mathbf{n})\, \mathbf{t}\) built from the nodal REACTION, NORMAL, FRICTION_COEFFICIENT and WEIGHTED_SLIP (l. 367-410, with // TODO: Check sign comments).
• The real work happens in the life-cycle hooks. Initialize (l. 78-120), InitializeSolutionStep (l. 125-172, skipped when the condition Is(BLOCKED)) and InitializeNonLinearIteration (l. 178-216, executed only when the condition Is(INTERACTION), which the MPC strategy sets when update_each_nl_iteration is true) share the same body:

if Has(CONSTRAINT_POINTER) and p_const != nullptr:
    dual_LM = MortarExplicitContributionUtilities<..., FRICTIONLESS_PENALTY, false, ...>::ComputePreviousMortarOperators(this, pinfo, D_M, order, false)
    relation_matrix = Zero(TDim*TNumNodes, TDim*TNumNodesMaster); constant_vector = Zero(TDim*TNumNodes)
    if Is(SLIP):     UpdateConstraintFrictional(D_M, relation_matrix, constant_vector, pinfo, dual_LM)     l. 716
    elif Is(RIGID):  UpdateConstraintTying(D_M, relation_matrix, constant_vector, pinfo, dual_LM)          l. 832
    else:            UpdateConstraintFrictionless(D_M, relation_matrix, constant_vector, pinfo, dual_LM)   l. 593
    ConstraintDofDatabaseUpdate(relation_matrix, constant_vector, pinfo)                                    l. 466
    p_const->SetLocalSystem(relation_matrix, constant_vector, pinfo)

• The relation matrix is built from \(\mathbf{D}^{-1}\mathbf{M}\). All three UpdateConstraint* methods invert \(\mathbf{D}\) (trivially, entry by entry, when dual_LM is true because the dual multipliers make \(\mathbf{D}\) diagonal; with MathUtils::InvertMatrix otherwise) and form D_inv_M = prod(inverse_DOperator, MOperator); rows with a negative entry below \(-10^{-8}\) are zeroed (l. 632-641). Then, for every slave node that IsActive(), with weight_coeff = 1 / node.GetValue(NODAL_PAUX):
  • frictionless (l. 593-713): only the normal component is constrained, relation(i,dim_a; j,dim_b) = weight_coeff * n_a n_b * D_inv_M(i,j), and the constant vector carries the gap, constant(i,dim) = weight_coeff * (WEIGHTED_GAP_i / NODAL_AREA_i + n_i . (u1_0 - D_inv_M u2_0)_i) * n_dim (the second term re-expresses the “initial” gap through the previous displacements);
  • frictional (l. 716-829): slip nodes get the normal-only relation, stick nodes get the full vector relation relation(i,dim; j,dim) = weight_coeff * D_inv_M(i,j) (both // TODO comments about nodal forces and slip in the constant vector are still present);
  • tying (l. 832-932): every active node gets the full vector relation, that is the mesh-tying constraint \(\mathbf{u}^{(1)} = \mathbf{D}^{-1}\mathbf{M}\mathbf{u}^{(2)}\) through constraints instead of multipliers.
• ConstraintDofDatabaseUpdate (l. 466-590) prunes rows and columns whose absolute sum is below \(10^{-4}\), rebuilds the reduced relation matrix and constant vector and calls p_const->SetDofList(...) with the surviving slave (DISPLACEMENT_X/Y/Z of active slave nodes) and master DoFs, so the constraint only involves the nodes that are really in contact.
• The frictional MPC also stores mPreviousMortarOperators (mpc_mortar_contact_condition.h:409-411), initialized in Initialize when Is(SLIP) and refreshed once per step by ComputePreviousMortarOperators (l. 935-943). AddExplicitContribution(const ProcessInfo&) (l. 447-462) updates WEIGHTED_GAP (and WEIGHTED_SLIP for SLIP) exactly like the penalty conditions, which is what MPCContactCriteria and the search read.

The CONSTRAINT_POINTER value is set by MPCContactSearchProcess when it creates the pair together with its ContactMasterSlaveConstraint; the SLIP / RIGID flags on the condition are inherited from the model part flags set by MPCContactProcess (mpc_contact_process.py:158-162). The dispatch flow is drawn in Frictional laws and MPC constraint.

What the conditions read from Properties and ProcessInfo

The conditions are configured indirectly: the Python processes push values into the pair Properties, into the ProcessInfo and into nodal non-historical values, and the conditions read them at assembly time. The full list of variables is in the Variables and flags reference; the ones that reach the conditions are:

Consistency checks and integration-point output

Check() (mortar_contact_condition.cpp:661-690 plus the family overrides) verifies at start-up that the slave nodes carry DISPLACEMENT, WEIGHTED_GAP, NORMAL and the displacement DoFs, that the pair geometry has been initialized (YOU HAVE NOT INITIALIZED THE PAIR GEOMETRY IN THE MortarContactCondition) and, per family, the multiplier variable and DoFs (VECTOR_LAGRANGE_MULTIPLIER and WEIGHTED_SLIP for the frictional condition, ALM_frictional_mortar_contact_condition_template.cpp:268-290). The solver adds these variables and DoFs according to mortar_type (Variables and flags reference).

The three CalculateOnIntegrationPoints overloads (double, array_1d<double,3> and Vector; mortar_contact_condition.cpp:592-658, mesh_tying_mortar_condition.cpp:746-818, mpc_mortar_contact_condition.cpp:945-1011) resize the output to the number of integration points of the slave geometry and fill it with zeros: the contact conditions do not carry integration-point results. Everything a user wants to post-process (multipliers, augmented pressures, weighted gap and slip, normals, flags) is stored on the slave nodes, see Output and post-processing. The Info() / PrintInfo() / PrintData() methods print the class name, the id and the two geometries (paired_condition.h:306-325).

Adding a new formulation family

The structure above makes the list of touch points for a new mortar family (or a new geometry pair) explicit. All of them must be kept consistent, otherwise the local sizes of the generated code and of the DoF lists disagree at run time:

1. Add the enumerator to FrictionalCase in contact_structural_mechanics_application_variables.h and extend the MatrixSize / IsFrictional / DerivativeDataType selectors of mortar_contact_condition.h:140-159.
2. Write the family header (copying ALM_frictionless_mortar_contact_condition.h is the simplest start): MatrixSize, EquationIdVector / GetDofList in the [master, slave, multiplier] order, GetActiveInactiveValue, Check, the Create overloads, AddExplicitContribution and any extra state (with save/load).
3. Write the sympy generator and the *_template.cpp in a new automatic_differentiation/<family>/ folder (Automatic differentiation); run it with sympy 1.2 and copy the output to custom_conditions/<family>.cpp. The explicit instantiations for the five geometry pairs times the two normal-variation values go at the bottom of the template.
4. Add the explicit template class MortarContactCondition<..., NEW_CASE, ...> lines to mortar_contact_condition.cpp (and, if the family needs them, to custom_frictional_laws/*.cpp, custom_utilities/derivatives_utilities.cpp and custom_utilities/mortar_explicit_contribution_utilities.cpp, which are instantiated per FrictionalCase).
5. Declare the prototypes in contact_structural_mechanics_application.h, build them in the constructor and register them in Register() following the naming decoder, so that BaseContactSearchProcess can find <prefix>Condition<dim>D<n>N[<m>N].
6. Teach the Python side: _get_condition_name of the corresponding process (alm_contact_process.py, penalty_contact_process.py), the mortar_type branches of AuxiliaryAddVariables / AuxiliaryAddDofs (auxiliary_methods_solvers.py:122-170) and, if a new multiplier variable is introduced, the MixedULMLinearSolver block classification and the convergence criteria (Strategies and convergence criteria).
7. Add a case to tests/cpp_tests/ or to the Python test factory (Test suite reference); the patch tests of tests/ALM_frictionless_contact_test_2D/ are the quickest regression check.

A new geometry pair (for instance quadratic faces) is heavier: it needs a PointBelongs* hash in kratos/includes/mortar_classes.h, support in ExactMortarIntegrationUtility, MortarUtilities and DerivativesUtilities, one more dim_combinations / nnodes_combinations entry in every generator, and the corresponding prototypes; the application currently supports linear geometries only.

Where to look next

Source files on GitHub: custom_conditions/, mortar_contact_condition.cpp, paired_condition.h, mesh_tying_mortar_condition.cpp, mpc_mortar_contact_condition.cpp, contact_structural_mechanics_application.cpp, kratos/includes/mortar_classes.h.