Handle prescribed dofs in RHS of affine constraints

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 has to be done at runtime and not during close!(::ConstraintHandler). 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.
Ferrite-FEM · Nov 7, 2022 · b659a9e · b659a9e
1 parent f0efd13
commit b659a9e
Show file tree

Hide file tree

Showing 4 changed files with 146 additions and 6 deletions.
diff --git a/CHANGELOG.md b/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

diff --git a/docs/src/literate/stokes-flow.jl b/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

diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
@@ -133,7 +133,7 @@ function apply_rhs!(data::RHSData, f::AbstractVector, ch::ConstraintHandler, app
 
     # TODO: Can the loops be combined or does the order matter?
     @inbounds for i in 1:length(ch.inhomogeneities)
-        v = ch.inhomogeneities[i]
+        v = computed_inhomogeneity(ch.dofmapping, ch.dofcoefficients, ch.inhomogeneities, i)
         if !applyzero && v != 0
             for j in nzrange(K, i)
                 f[K.rowval[j]] -= v * K.nzval[j]
@@ -142,7 +142,7 @@ function apply_rhs!(data::RHSData, f::AbstractVector, ch::ConstraintHandler, app
     end
     @inbounds for (i, pdof) in pairs(ch.prescribed_dofs)
         dofcoef = ch.dofcoefficients[i]
-        b = ch.inhomogeneities[i]
+        b = computed_inhomogeneity(ch.dofmapping, ch.dofcoefficients, ch.inhomogeneities, i)
         if dofcoef !== nothing # if affine constraint
             for (d, v) in dofcoef
                 f[d] += f[pdof] * v
@@ -198,6 +198,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). See comments in
+    #   `computed_inhomogeneity`.
+    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
@@ -260,6 +279,9 @@ function add_prescribed_dof!(ch::ConstraintHandler, constrained_dof::Int, inhomo
     i = get(ch.dofmapping, constrained_dof, 0)
     if i != 0
         @debug @warn "dof $i already prescribed, overriding the old constraint"
+        # TODO: If this is a affine-prescribed dof we could potentially rewrite the
+        #       constraint to fulfill both, i.e. if we have `u1 = u2` and add `u1 = f(t)` we
+        #       can rewrite the first to `u2 = u1 = f(t)` automatically.
         ch.prescribed_dofs[i] = constrained_dof
         ch.inhomogeneities[i] = inhomogeneity
         ch.dofcoefficients[i] = dofcoefficients
@@ -586,7 +608,7 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con
     if !applyzero
         @inbounds for i in 1:length(ch.inhomogeneities)
             d = ch.prescribed_dofs[i]
-            v = ch.inhomogeneities[i]
+            v = computed_inhomogeneity(ch.dofmapping, ch.dofcoefficients, ch.inhomogeneities, i)
             if v != 0
                 for j in nzrange(K, d)
                     f[K.rowval[j]] -= v * K.nzval[j]
@@ -607,7 +629,8 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con
         d = ch.prescribed_dofs[i]
         K[d, d] = m
         if length(f) != 0
-            vz = applyzero ? zero(eltype(f)) : ch.inhomogeneities[i]
+            vz = applyzero ? zero(eltype(f)) :
+                 computed_inhomogeneity(ch.dofmapping, ch.dofcoefficients, ch.inhomogeneities, i)
             f[d] = vz * m
         end
     end
@@ -621,6 +644,27 @@ end
     return dofcoeffs[idx]
 end
 
+# Compute inhomogenity for a dof. For a Dirichlet constraint u1 = f1(t) this will just
+# return f1(t), but for affine constraints, e.g. u2 = w3 * u3 + w4 * u4 + b2 we allow master
+# dofs to also be prescribed by Dirichlet, e.g. u3 = f3(t). Since f3(t) is not known in
+# `close!` this needs to be resolved at runtime. For the example above, the computed
+# inhomogenity h2 for u2 becomes h2 = w3 * u3 + b2 = w3 * f3(t) + b2.
+@inline function computed_inhomogeneity(dofmapping, dofcoeffs, inhomogeneities, i)
+    # Note that `i` is not the prescribed dof, it is the index into ch.inhomogeneities/
+    # ch.dofcoefficients/ch.prescribed_dofs.
+    h = inhomogeneities[i]
+    coeffs = dofcoeffs[i]
+    coeffs === nothing && return h
+    for (d, w) in coeffs
+        j = get(dofmapping, d, 0)
+        j == 0 && continue
+        # If this dof is prescribed it must not have any master dofs (verified in close!)
+        @assert (jcoeffs = dofcoeffs[j]; jcoeffs === nothing || isempty(jcoeffs))
+        h += inhomogeneities[j] * w
+    end
+    return h
+end
+
 # Similar to Ferrite._condense!(K, ch), but only add the non-zero entries to K (that arises from the condensation process)
 function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, dofmapping::Dict{Int,Int}) where T
     ndofs = size(K, 1)

diff --git a/test/test_constraints.jl b/test/test_constraints.jl
@@ -1032,5 +1032,94 @@ 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
+
+    # 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)
+
+    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)
+    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 K1 \ f1 ≈ K2 \ f2
+
+    # Error paths
+    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