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)` only from abstract data without knowledge of the grid or the FES.
+    It is assumed that the coupling is given between the source and target region and that these regions are parallel hyperplanes.
+
 ## 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; 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/ExtendableFEM.jl

@@ -40,7 +40,7 @@ using ExtendableFEMBase: ExtendableFEMBase, AbstractFiniteElement,
     norms, unicode_gridplot, unicode_scalarplot,
     update_basis!, SymmetricGradient
 using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
-    Adjacency, AssemblyType, BEdgeNodes, BFaceFaces,
+    Adjacency, AssemblyType, BEdgeNodes, BFaceFaces, BFaceNormals,
     BFaceNodes, BFaceRegions, CellAssemblyGroups,
     CellFaceOrientations, CellFaces, CellGeometries,
     CellNodes, CellRegions, Coordinates, EdgeNodes,

src/common_restrictions/coupled_dofs_restriction.jl

@@ -45,8 +45,77 @@ 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_region: source boundary region
+    - target_region: target boundary region
+
+    This constructor assumes that
+    - the source and target boundary regions form two parallel hyperplanes
+    - the coupling is performed orthogonally between the given hyperplanes
+
+    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; 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,
+            :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]
+        R_kwargs = R.parameters[:kwargs]
+
+        grid = FES.dofgrid
+
+        # at first, we select one face on source/target region
+        source_bface = findfirst(==(source_region), grid[BFaceRegions])
+        target_bface = findfirst(==(target_region), grid[BFaceRegions])
+
+        # the normal (it should really be opposite to the normal on the other side)
+        normal = grid[BFaceNormals][:, source_bface]
+        @assert normal ≈ -grid[BFaceNormals][:, target_bface]
+
+        # pick a coordinate on each boundary region
+        source_coord = grid[Coordinates][:, grid[BFaceNodes][1, source_bface]]
+        target_coord = grid[Coordinates][:, grid[BFaceNodes][1, target_bface]]
+
+        # compute the sum of scalar product with the normals (reflection point)
+        γ = source_coord'normal + target_coord'normal
+
+        function give_opposite!(y, x)
+            σ = 2.0 * normal'x
+            @. y = x + (γ - σ) * normal # then x ⇔ y are opposite along the normal vector
+            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 +132,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