get_cell_dof_ids skips non-zero offset for field 1 in BlockMultiFieldStyle with permuted blocks, causing silent DOF corruption
- Dominant language
- Julia
- Stars
- 879
- Forks
- 119
- Avg merge
- 6d 3h
- Merged PRs (30d)
- 4
Description
**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.
Contributor guide
Research direction
Start in src/MultiField/MultiFieldFESpaces.jl at FESpaces.get_cell_dof_ids and compare its offset handling with compute_field_offsets for BlockMultiFieldStyle. Reproduce the case with a non-identity permutation and verify that field 1 receives its non-zero offset without overlapping or swapping DOF slices.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- backend
- Issue type
- Bug
- Difficulty
- 2/5
- Estimated time
- 1-3 hours
- Activity status
- Stale
- Clarity
- Clearly specified
- Newbie friendliness
- 58/100