Handle prescribed dofs in RHS of affine constraints

CHANGELOG.md

@@ -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
@@ -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

docs/src/literate/stokes-flow.jl

@@ -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
@@ -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

src/Dofs/ConstraintHandler.jl

@@ -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
@@ -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
@@ -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))
@@ -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
@@ -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
@@ -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
@@ -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

test/test_constraints.jl

@@ -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