Lagrange Restrictions (#79)

adds restrictions to ProblemDescription to enforce constraints like zero mean value, coupled dofs and boundary data by Lagrange multipliers

Co-authored-by: chmerdon <[email protected]>

Author: pjaap

CHANGELOG.md

@@ -5,6 +5,9 @@
 ### Added
 
   - `compute_periodic_coupling_matrix` can be computed thread parallel
+  - new "Restriction" objects using Lagrange multipliers internally
+  - user can apply linear `LinearFunctionalRestriction`, `ZeroMeanValueRestriction`, `MassRestriction`, `CoupledDofsRestriction` or `BoundaryDataRestriction` to the `ProblemDescription` by `apply_restriction!()`
+  - Incorporated new restriction objects in examples 201, 252, 312 and 212.
 
 ## v1.4.0
 

Project.toml

@@ -4,13 +4,15 @@ version = "1.5.0"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 
 [deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -26,6 +28,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 Aqua = "0.8"
+BlockArrays = "1.7.0"
 ChunkSplitters = "3.1.2"
 CommonSolve = "0.2"
 DiffResults = "1"
@@ -35,6 +38,7 @@ ExplicitImports = "1"
 ExtendableFEMBase = "1.3.0"
 ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
+FillArrays = "1.13.0"
 ForwardDiff = "0.10.35,1"
 GridVisualize = "1.8.1"
 IncompleteLU = "0.2.1"

docs/make.jl

@@ -82,6 +82,7 @@ function make_all(; with_examples::Bool = true, modules = :all, run_examples::Bo
                 "nonlinearoperator.md",
                 "bilinearoperator.md",
                 "linearoperator.md",
+                "restrictions.md",
                 "interpolateboundarydata.md",
                 "homogeneousdata.md",
                 "fixdofs.md",

docs/src/restrictions.md

@@ -0,0 +1,81 @@
+# Restrictions
+
+ExtendableFEM.jl provides functionality to apply various restrictions to your finite element problems through the `AbstractRestriction` type system. Restrictions enforce additional constraints on your solution that cannot be easily expressed through standard boundary conditions or the weak formulation.
+
+## Built-in Restrictions
+
+### Linear Functional Restriction
+
+`LinearFunctionalRestriction` allows you to constrain a linear functional of a finite element unknown to a specific value. This is useful, for example, for enforcing that the solution has mean zero, or that an integral constraint is satisfied.
+
+Documentation:
+```@autodocs
+Modules = [ExtendableFEM]
+Pages = ["common_restrictions/linear_functional_restriction.jl"]
+Order   = [:type, :function]
+```
+
+#### Example Usage
+```julia
+# Restrict the mean value of unknown u to 0
+restriction = ZeroMeanValueRestriction(u)
+
+# Restrict the mean value to 1.0 with a specific operator
+restriction = MassRestriction(u, value = 1.0, operator = MyCustomOperator)
+
+# Assign to the problem
+assign_restriction!(PD, restriction)
+```
+
+### Coupled DOFs Restriction
+
+`CoupledDofsRestriction` enables coupling between different degrees of freedom (DOFs) in your system. This is particularly useful for implementing periodic boundary conditions or other constraints where DOFs need to be related to each other. Compared to manipulating the system matrix directly via operators (e.g. `CombineDofs`), `CoupledDofsRestriction` is much faster for large systems.
+
+Documentation:
+```@autodocs
+Modules = [ExtendableFEM]
+Pages = ["common_restrictions/coupled_dofs_restriction.jl"]
+Order   = [:type, :function]
+```
+
+#### Example Usage
+```julia
+# Create a coupling matrix for periodic boundary conditions
+coupling_matrix = get_periodic_coupling_matrix(...)
+restriction = CoupledDofsRestriction(coupling_matrix)
+assign_restriction!(PD, restriction)
+```
+
+### Boundary Data Restriction
+
+`BoundaryDataRestriction` enforces prescribed boundary values for a finite element unknown. It can be used for both homogeneous (zero) and non-homogeneous boundary conditions, interpolating the boundary data as needed.
+
+Documentation:
+```@autodocs
+Modules = [ExtendableFEM]
+Pages = ["common_restrictions/boundarydata_restriction.jl"]
+Order   = [:type, :function]
+```
+
+#### Example Usage
+```julia
+# Homogeneous boundary data (default)
+restriction = BoundaryDataRestriction(u)
+
+# Inhomogeneous boundary data using a function
+restriction = BoundaryDataRestriction(u, data = (result, qpinfo) -> result .= sin(qpinfo.x[1]))
+
+assign_restriction!(PD, restriction)
+```
+
+## AbstractRestriction API
+
+All restrictions are subtypes of `AbstractRestriction`. The following functions are available for interacting with restrictions:
+
+- `assemble!(restriction, solution, SC; kwargs...)`: Assemble the internal matrices and vectors for the restriction.
+- `restriction_matrix(restriction)`: Get the restriction matrix.
+- `restriction_rhs(restriction)`: Get the right-hand side vector for the restriction.
+- `fixed_dofs(restriction)`: Get the list of fixed degrees of freedom affected by the restriction.
+- `is_timedependent(restriction)`: Returns whether the restriction is time-dependent.
+
+---

examples/Example201_PoissonProblem.jl

@@ -35,13 +35,27 @@ function default_f!(fval, qpinfo)
     return nothing
 end
 
-function main(; μ = 1.0, nrefs = 4, order = 2, Plotter = nothing, parallel = false, f! = default_f!, npart = parallel ? 8 : 1, kwargs...)
+function main(;
+        μ = 1.0,
+        nrefs = 4,
+        order = 2,
+        f! = default_f!,
+        use_restriction = true,
+        parallel = false,
+        npart = parallel ? 8 : 1,
+        Plotter = nothing,
+        kwargs...
+    )
     ## problem description
     PD = ProblemDescription()
     assign_unknown!(PD, u)
     assign_operator!(PD, BilinearOperator([grad(u)]; parallel = parallel, factor = μ, kwargs...))
     assign_operator!(PD, LinearOperator(f!, [id(u)]; parallel = parallel, kwargs...))
-    assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4))
+    if use_restriction
+        assign_restriction!(PD, BoundaryDataRestriction(u; regions = 1:4, value = 0))
+    else
+        assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4))
+    end
 
     ## discretize
     xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
@@ -65,6 +79,8 @@ function runtests() #hide
     @test sum(sol.entries) ≈ 1.1140313632246377 #hide
     sol_parallel, plt = main(; μ = 1.0, nrefs = 2, order = 2, parallel = true, npart = 2) #hide
     @assert sum((sol_parallel.entries .- sol.entries) .^ 2) ≈ 0.0 #hide
+    sol_restrict, plt = main(; μ = 1.0, nrefs = 2, use_restriction = false, parallel = false, order = 2, npart = 2) #hide
+    @assert sqrt(sum((sol_restrict.entries .- sol.entries) .^ 2)) < 1.0e-15 #hide
     return nothing #hide
 end #hide
 end # module

examples/Example212_PeriodicBoundary2D.jl

@@ -116,6 +116,7 @@ end
 function main(;
         order = 1,
         periodic = true,
+        use_LM_restrictions = true,
         Plotter = nothing,
         force = 10.0,
         h = 1.0e-2,
@@ -154,8 +155,12 @@ function main(;
             return nothing
         end
 
-        @time coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
+        if use_LM_restrictions
+            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
+        else
+            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        end
     end
 
     sol = solve(PD, FES)
@@ -172,10 +177,15 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide
 function runtests()                                                                                  #hide
-    sol, _ = main()                                                                                  #hide
-    @test abs(maximum(view(sol[1])) - 1.3447465095618172) < 1.0e-3                                   #hide
-    sol2, _ = main(threads = 4)                                                                      #hide
-    @test sol.entries ≈ sol2.entries                                                                 #hide
+    sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
+    @test abs(maximum(view(sol1[1])) - 1.3447465095618172) < 1.0e-3                                  #hide
+
+    sol2, _ = main(use_LM_restrictions = false, threads = 4)                                         #hide
+    @test sol1.entries ≈ sol2.entries                                                                #hide
+
+    sol3, _ = main(use_LM_restrictions = true, threads = 4)                                          #hide
+    @test sol1.entries ≈ sol3.entries                                                                #hide
+
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

examples/Example252_NSEPlanarLatticeFlow.jl

@@ -83,7 +83,14 @@ function exact_error!(result, u, qpinfo)
     return result .= result .^ 2
 end
 
-function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kwargs...)
+function main(;
+        μ = 0.001,
+        nrefs = 5,
+        reconstruct = true,
+        Plotter = nothing,
+        use_LM_restrictions = true,
+        kwargs...
+    )
 
     ## problem description
     PD = ProblemDescription()
@@ -96,6 +103,11 @@ function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kw
     assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id_u, grad(u), id(p)]; params = [μ], kwargs...))
     assign_operator!(PD, LinearOperator(f!(μ), [id_u]; kwargs...))
     assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:4))
+    if use_LM_restrictions
+        assign_restriction!(PD, ZeroMeanValueRestriction(p))
+    else
+        assign_operator!(PD, FixDofs(p; dofs = [1]))
+    end
 
     ## grid
     xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
@@ -106,10 +118,14 @@ function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kw
     ## solve
     sol = solve(PD, FES; kwargs...)
 
-    ## move integral mean of pressure
     pintegrate = ItemIntegrator([id(p)])
     pmean = sum(evaluate(pintegrate, sol)) / sum(xgrid[CellVolumes])
-    view(sol[p]) .-= pmean
+    if use_LM_restrictions
+        @show pmean
+    else
+        ## move integral mean of pressure
+        view(sol[p]) .-= pmean
+    end
 
     ## error calculation
     ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u), id(p)]; quadorder = 4, params = [μ], kwargs...)
@@ -127,7 +143,9 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example252_NSEPlanarLatticeFlow, "example252.png") #hide
 function runtests() #hide
-    L2errorU, plt = main(; nrefs = 4) #hide
+    L2errorU, _ = main(; use_LM_restrictions = false, nrefs = 4) #hide
+    @test L2errorU ≈ 0.010616923333947861 #hide
+    L2errorU, _ = main(; use_LM_restrictions = true, nrefs = 4) #hide
     @test L2errorU ≈ 0.010616923333947861 #hide
     return nothing #hide
 end #hide

examples/Example312_PeriodicBoundary3D.jl

@@ -130,6 +130,7 @@ end
 function main(;
         order = 1,
         periodic = true,
+        use_LM_restrictions = true,
         Plotter = nothing,
         force = 1.0,
         h = 1.0e-4,
@@ -169,8 +170,12 @@ function main(;
             y[1] = width - x[1]
             return nothing
         end
-        @time coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
+        if use_LM_restrictions
+            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
+        else
+            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        end
     end
 
     ## solve
@@ -185,10 +190,15 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
 function runtests()                                                                                  #hide
-    sol, _ = main()                                                                                  #hide
-    @test abs(maximum(view(sol[1])) - 1.8004602502175202) < 2.0e-3                                   #hide
-    sol2, _ = main(threads = 4)                                                                      #hide
-    @test sol.entries ≈ sol2.entries                                                                                 #hide
+    sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
+    @test abs(maximum(view(sol1[1])) - 1.8004602502175202) < 2.0e-3                                  #hide
+
+    sol2, _ = main(use_LM_restrictions = false, threads = 4)                                         #hide
+    @test sol1.entries ≈ sol2.entries                                                                #hide
+
+    sol3, _ = main(use_LM_restrictions = true, threads = 4)                                          #hide
+    @test sol1.entries ≈ sol3.entries                                                                #hide
+
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

src/ExtendableFEM.jl

@@ -6,6 +6,7 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 module ExtendableFEM
 
 using ChunkSplitters: chunks
+using BlockArrays: BlockMatrix, BlockVector, Block, blocks, axes, mortar
 using CommonSolve: CommonSolve
 using DiffResults: DiffResults
 using DocStringExtensions: DocStringExtensions, TYPEDEF, TYPEDSIGNATURES
@@ -60,6 +61,7 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     MTExtendableSparseMatrixCSC, findindex,
     rawupdateindex!
+using FillArrays: Zeros
 using ForwardDiff: ForwardDiff
 using GridVisualize: GridVisualize, GridVisualizer, gridplot!, reveal, save,
     scalarplot!, vectorplot!
@@ -68,7 +70,7 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
 using Printf: Printf, @printf, @sprintf
 using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz,
-    nzrange, rowvals, sparse, SparseVector
+    nzrange, rowvals, sparse, SparseVector, spzeros
 using StaticArrays: @MArray
 using SparseDiffTools: SparseDiffTools, ForwardColorJacCache,
     forwarddiff_color_jacobian!, matrix_colors
@@ -123,11 +125,14 @@ include("common_operators/reduction_operator.jl")
 #export AbstractReductionOperator
 #export FixbyInterpolation
 
+include("restrictions.jl")
+
 include("problemdescription.jl")
 export ProblemDescription
 export assign_unknown!
 export assign_operator!
 export replace_operator!
+export assign_restriction!
 
 include("helper_functions.jl")
 export get_periodic_coupling_info
@@ -196,6 +201,15 @@ export FixDofs
 include("common_operators/discface_interpolator.jl")
 export FaceInterpolator
 
+include("common_restrictions/boundarydata_restriction.jl")
+export BoundaryDataRestriction
+include("common_restrictions/coupled_dofs_restriction.jl")
+export CoupledDofsRestriction
+include("common_restrictions/linear_functional_restriction.jl")
+export LinearFunctionalRestriction
+export ZeroMeanValueRestriction
+export MassRestriction
+
 include("plots.jl")
 export plot_convergencehistory!
 export scalarplot!

src/common_operators/homogeneousdata_operator.jl

@@ -7,7 +7,7 @@ mutable struct HomogeneousData{UT, AT} <: AbstractOperator
 end
 
 fixed_dofs(O::HomogeneousData) = O.bdofs
-fixed_vals(O::HomogeneousData) = O.parameters[:value]
+fixed_vals(O::HomogeneousData) = O.parameters[:value] * ones(length(fixed_dofs(O)))
 
 
 default_homdata_kwargs() = Dict{Symbol, Tuple{Any, String}}(
@@ -143,6 +143,7 @@ function assemble!(O::HomogeneousData{UT, AT}, FES = O.FES; offset = 0, kwargs..
                 end
             end
         end
+        unique!(bdofs)
         if O.parameters[:verbosity] > 0
             @info "$(O.parameters[:name]) : penalizing $(length(bdofs)) dofs of '$(O.u.name)' ($AT, regions = $(O.parameters[:regions]))"
         end