Ferrite-FEM / Ferrite-FEM/Ferrite.jl
:robot: `keep_constrained = false` misses sparsity entries needed by `apply_assemble!` when an affine master is in a non-adjacent cell
- Dominant language
- Julia
- Stars
- 453
- Forks
- 115
- Avg merge
- 1d 17h
- Merged PRs (30d)
- 11
Description
With `keep_constrained = false`, `add_cell_entries!` skips all entries touching
constrained dofs, and `add_constraint_entries!` only expands entries that are already
stored in the pattern. The eliminated entries are therefore invisible to the expansion,
and the master-dof entries that local condensation writes to are never created when the
master sits in a cell that does not share the relevant dofs. `apply_assemble!` then fails
with `SparsityError`.
MWE (Ferrite v1.6.0, also reproduces on current master): 1D Poisson on 3 line cells with
a periodic-style constraint u(left) = u(right). Condensing the slave row of cell 1 writes
to `K[master, j]` for the other dof `j` of cell 1, but no cell contains both `master` and
`j`, and the eliminated slave entries `(slave, j)` that would have been expanded to
`(master, j)` were never inserted:
```julia
using Ferrite
grid = generate_grid(Line, (3,))
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefLine, 1}())
close!(dh)
# Slave: the boundary dof of cell 1; master: the boundary dof of cell 3
slave = only(setdiff(celldofs(dh, 1), celldofs(dh, 2)))
master = only(setdiff(celldofs(dh, 3), celldofs(dh, 2)))
ch = ConstraintHandler(dh)
add!(ch, AffineConstraint(slave, [master => 1.0], 0.0))
close!(ch)
K = allocate_matrix(dh, ch; keep_constrained = false)
f = zeros(ndofs(dh))
cv = CellValues(QuadratureRule{RefLine}(2), Lagrange{RefLine, 1}())
assembler = start_assemble(K, f)
Ke = zeros(2, 2)
fe = zeros(2)
for cell in CellIterator(dh)
reinit!(cv, cell)
fill!(Ke, 0)
fill!(fe, 0)
for qp in 1:getnquadpoints(cv)
dΩ = getdetJdV(cv, qp)
for i in 1:2, j in 1:2
Ke[i, j] += shape_gradient(cv, qp, i) ⋅ shape_gradient(cv, qp, j) * dΩ
end
for i in 1:2
fe[i] += shape_value(cv, qp, i) * dΩ
end
end
apply_assemble!(assembler, ch, celldofs(cell), Ke, fe)
end
```
```
ERROR: SparsityError: writing to an index outside the sparsity pattern is not allowed
Stacktrace:
[1] addindex!
[2] _condense_local!(...)
[3] _apply_local!(...)
[4] apply_assemble!(...)
```
The identical loop works with `keep_constrained = true`: there the slave entries are
stored, `add_constraint_entries!` expands them to the master, and condensation finds its
target entries.
In 2D/3D the same happens for e.g. `PeriodicDirichlet`, where the image dofs' cells
generally do not contain all dofs of the mirror dofs' cells — so
`keep_constrained = false` (whose docs point to condensed assembly as the intended usage)
currently cannot be combined with cross-cell affine constraints at all.
A possible fix is to distribute eliminated entries through the constraint coefficients at
insertion time in `_add_cell_entries!`: for an entry `(i, j)` where `i` and/or `j` is
constrained, add the entries `(i′, j′)` for the (unconstrained) master dofs `i′`/`j′`
instead of dropping the entry, and drop only entries of masterless prescribed dofs. These
are exactly the entries `_condense_local!` writes to. Note that #1397 reimplements the
pattern build (preserving the current semantics), so a fix should probably be coordinated
with it.
Contributor guide
Assessment
This issue has not been assessed yet.