Lagrange Restrictions

Project.toml

@@ -1,15 +1,17 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.4"
+version = "1.4.0"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 
 [deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 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"
@@ -24,6 +26,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 Aqua = "0.8"
+BlockArrays = "1.7.0"
 CommonSolve = "0.2"
 DiffResults = "1"
 DocStringExtensions = "0.8,0.9"
@@ -32,6 +35,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"

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(;
         end
 
         coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!)
-        display(coupling_matrix)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        # display(coupling_matrix)
+        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,8 +177,10 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide
 function runtests()                                                                                  #hide
-    sol, plt = main()                                                                                #hide
+    sol, plt = main(use_LM_restrictions = true)                                                      #hide
     @test abs(maximum(view(sol[1])) - 1.3447465095618172) < 1.0e-3                                   #hide
+    sol2, plt = main(use_LM_restrictions = false)                                                    #hide
+    @test sol.entries ≈ sol2.entries                                                                 #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(;
             return nothing
         end
         coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!)
-        display(coupling_matrix)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        # display(coupling_matrix)
+        if use_LM_restrictions
+            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
+        else
+            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        end
     end
 
     ## solve
@@ -185,8 +190,10 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
 function runtests()                                                                                  #hide
-    sol, plt = main()                                                                                #hide
-    @test abs(maximum(view(sol[1])) - 1.8004602502175202) < 2.0e-3                                  #hide
+    sol, plt = main(use_LM_restrictions = true)                                                      #hide
+    @test abs(maximum(view(sol[1])) - 1.8004602502175202) < 2.0e-3                                   #hide
+    sol2, plt = main(use_LM_restrictions = false)                                                    #hide
+    @test sol.entries ≈ sol2.entries                                                                 #hide
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

src/ExtendableFEM.jl

@@ -5,6 +5,7 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 """
 module ExtendableFEM
 
+using BlockArrays: BlockMatrix, BlockVector, Block, BlockedMatrix, blocks, axes
 using CommonSolve: CommonSolve
 using DiffResults: DiffResults
 using DocStringExtensions: DocStringExtensions, TYPEDEF, TYPEDSIGNATURES
@@ -59,6 +60,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!
@@ -67,7 +69,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 SparseDiffTools: SparseDiffTools, ForwardColorJacCache,
     forwarddiff_color_jacobian!, matrix_colors
 using Symbolics: Symbolics
@@ -121,11 +123,16 @@ include("common_operators/reduction_operator.jl")
 #export AbstractReductionOperator
 #export FixbyInterpolation
 
+include("restrictions.jl")
+include("common_restrictions/coupled_dofs_restriction.jl")
+export CoupledDofsRestriction
+
 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

src/common_restrictions/coupled_dofs_restriction.jl

@@ -0,0 +1,38 @@
+struct CoupledDofsRestriction{TM} <: AbstractRestriction
+    coupling_matrix::TM
+    parameters::Dict{Symbol, Any}
+end
+
+
+"""
+    CoupledDofsRestriction(matrix::AbstractMatrix)
+
+    Creates an restriction that couples dofs together.
+
+    The coupling is stored in a CSC Matrix `matrix`, s.t.,
+
+    dofᵢ = Σⱼ Aⱼᵢ dofⱼ (col-wise)
+
+    The matrix can be obtained from, e.g., `get_periodic_coupling_matrix`.
+"""
+function CoupledDofsRestriction(matrix::AbstractMatrix)
+    return CoupledDofsRestriction(matrix, Dict{Symbol, Any}(:name => "CoupledDofsRestriction"))
+end
+
+
+function assemble!(R::CoupledDofsRestriction, SC; kwargs...)
+
+    # extract all col indices
+    _, J, _ = findnz(R.coupling_matrix)
+
+    # remove duplicates
+    unique_cols = unique(J)
+
+    # subtract diagonal and shrink matrix to non-empty cols
+    B = (R.coupling_matrix - LinearAlgebra.I)[:, unique_cols]
+
+    R.parameters[:matrix] = B
+    R.parameters[:rhs] = Zeros(length(unique_cols))
+
+    return nothing
+end

src/problemdescription.jl

@@ -26,6 +26,11 @@ mutable struct ProblemDescription
     """
     operators::Array{AbstractOperator, 1}
     #reduction_operators::Array{AbstractReductionOperator,1}
+
+    """
+    A vector of Lagrange restrictions that are involved in the problem.
+    """
+    restrictions::Vector{AbstractRestriction}
 end
 
 """
@@ -41,7 +46,7 @@ Create an empty `ProblemDescription` with the given name.
 
 """
 function ProblemDescription(name = "My problem")
-    return ProblemDescription(name, Array{Unknown, 1}(undef, 0), Array{AbstractOperator, 1}(undef, 0))
+    return ProblemDescription(name, [], [], [])
 end
 
 
@@ -91,6 +96,23 @@ function assign_operator!(PD::ProblemDescription, o::AbstractOperator)
     return length(PD.operators)
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Adds a restrction to a `ProblemDescription`.
+
+# Arguments
+- `PD::ProblemDescription`: The problem description to which the operator should be added.
+- `r::AbstractRestriction`: The restriction to add.
+
+# Returns
+- The index (position) of the restriction in the `restrictions` array of `PD`.
+"""
+function assign_restriction!(PD::ProblemDescription, r::AbstractRestriction)
+    push!(PD.restrictions, r)
+    return length(PD.restrictions)
+end
+
 
 """
 $(TYPEDSIGNATURES)

src/restrictions.jl

@@ -0,0 +1,22 @@
+"""
+	AbstractRestriction
+
+Root type for all restrictions
+"""
+abstract type AbstractRestriction end
+
+
+function Base.show(io::IO, R::AbstractRestriction)
+    print(io, "AbstractRestriction")
+    return nothing
+end
+
+# informs solver when operator needs reassembly in a time dependent setting
+function is_timedependent(R::AbstractRestriction)
+    return false
+end
+
+function assemble!(R::AbstractRestriction, SC; kwargs...)
+    ## assembles internal restriction matrix in R
+    return nothing
+end