Summary

ResidualBasedEliminationBuilderAndSolverWithConstraints::AssembleRHSWithoutConstraints
(in kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h)
contains a heap-buffer-overflow that triggers whenever all three of the following are true
in the same analysis:

1. The model has one or more MasterSlaveConstraints.
2. At least one DOF is genuinely fixed (Dof::IsFixed() == true).
3. mCalculateReactionsFlag is true (i.e. reactions/REACTION are being computed).

Confirmed with AddressSanitizer:

==305058==ERROR: AddressSanitizer: heap-buffer-overflow on address 0x7be9f77bae88
READ of size 8 at 0x7be9f77bae88 thread T13
    #0 in void Kratos::AtomicAdd<double>(double&, double const&)
       kratos/utilities/atomic_utilities.h:60
    #1 in Kratos::ResidualBasedEliminationBuilderAndSolverWithConstraints<...>::AssembleRHSWithoutConstraints(...)
       kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h:1948

In a real project this manifests as a hard crash (SIGABRT / double free or corruption (!prev) / free(): invalid next size under plain glibc malloc, depending on how the
corrupted memory is later reused) rather than a clean exception — it can appear far from
the actual bug (e.g. inside KratosContactStructuralMechanicsApplication's contact search,
one or more solve steps later), which made it very difficult to trace without ASan.

Root cause

In AssembleRHSWithoutConstraints, the mCalculateReactionsFlag == true branch does:

} else {
    TSystemVectorType& r_reactions_vector = *BaseType::mpReactionsVector;
    for (IndexType i_local = 0; i_local < local_size; ++i_local) {
        const IndexType i_global = rEquationId[i_local];
        auto it_dof = BaseType::mDofSet.begin() + i_global;   // <-- BUG

        const bool is_master_fixed = mDoFMasterFixedSet.find(*it_dof) == mDoFMasterFixedSet.end() ? false : true;
        const bool is_slave = mDoFSlaveSet.find(*it_dof) == mDoFSlaveSet.end() ? false : true;
        if (is_master_fixed || is_slave) { // Fixed or MPC dof
            double& r_b_value = r_reactions_vector[mReactionEquationIdMap[i_global]];
            ...

BaseType::mDofSet.begin() + i_global assumes a DOF's position in the mDofSet
container equals its equation ID (i_global). That assumption is false in this class:

• SetUpDofSetWithConstraints builds mDofSet by collecting all DOFs into a set, then
  calls dof_temp_all.Sort() (sorts by DOF identity — node id + variable) and assigns
  BaseType::mDofSet = dof_temp_all.
• SetUpSystemWithConstraints (called afterwards) assigns equation IDs by iterating
  BaseType::mDofSet in that same identity-sorted order, but with two independent
  counters moving in opposite directions:

  int free_id = 0;
  int fix_id = BaseType::mDofSet.size();
  for (auto& dof : BaseType::mDofSet) {
      if (dof.IsFixed()) {
          auto it = mDoFMasterFixedSet.find(dof);
          if (it == mDoFMasterFixedSet.end()) {
              dof.SetEquationId(--fix_id);
          } else {
              dof.SetEquationId(free_id++);
          }
      } else {
          dof.SetEquationId(free_id++);
      }
  }

• The resulting equation ID for a given DOF has no fixed relationship to that DOF's
  position in the mDofSet container — it depends on how many free vs. fixed DOFs were
  seen before it in the iteration, which is unrelated to mDofSet's storage order.

So mDofSet.begin() + i_global fetches an essentially arbitrary, usually incorrect Dof.
This corrupts:

• the is_master_fixed / is_slave classification of the DOF actually being processed,
  and
• consequently the index passed to r_reactions_vector[mReactionEquationIdMap[i_global]]
  can be wildly out of range (we observed indices corresponding to reads ~2400 bytes
  before a ~4.4MB heap allocation, i.e. a large/negative-wrapped IndexType).

Notably, elsewhere in the very same file (FinalizeAndClearProcessInfo /
FinalizeSolutionStep-adjacent reaction-accumulation code, ~line 556), the correct pattern
is already used — a direct range-based iteration over BaseType::mDofSet, which never
suffers this issue:

for (auto& r_dof : BaseType::mDofSet) {
    const bool is_master_fixed = mDoFMasterFixedSet.find(r_dof) == mDoFMasterFixedSet.end() ? false : true;
    const bool is_slave = mDoFSlaveSet.find(r_dof) == mDoFSlaveSet.end() ? false : true;
    if (is_master_fixed || is_slave) {
        const IndexType equation_id = r_dof.EquationId();
        r_reactions_vector[mReactionEquationIdMap[equation_id]] += rb[equation_id];
    }
}

AssembleRHSWithoutConstraints, however, is called per-element/per-condition from inside
Build/BuildAndSolve's parallel assembly loop, where only rEquationId (a plain vector
of equation IDs for the current element's DOFs) is available — not direct Dof pointers —
which is presumably why the (incorrect) reverse lookup via mDofSet.begin() + i_global was
attempted instead.

Suggested fix

Build a lookup keyed by equation ID (not container position) at the same point
mReactionEquationIdMap is populated in SetUpSystemWithConstraints, and use that in
AssembleRHSWithoutConstraints instead of the position-based mDofSet.begin() + i_global.

Patch (tested locally, see Verification below):

// New member, alongside mReactionEquationIdMap:
IndexMapType mReactionCategoryByEquationId; // equation_id -> {1: is_master_fixed, 2: is_slave}

// In SetUpSystemWithConstraints, where mReactionEquationIdMap is currently populated:
counter = 0;
mReactionCategoryByEquationId.clear();
for (auto& r_dof : BaseType::mDofSet) {
    const bool is_master_fixed = mDoFMasterFixedSet.find(r_dof) == mDoFMasterFixedSet.end() ? false : true;
    const bool is_slave = mDoFSlaveSet.find(r_dof) == mDoFSlaveSet.end() ? false : true;
    if (is_master_fixed || is_slave) {
        mReactionEquationIdMap.insert({r_dof.EquationId(), counter});
        mReactionCategoryByEquationId.insert({r_dof.EquationId(), is_master_fixed ? 1 : 2});
        ++counter;
    }
}

// In AssembleRHSWithoutConstraints, mCalculateReactionsFlag == true branch:
TSystemVectorType& r_reactions_vector = *BaseType::mpReactionsVector;
for (IndexType i_local = 0; i_local < local_size; ++i_local) {
    const IndexType i_global = rEquationId[i_local];

    auto it_category = mReactionCategoryByEquationId.find(i_global);
    if (it_category != mReactionCategoryByEquationId.end()) { // Fixed or MPC dof
        double& r_b_value = r_reactions_vector[mReactionEquationIdMap[i_global]];
        const double rhs_value = rRHSContribution[i_local];
        AtomicAdd(r_b_value, rhs_value);
    } else if (i_global < BaseType::mEquationSystemSize) {  // Free dof not in the MPC
        double& r_b_value = rb[i_global];
        const double& rhs_value = rRHSContribution[i_local];
        AtomicAdd(r_b_value, rhs_value);
    }
}

Also clear mReactionCategoryByEquationId in Clear() alongside mReactionEquationIdMap.

Steps to reproduce

1. Build a model with:
   • At least one MasterSlaveConstraint (e.g. a tie constraint between two node groups).
   • At least one node group with a genuinely fixed DOF (e.g. via
     AssignVectorVariableProcess/AssignScalarVariableProcess with constrained: true,
     the default) — the fixed group should be reasonably large (thousands of DOFs) to make
     the corrupted index likely to land somewhere that crashes rather than silently
     succeeding.
   • calculate_reactions: true (or whatever the solver setting is that sets
     mCalculateReactionsFlag), needed to output REACTION.
   • builder_and_solver_settings.type: "elimination".
2. Run a non-linear structural analysis (we reproduced with
   ResidualBasedNewtonRaphsonContactStrategy +
   ResidualBasedEliminationBuilderAndSolverWithConstraints, but the bug is in the
   builder/solver, not the contact strategy — likely reproducible without contact too).
3. Build with -fsanitize=address to see the ASan report directly, or run under plain
   glibc and observe a crash (double free or corruption / SIGABRT) a short time after
   the first step where the newly-fixed DOFs are actually processed.

Environment

• Kratos version: source checkout at commit 6c03e65759 (tag-ish description
  DGeoFlow-2026.01-302-g6c03e65759), applications built: LinearSolversApplication,
  StructuralMechanicsApplication, ConstitutiveLawsApplication,
  ContactStructuralMechanicsApplication.
• Also confirmed present in the pip-distributed wheel KratosMultiphysics==10.4.2 (same
  crash reproduced against the prebuilt wheel before rebuilding from source to debug it).
• OS: Fedora 44, GCC 16.1.1, Python 3.14.6, x86_64.
• Compiler flags for the confirming ASan build: -fsanitize=address -fno-omit-frame-pointer -g -std=c++20, CMAKE_BUILD_TYPE=Debug.

Verification

Applied the suggested fix locally and confirmed it two ways:

1. Debug + ASan build: reran the originally-failing case; the heap-buffer-overflow no
   longer occurs.
2. Plain Release build (no sanitizer, normal optimized build a real user would run):
   before the fix, the same test case reliably crashed with double free or corruption (!prev) / free(): invalid next size at the exact same point in every one of 5+
   isolated repro runs (varying whether the fix-on-DOF was applied via the JSON
   assign_*_process, directly via Python Dof::Fix(), or with/without other unrelated
   changes present) — always right at the first assembly of reactions for the newly-fixed
   DOF group, one step after the initial factorization. After the fix, the identical test
   case (large MasterSlaveConstraint-tied model, ~3600-node fixed boundary condition,
   calculate_reactions: true) ran multiple time steps to completion (Analysis -END-)
   without any crash.

Suggested patch

diff --git a/kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h b/kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h
index f4e45f848f..d8612a445d 100644
--- a/kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h
+++ b/kratos/solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_with_constraints.h
@@ -389,6 +389,13 @@ protected:
     DofsArrayType mDoFSlaveSet;                            /// The set containing the slave DoF of the system
     SizeType mDoFToSolveSystemSize = 0;                    /// Number of degrees of freedom of the problem to actually be solved
     IndexMapType mReactionEquationIdMap;                   /// In order to know the corresponding EquaionId for each component of the reaction vector
+    IndexMapType mReactionCategoryByEquationId;            /// BUGFIX: maps equation_id -> {0: neither, 1: is_master_fixed, 2: is_slave} so
+                                                            /// AssembleRHSWithoutConstraints can classify a DOF by its equation_id directly,
+                                                            /// without the incorrect "BaseType::mDofSet.begin() + i_global" pointer arithmetic
+                                                            /// (mDofSet is ordered by DOF identity, NOT by equation_id, so that lookup fetched
+                                                            /// an arbitrary wrong DOF and corrupted the reactions-vector index, causing a
+                                                            /// heap-buffer-overflow whenever constraints + genuinely fixed DOFs + reactions
+                                                            /// were all in use simultaneously).
 
     bool mCheckConstraintRelation = false;                 /// If we do a constraint check relation
     bool mResetRelationMatrixEachIteration = false;        /// If we reset the relation matrix at each iteration
@@ -1434,6 +1441,7 @@ protected:
 
         // Clearing the relation map
         mReactionEquationIdMap.clear();
+        mReactionCategoryByEquationId.clear();
 
         // Clear constraint system
         if (mpTMatrix != nullptr)
@@ -1593,11 +1601,16 @@ private:
 
         // Finally we build the relation between the EquationID and the component of the reaction
         counter = 0;
+        mReactionCategoryByEquationId.clear();
         for (auto& r_dof : BaseType::mDofSet) {
             const bool is_master_fixed = mDoFMasterFixedSet.find(r_dof) == mDoFMasterFixedSet.end() ? false : true;
             const bool is_slave = mDoFSlaveSet.find(r_dof) == mDoFSlaveSet.end() ? false : true;
             if (is_master_fixed || is_slave) { // Fixed or MPC dof
                 mReactionEquationIdMap.insert({r_dof.EquationId(), counter});
+                // BUGFIX: record the classification keyed by equation_id (not container
+                // position), so AssembleRHSWithoutConstraints can look it up correctly
+                // instead of using the incorrect "mDofSet.begin() + i_global" arithmetic.
+                mReactionCategoryByEquationId.insert({r_dof.EquationId(), is_master_fixed ? 1 : 2});
                 ++counter;
             }
         }
@@ -1931,19 +1944,31 @@ private:
                 }
             }
         } else {
+            // BUGFIX: this branch used to look up the Dof via
+            // "BaseType::mDofSet.begin() + i_global", treating the equation_id as if it
+            // were a container position in mDofSet. mDofSet is sorted by DOF identity
+            // (see SetUpDofSetWithConstraints), while equation IDs are assigned via two
+            // independent counters (ascending for free DOFs, descending for master-fixed
+            // ones) in that same identity-sorted iteration order — so a DOF's equation_id
+            // has no fixed relationship to its position in mDofSet. The old code therefore
+            // fetched an arbitrary, usually wrong Dof, corrupting the is_master_fixed/
+            // is_slave classification and the mReactionEquationIdMap lookup, causing a
+            // heap-buffer-overflow (out-of-bounds read/write into the reactions vector)
+            // whenever MasterSlaveConstraints + genuinely fixed DOFs + reactions were all
+            // in use at once. Fixed by classifying purely via equation_id, using the
+            // mReactionCategoryByEquationId map built alongside mReactionEquationIdMap in
+            // SetUpSystemWithConstraints (which iterates real Dof objects correctly).
             TSystemVectorType& r_reactions_vector = *BaseType::mpReactionsVector;
             for (IndexType i_local = 0; i_local < local_size; ++i_local) {
                 const IndexType i_global = rEquationId[i_local];
-                auto it_dof = BaseType::mDofSet.begin() + i_global;
 
-                const bool is_master_fixed = mDoFMasterFixedSet.find(*it_dof) == mDoFMasterFixedSet.end() ? false : true;
-                const bool is_slave = mDoFSlaveSet.find(*it_dof) == mDoFSlaveSet.end() ? false : true;
-                if (is_master_fixed || is_slave) { // Fixed or MPC dof
+                auto it_category = mReactionCategoryByEquationId.find(i_global);
+                if (it_category != mReactionCategoryByEquationId.end()) { // Fixed or MPC dof
                     double& r_b_value = r_reactions_vector[mReactionEquationIdMap[i_global]];
                     const double rhs_value = rRHSContribution[i_local];
 
                     AtomicAdd(r_b_value, rhs_value);
-                } else if (it_dof->IsFree()) {  // Free dof not in the MPC
+                } else if (i_global < BaseType::mEquationSystemSize) {  // Free dof not in the MPC
                     // ASSEMBLING THE SYSTEM VECTOR
                     double& r_b_value = rb[i_global];
                     const double& rhs_value = rRHSContribution[i_local];