get_cell_dof_ids skips non-zero offset for field 1 in BlockMultiFieldStyle with permuted blocks, causing silent DOF corruption

[anusha19murthy]

Describe the bug

get_cell_dof_ids in MultiFieldFESpaces.jl incorrectly skips applying the DOF offset for field 1 when using BlockMultiFieldStyle with a non-identity permutation P.

The implementation assumes that field 1 always has zero offset:

if i == 1
push!(active_block_data,cell_dofs_i)

This assumption is valid for:

ConsecutiveMultiFieldStyle

BlockMultiFieldStyle with identity permutation

However, for BlockMultiFieldStyle with a non-identity permutation, field 1 may not be first within its block, and therefore offsets[1] may be non-zero.

Since compute_field_offsets correctly computes within-block offsets accounting for permutation, skipping the offset application for i == 1 causes incorrect DOF numbering.

As a consequence, DOF slices belonging to different fields may overlap or be swapped.

Impact

Field offsets are later used to extract each field’s DOF slice from the global block vector.

When offsets are incorrect:

• fields may read DOF ranges belonging to other fields
• DOF slices can overlap or become swapped
• matrix and vector assembly complete without errors
• linear and nonlinear solvers converge normally
• resulting solutions silently contain corrupted or swapped field values

This produces incorrect numerical results in multiphysics simulations using BlockMultiFieldStyle with explicit permutations.

The failure is silent and difficult to detect because:

• no bounds errors occur
• matrix dimensions are correct
• solver convergence is unaffected

This can silently corrupt multiphysics block-preconditioned simulations (e.g., Stokes, Darcy, coupled problems).

Affected file
src/MultiField/MultiFieldFESpaces.jl

Function:

FESpaces.get_cell_dof_ids(
  f::MultiFieldFESpace,
  trian::Triangulation,
  ::Union{<:ConsecutiveMultiFieldStyle,<:BlockMultiFieldStyle}
)

Related:

compute_field_offsets(
  f::MultiFieldFESpace,
  ::BlockMultiFieldStyle{NB,SB,P}
)

Exact problematic code

for i in active_block_ids
  cell_dofs_i = get_cell_dof_ids(f.spaces[i],trian)
  if i == 1
    push!(active_block_data,cell_dofs_i)
  else
    offset = Int32(offsets[i])
    o = Fill(offset,length(cell_dofs_i))
    cell_dofs_i_b =
      lazy_map(Broadcasting(_sum_if_first_positive),cell_dofs_i,o)
    push!(active_block_data,cell_dofs_i_b)
  end
end

The condition if i == 1 assumes offsets[1] == 0, which is not guaranteed for non-identity permutations.

Possible fix

Apply the offset based on its value rather than the field index.

Replace:

if i == 1

with:

if offsets[i] == 0

Full corrected loop:

for i in active_block_ids
  cell_dofs_i = get_cell_dof_ids(f.spaces[i],trian)
  if offsets[i] == 0
    push!(active_block_data,cell_dofs_i)
  else
    offset = Int32(offsets[i])
    o = Fill(offset,length(cell_dofs_i))
    cell_dofs_i_b =
      lazy_map(Broadcasting(_sum_if_first_positive),cell_dofs_i,o)
    push!(active_block_data,cell_dofs_i_b)
  end
end

This ensures offsets are respected for all field orderings and prevents silent DOF misassignment.