gridap / gridap/Gridap.jl

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

Open
#1,233 0 comments 0 reactions 0 assignees View on GitHub
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

Open the contributing 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.