Skip to content

Commit

Permalink
Handle prescribed dofs in RHS of affine constraints
Browse files Browse the repository at this point in the history 
This patch fixes affine constraints with prescribed dofs in the RHS. In
particular, we allow dofs that are prescribed by just an inhomogeneity
(i.e. DBC) but disallow "nesting" affine constraints.

Concretely, consider e.g. the following two constraints:

    u2 = f(t)
    u3 = u2 + b3

Before this patch this was not handled correctly since the inhomogeneity
for u3 was taken as b3, but it should really be b3 + f(t) by
substituting u2 for f(t). Since we allow for time-dependent
inhomogeneities this substitution can not be done in
close!(::ConstraintHandler) and instead the effective inhomogeneity is
computed in update!.

Nested constraints, e.g.

    u2 = u3
    u3 = u5

are still not allowed but in the future this can be resolved in
close!(::ConstraintHandler) to

    u2 = u5
    u3 = u5

However, this patch checks for such nesting and raises an error instead
of resulting in incorrect answers as is the case before.

Fixes #530.
  • Loading branch information
@fredrikekre
fredrikekre committed Nov 9, 2022
1 parent dd88d1e commit d68f9a1
Show file tree
Hide file tree
Showing 4 changed files with 191 additions and 5 deletions.
7 changes: 7 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
always override any previous constraints. Conflicting constraints could previously cause
problems when a DoF where prescribed by both `Dirichlet` and `AffineConstraint`.
([#529][github-529])
### Fixed
- Fix affine constraints with prescribed dofs in the right-hand-side. In particular, dofs
that are prescribed by just an inhomogeneity are now handled correctly, and nested affine
constraints now give an error instead of silently giving the wrong result.
([#530][github-530], [#535][github-535])

## [0.3.9] - 2022-10-19
### Added
Expand Down Expand Up @@ -153,6 +158,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
[github-514]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/514
[github-524]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/524
[github-529]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/529
[github-530]: https://github.com/Ferrite-FEM/Ferrite.jl/issues/530
[github-535]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/535

[Unreleased]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.9...HEAD
[0.3.9]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.8...v0.3.9
Expand Down
4 changes: 2 additions & 2 deletions docs/src/literate/stokes-flow.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,7 +319,7 @@ function setup_mean_constraint(dh, fvp)
V ./= V[constrained_dof]
mean_value_constraint = AffineConstraint(
constrained_dof,
Pair{Int,Float64}[J[i] => -V[i] for i in 1:length(J) if i != constrained_dof],
Pair{Int,Float64}[J[i] => -V[i] for i in 1:length(J) if J[i] != constrained_dof],
0.0,
)
return mean_value_constraint
Expand Down Expand Up @@ -460,7 +460,7 @@ function main()
vtk_grid("stokes-flow", grid) do vtk
vtk_point_data(vtk, dh, u)
end
Sys.isapple() || @test norm(u) ≈ 0.32353713318639005 #src
Sys.isapple() || @test norm(u) ≈ 0.32254330524111213 #src
return
end
#md nothing #hide
Expand Down
53 changes: 50 additions & 3 deletions src/Dofs/ConstraintHandler.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ struct ConstraintHandler{DH<:AbstractDofHandler,T}
prescribed_dofs::Vector{Int}
free_dofs::Vector{Int}
inhomogeneities::Vector{T}
# Store the original constant inhomogeneities for affine constraints used to compute
# "effective" inhomogeneities in `update!` and then stored in .inhomogeneities.
affine_inhomogeneities::Vector{Union{Nothing,T}}
# `nothing` for pure DBC constraint, otherwise affine constraint
dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}
# global dof -> index into dofs and inhomogeneities and dofcoefficients
Expand All @@ -88,7 +91,7 @@ end
function ConstraintHandler(dh::AbstractDofHandler)
@assert isclosed(dh)
ConstraintHandler(
Dirichlet[], Int[], Int[], Float64[], Union{Nothing,DofCoefficients{Float64}}[],
Dirichlet[], Int[], Int[], Float64[], Union{Nothing,Float64}[], Union{Nothing,DofCoefficients{Float64}}[],
Dict{Int,Int}(), BCValues{Float64}[], dh, ScalarWrapper(false),
)
end
Expand Down Expand Up @@ -181,6 +184,7 @@ function close!(ch::ConstraintHandler)
I = sortperm(ch.prescribed_dofs)
ch.prescribed_dofs .= ch.prescribed_dofs[I]
ch.inhomogeneities .= ch.inhomogeneities[I]
ch.affine_inhomogeneities .= ch.affine_inhomogeneities[I]
ch.dofcoefficients .= ch.dofcoefficients[I]

copy!(ch.free_dofs, setdiff(1:ndofs(ch.dh), ch.prescribed_dofs))
Expand All @@ -198,6 +202,25 @@ function close!(ch::ConstraintHandler)
# such that they become independent. However, at this point, it is left to
# the user to assure this.

# Basic verification of constraints:
# - `add_prescribed_dof` make sure all prescribed dofs are unique by overwriting the old
# constraint when adding a new (TODO: Might change in the future, see comment in
# `add_prescribed_dof`.)
# - We allow affine constraints to have prescribed dofs as master dofs iff those master
# dofs are constrained with just an inhomogeneity (i.e. DBC). The effective
# inhomogeneity is computed in `update!`.
for coeffs in ch.dofcoefficients
coeffs === nothing && continue
for (d, _) in coeffs
i = get(ch.dofmapping, d, 0)
i == 0 && continue
icoeffs = ch.dofcoefficients[i]
if !(icoeffs === nothing || isempty(icoeffs))
error("nested affine constraints currently not supported")
end
end
end

ch.closed[] = true
return ch
end
Expand Down Expand Up @@ -262,11 +285,13 @@ function add_prescribed_dof!(ch::ConstraintHandler, constrained_dof::Int, inhomo
@debug @warn "dof $i already prescribed, overriding the old constraint"
ch.prescribed_dofs[i] = constrained_dof
ch.inhomogeneities[i] = inhomogeneity
ch.affine_inhomogeneities[i] = dofcoefficients === nothing ? nothing : inhomogeneity
ch.dofcoefficients[i] = dofcoefficients
else
N = length(ch.dofmapping)
push!(ch.prescribed_dofs, constrained_dof)
push!(ch.inhomogeneities, inhomogeneity)
push!(ch.affine_inhomogeneities, dofcoefficients === nothing ? nothing : inhomogeneity)
push!(ch.dofcoefficients, dofcoefficients)
ch.dofmapping[constrained_dof] = N + 1
end
Expand Down Expand Up @@ -377,6 +402,26 @@ function update!(ch::ConstraintHandler, time::Real=0.0)
_update!(ch.inhomogeneities, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset,
dbc.components, ch.dh, ch.bcvalues[i], ch.dofmapping, ch.dofcoefficients, convert(Float64, time))
end
# Compute effective inhomogeneity for affine constraints with prescribed dofs in the
# RHS. For example, in u2 = w3 * u3 + w4 * u4 + b2 we allow e.g. u3 to be prescribed by
# a trivial constraint with just an inhomogeneity (e.g. DBC), for example u3 = f(t).
# This value have now been computed in _update! and we can compute the effective
# inhomogeneity h2 for u2 which becomes h2 = w3 * u3 + b2 = w3 * f3(t) + b2.
for i in eachindex(ch.prescribed_dofs, ch.dofcoefficients, ch.inhomogeneities)
coeffs = ch.dofcoefficients[i]
coeffs === nothing && continue
h = ch.affine_inhomogeneities[i]
@assert h !== nothing
for (d, w) in coeffs
j = get(ch.dofmapping, d, 0)
j == 0 && continue
# If this dof is prescribed it must only have an inhomogeneity (verified in close!)
@assert (jcoeffs = ch.dofcoefficients[j]; jcoeffs === nothing || isempty(jcoeffs))
h += ch.inhomogeneities[j] * w
end
ch.inhomogeneities[i] = h
end
return nothing
end

# for faces
Expand Down Expand Up @@ -550,9 +595,11 @@ apply!( v::AbstractVector, ch::ConstraintHandler) = _apply_v(v, ch, false)
function _apply_v(v::AbstractVector, ch::ConstraintHandler, apply_zero::Bool)
@assert length(v) >= ndofs(ch.dh)
v[ch.prescribed_dofs] .= apply_zero ? 0.0 : ch.inhomogeneities
# Apply affine constraints, e.g u2 = s6*u6 + s3*u3 (inhomogenity added above)
for (dof, dofcoef) in zip(ch.prescribed_dofs, ch.dofcoefficients)
# Apply affine constraints, e.g u2 = s6*u6 + s3*u3 + h2
for (dof, dofcoef, h) in zip(ch.prescribed_dofs, ch.dofcoefficients, ch.affine_inhomogeneities)
dofcoef === nothing && continue
@assert h !== nothing
v[dof] = h
for (d, s) in dofcoef
v[dof] += s * v[d]
end
Expand Down
132 changes: 132 additions & 0 deletions test/test_constraints.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1032,5 +1032,137 @@ end # testset
end
end

end # testset

@testset "Affine constraints with master dofs that are prescribed" begin
grid = generate_grid(Quadrilateral, (2, 2))
dh = DofHandler(grid); push!(dh, :u, 1); close!(dh)

# 8───7───9
# │ │ │
# 4───3───6
# │ │ │
# 1───2───5

ke = [ 1.0 -0.25 -0.5 -0.25
-0.25 1.0 -0.25 -0.5
-0.5 -0.25 1.0 -0.25
-0.25 -0.5 -0.25 1.0 ]
fe = rand(4)

@testset "affine constraints before/after Dirichlet" begin
# Construct two ConstraintHandler's which should result in the same end result.

## Ordering of constraints for first ConstraintHandler:
## 1. DBC left: u1 = u4 = u8 = 0
## 2. DBC right: u5 = u6 = u9 = 1
## 3. Periodic bottom/top: u1 = u8, u2 = u7, u5 = u9
## meaning that u1 = 0 and u5 = 1 are overwritten by 3 and we end up with
## u1 = u8 = 0
## u2 = u7
## u4 = 0
## u5 = u9 = 1
## u6 = 1
## u8 = 0
## u9 = 1
## where the inhomogeneity of u1 and u5 have to be resolved at runtime.
ch1 = ConstraintHandler(dh)
add!(ch1, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0))
add!(ch1, Dirichlet(:u, getfaceset(grid, "right"), (x, t) -> 1))
add!(ch1, PeriodicDirichlet(:u, collect_periodic_faces(grid, "bottom", "top")))
close!(ch1)
update!(ch1, 0)

## Ordering of constraints for second ConstraintHandler:
## 1. Periodic bottom/top: u1 = u8, u2 = u7, u5 = u9
## 2. DBC left: u1 = u4 = u8 = 0
## 3. DBC right: u5 = u6 = u9 = 1
## meaning that u1 = u8 and u5 = u9 are overwritten by 2 and 3 and we end up with
## u1 = 0
## u2 = u7
## u4 = 0
## u5 = 1
## u6 = 1
## u8 = 0
## u9 = 1
ch2 = ConstraintHandler(dh)
add!(ch2, PeriodicDirichlet(:u, collect_periodic_faces(grid, "bottom", "top")))
add!(ch2, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0))
add!(ch2, Dirichlet(:u, getfaceset(grid, "right"), (x, t) -> 1))
close!(ch2)
update!(ch2, 0)

K1 = create_sparsity_pattern(dh, ch1)
f1 = zeros(ndofs(dh))
a1 = start_assemble(K1, f1)
K2 = create_sparsity_pattern(dh, ch2)
f2 = zeros(ndofs(dh))
a2 = start_assemble(K2, f2)

for cell in CellIterator(dh)
assemble!(a1, celldofs(cell), ke, fe)
assemble!(a2, celldofs(cell), ke, fe)
end

# Equivalent assembly
@test K1 == K2
@test f1 == f2

# Equivalence after apply!
apply!(K1, f1, ch1)
apply!(K2, f2, ch2)
@test K1 == K2
@test f1 == f2
@test apply!(K1 \ f1, ch1) ≈ apply!(K2 \ f2, ch2)
end # subtestset

@testset "time dependence" begin
## Pure Dirichlet
ch1 = ConstraintHandler(dh)
add!(ch1, Dirichlet(:u, getfaceset(grid, "top"), (x, t) -> 3.0t + 2.0))
add!(ch1, Dirichlet(:u, getfaceset(grid, "bottom"), (x, t) -> 1.5t + 1.0))
add!(ch1, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 1.0t))
add!(ch1, Dirichlet(:u, getfaceset(grid, "right"), (x, t) -> 2.0t))
close!(ch1)
## Dirichlet with corresponding AffineConstraint on dof 2 and 7
ch2 = ConstraintHandler(dh)
add!(ch2, AffineConstraint(7, [8 => 1.0, 9 => 1.0], 2.0))
add!(ch2, AffineConstraint(2, [1 => 0.5, 5 => 0.5], 1.0))
add!(ch2, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 1.0t))
add!(ch2, Dirichlet(:u, getfaceset(grid, "right"), (x, t) -> 2.0t))
close!(ch2)

K1 = create_sparsity_pattern(dh, ch1)
f1 = zeros(ndofs(dh))
K2 = create_sparsity_pattern(dh, ch2)
f2 = zeros(ndofs(dh))

for t in (1.0, 2.0)
update!(ch1, t)
update!(ch2, t)
a1 = start_assemble(K1, f1)
a2 = start_assemble(K2, f2)
for cell in CellIterator(dh)
assemble!(a1, celldofs(cell), ke, fe)
assemble!(a2, celldofs(cell), ke, fe)
end
@test K1 == K2
@test f1 == f2
apply!(K1, f1, ch1)
apply!(K2, f2, ch2)
@test K1 == K2
@test f1 == f2
@test K1 \ f1 ≈ K2 \ f2
end

end # subtestset


@testset "error paths" begin
ch = ConstraintHandler(dh)
add!(ch, AffineConstraint(1, [2 => 1.0], 0.0))
add!(ch, AffineConstraint(2, [3 => 1.0], 0.0))
@test_throws ErrorException("nested affine constraints currently not supported") close!(ch)
end # subtestset

end # testset

0 comments on commit d68f9a1

Please sign in to comment.