Add high-level abstract coupling restriction

This may be much more user friendly

Author: pjaap

CHANGELOG.md

@@ -1,5 +1,11 @@
 # CHANGES
 
+## v1.9.0
+
+### Added
+  - It is now possible to create a `CoupledDofsRestriction(unknown, source_region, target_region, axis)` only from abstract data without knowledge of the grid or the FES.
+    It is assumed that the coupling is given between the extreme coordinates of the grid along the given axis.
+
 ## v1.8.0
 
 ### Added

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.8.0"
+version = "1.9.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

examples/Example312_PeriodicBoundary3D.jl

@@ -129,8 +129,7 @@ end
 
 function main(;
         order = 1,
-        periodic = true,
-        use_LM_restrictions = true,
+        periodic_coupling = :none, # :restriction, :operator, :high_level_restriction
         Plotter = nothing,
         force = 1.0,
         h = 1.0e-4,
@@ -164,16 +163,22 @@ function main(;
 
     assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet]))
 
-    if periodic
+    if periodic_coupling == :high_level_restriction
+
+        # new high-level call
+        assign_restriction!(PD, CoupledDofsRestriction(u, reg_left, reg_right, 1; parallel = threads > 1, threads))
+
+    elseif periodic_coupling != :none
+
         function give_opposite!(y, x)
             y .= x
             y[1] = width - x[1]
             return nothing
         end
         @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
-        if use_LM_restrictions
+        if periodic_coupling == :restriction
             assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
-        else
+        else # :operator
             assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
         end
     end
@@ -191,15 +196,18 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
 function runtests()                                                                                  #hide
-    sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
+    sol1, _ = main(periodic_coupling = :operator, threads = 1)                                       #hide
     @test abs(maximum(view(sol1[1])) - 1.8004602502175202) < 2.0e-3                                  #hide
 
-    sol2, _ = main(use_LM_restrictions = false, threads = 4)                                         #hide
+    sol2, _ = main(periodic_coupling = :operator, threads = 4)                                       #hide
     @test sol1.entries ≈ sol2.entries                                                                #hide
 
-    sol3, _ = main(use_LM_restrictions = true, threads = 4)                                          #hide
+    sol3, _ = main(periodic_coupling = :restriction, threads = 4)                                    #hide
     @test sol1.entries ≈ sol3.entries                                                                #hide
 
+    sol4, _ = main(periodic_coupling = :high_level_restriction, threads = 4)                         #hide
+    @test sol1.entries ≈ sol4.entries                                                                #hide
+
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

src/common_restrictions/coupled_dofs_restriction.jl

@@ -45,8 +45,65 @@ function CoupledDofsRestriction(matrices::Vector{AM}) where {AM <: AbstractMatri
 end
 
 
+"""
+    CoupledDofsRestriction(source::Ti, target::Ti, axis::Ti) where Ti
+
+    Create an incomplete CoupledDofsRestriction by only providing abstract information:
+    - unknown: the unknown to be coupled
+    - source: source boundary region
+    - target: target boundary region
+    - axis: axis along which the coupling is given
+
+    This constructor assumes that the coupling is given between the extrema of the grid coordinates
+    along the given axis.
+
+    The concrete coupling matrix will be computed in the solver, when the grid geometry and the FES is known.
+"""
+function CoupledDofsRestriction(unknown::Unknown, source_region::Ti, target_region::Ti, axis::Ti; kwargs...) where {Ti}
+    return CoupledDofsRestriction(
+        [nothing],
+        Dict{Symbol, Any}(
+            :name => "CoupledDofsRestriction",
+            :reduce_col_space => false,
+            :unknown => unknown,
+            :source_region => source_region,
+            :target_region => target_region,
+            :axis => axis,
+            :kwargs => kwargs
+        )
+    )
+end
+
+
 function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
 
+    # remember original R
+    R_orig = R
+
+    # compute the coupling matrix first, if necessary
+    if R.coupling_matrices[begin] === nothing
+        FES = sol[R.parameters[:unknown]].FES
+        source_region = R.parameters[:source_region]
+        target_region = R.parameters[:target_region]
+        axis = R.parameters[:axis]
+        R_kwargs = R.parameters[:kwargs]
+
+        xrange = @views extrema(FES.dofgrid[Coordinates][axis, :])
+
+        temp = sum(xrange)
+        function give_opposite!(a, b)
+            a .= b
+            a[axis] = temp - b[axis] # then, x ↔ y are opposite
+            return nothing
+        end
+
+        coupling_matrix = get_periodic_coupling_matrix(FES, source_region, target_region, give_opposite!; R_kwargs...)
+
+        # replace R (in this scope)
+        R = CoupledDofsRestriction(coupling_matrix)
+    end
+
+
     # extract all col indices and remove duplicates
     # subtract diagonal and shrink matrix to non-empty cols
     Bs = [ (matrix - LinearAlgebra.I)[:, unique(findnz(matrix)[2])] for matrix in R.coupling_matrices ]
@@ -63,11 +120,11 @@ function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
         B = B[:, cols_of_interest]
     end
 
-    R.parameters[:matrix] = B
-    R.parameters[:rhs] = zeros(size(B, 2))
+    R_orig.parameters[:matrix] = B
+    R_orig.parameters[:rhs] = zeros(size(B, 2))
 
     # fixed dofs are all active rows of B
-    R.parameters[:fixed_dofs] = unique(findnz(B)[1])
+    R_orig.parameters[:fixed_dofs] = unique(findnz(B)[1])
 
     return nothing
 end