Merge pull request #101 from WIAS-PDELib/feature/high-level-coupling

Add high-level abstract coupling restriction

Author: pjaap

CHANGELOG.md

@@ -1,5 +1,13 @@
 # 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.
+  - New `SolverConfig` key word argument `:compress_restrictions` (default = `true`) to call a QR decomposition on all restriction matrices to eliminate linearly dependent columns.
+  - New example `Example313_PeriodicPoisson` to cover corner cases of multiple restrictions at once.
+
 ## 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]

docs/make.jl

@@ -27,7 +27,7 @@ function make_all(; with_examples::Bool = true, modules = :all, run_examples::Bo
                 "Example207_AdvectionUpwindDG.jl",
                 "Example210_LshapeAdaptivePoissonProblem.jl",
                 "Example211_LshapeAdaptiveEQPoissonProblem.jl",
-                "Example212_PeriodicBoundary2D.jl",
+                "Example212_PeriodicElasticity2D.jl",
                 "Example220_ReactionConvectionDiffusion.jl",
                 "Example225_ObstacleProblem.jl",
                 "Example226_Thermoforming.jl",
@@ -51,7 +51,8 @@ function make_all(; with_examples::Bool = true, modules = :all, run_examples::Bo
                 "Example295_SlidingDroplet.jl",
                 "Example301_PoissonProblem.jl",
                 "Example310_DivFreeBasis.jl",
-                "Example312_PeriodicBoundary3D.jl",
+                "Example312_PeriodicElasticity3D.jl",
+                "Example313_PeriodicPoisson.jl",
                 "Example330_HyperElasticity.jl",
             ]
         end

docs/src/combinedofs.md

@@ -19,7 +19,7 @@ get_periodic_coupling_matrix
 
 ## Example: Periodic Boundary Coupling (extract from Example212)
 
-Suppose you want to enforce periodicity between the left and right boundaries of a 2D domain. You can use `get_periodic_coupling_matrix` to find the corresponding DOFs, and then use `CombineDofs` to couple them in the system. The following code is adapted from [Example212_PeriodicBoundary2D.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl/blob/main/examples/Example212_PeriodicBoundary2D.jl):
+Suppose you want to enforce periodicity between the left and right boundaries of a 2D domain. You can use `get_periodic_coupling_matrix` to find the corresponding DOFs, and then use `CombineDofs` to couple them in the system. The following code is adapted from [Example212_PeriodicElasticity2D.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl/blob/main/examples/Example212_PeriodicElasticity2D.jl):
 
 ```julia
 # Define a function to map points on the left to the right boundary

docs/src/examples_intro.md

@@ -11,19 +11,19 @@ The examples in this package are designed to be practical, reproducible, and edu
 
 ## Running the Examples
 
-To run an example (e.g., `Example212_PeriodicBoundary2D`):
+To run an example (e.g., `Example212_PeriodicElasticity2D`):
 
 1. Download the example file (see the source code link at the top of the example page).
 2. Ensure all required packages are installed in your Julia environment.
 3. In the Julia REPL:
 
     ```julia
-    julia> include("Example212_PeriodicBoundary2D.jl")
-    julia> Example212_PeriodicBoundary2D.main()
+    julia> include("Example212_PeriodicElasticity2D.jl")
+    julia> Example212_PeriodicElasticity2D.main()
     ```
 
 4. Some examples offer visual output via the optional argument `Plotter = PyPlot` or `Plotter = GLMakie` (provided the package is installed and loaded):
 
     ```julia
-    julia> Example212_PeriodicBoundary2D.main(Plotter = PyPlot)
+    julia> Example212_PeriodicElasticity2D.main(Plotter = PyPlot)
     ```

examples/Example212_PeriodicBoundary2D.jl …mples/Example212_PeriodicElasticity2D.jlexamples/Example212_PeriodicBoundary2D.jl renamed to examples/Example212_PeriodicElasticity2D.jl examples/Example212_PeriodicBoundary2D.jl renamed to examples/Example212_PeriodicElasticity2D.jl

@@ -11,7 +11,7 @@ with periodic boundary coupling along the x-axis.
 ![](example212.png)
 =#
 
-module Example212_PeriodicBoundary2D
+module Example212_PeriodicElasticity2D
 
 using ExtendableFEM
 using ExtendableGrids
@@ -176,7 +176,7 @@ function main(;
     return sol, plt
 end
 
-generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide
+generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicElasticity2D, "example212.png") #hide
 function runtests()                                                                                  #hide
     sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
     @test abs(maximum(view(sol1[1])) - 1.3447465095618172) < 1.0e-3                                  #hide

examples/Example312_PeriodicBoundary3D.jl …mples/Example312_PeriodicElasticity3D.jlexamples/Example312_PeriodicBoundary3D.jl renamed to examples/Example312_PeriodicElasticity3D.jl examples/Example312_PeriodicBoundary3D.jl renamed to examples/Example312_PeriodicElasticity3D.jl

@@ -11,7 +11,7 @@ with periodic coupling along the x-axis.
 ![](example312.png)
 =#
 
-module Example312_PeriodicBoundary3D
+module Example312_PeriodicElasticity3D
 
 using ExtendableFEM
 using ExtendableGrids
@@ -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
@@ -189,17 +194,20 @@ function main(;
 
 end
 
-generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
+generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicElasticity3D, "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
 

examples/Example313_PeriodicPoisson.jl

@@ -0,0 +1,144 @@
+#=
+
+# 312 : Periodic Poisson 3D
+([source code](@__SOURCE_URL__))
+
+This is a simple demonstration and validation of the new restriction based periodic boundary operator.
+
+An unstructured cube grid is coupled along two axes periodically and restricted to non-zero Dirichlet values at the other to faces.
+
+The result is a linear function and the error is measured directly.
+
+Note that the result is independent of the periodic coupling! Therefore the correctness of the
+new QR based column compression is covered by this example.
+
+![](example313.png)
+=#
+
+module Example313_PeriodicPoisson
+
+using ExtendableFEM
+using ExtendableGrids
+using SimplexGridFactory
+using GridVisualize
+using TetGen
+using UnicodePlots
+using StaticArrays
+using LinearAlgebra
+using Test #hide
+
+const reg_left = 1
+const reg_right = 2
+const reg_front = 3
+const reg_back = 4
+const reg_bottom = 5
+const reg_top = 6
+
+
+function create_grid(; h)
+    builder = SimplexGridBuilder(; Generator = TetGen)
+
+    ## bottom points
+    b00 = point!(builder, 0, 0, 0)
+    b01 = point!(builder, 0, 1, 0)
+    b10 = point!(builder, 1, 0, 0)
+    b11 = point!(builder, 1, 1, 0)
+
+    ## top points
+    t00 = point!(builder, 0, 0, 1)
+    t01 = point!(builder, 0, 1, 1)
+    t10 = point!(builder, 1, 0, 1)
+    t11 = point!(builder, 1, 1, 1)
+
+    ## left face
+    facetregion!(builder, reg_left)
+    facet!(builder, b00, b01, t01, t00)
+
+    ## right face
+    facetregion!(builder, reg_right)
+    facet!(builder, b10, b11, t11, t10)
+
+    ## front face
+    facetregion!(builder, reg_front)
+    facet!(builder, b00, b10, t10, t00)
+
+    ## back face
+    facetregion!(builder, reg_back)
+    facet!(builder, b01, b11, t11, t01)
+
+    ## top face
+    facetregion!(builder, reg_top)
+    facet!(builder, t00, t10, t11, t01)
+
+    ## bottom face
+    facetregion!(builder, reg_bottom)
+    facet!(builder, b00, b10, b11, b01)
+
+
+    cellregion!(builder, 1)
+    maxvolume!(builder, h)
+    regionpoint!(builder, 0.5, 0.5, 0.5)
+
+    return simplexgrid(builder)
+end
+
+function main(;
+        Plotter = nothing,
+        h = 1.0e-3,
+        periodic = true,
+        kwargs...
+    )
+
+    xgrid = create_grid(; h)
+
+    FES = FESpace{H1P1{1}}(xgrid)
+
+    ## problem description
+    PD = ProblemDescription("Periodic Poisson Problem")
+    u = Unknown("u"; name = "temperature")
+    assign_unknown!(PD, u)
+
+    assign_operator!(PD, BilinearOperator([grad(u)]; kwargs...))
+
+    assign_restriction!(PD, BoundaryDataRestriction(u; value = -1, regions = [reg_bottom]))
+    assign_restriction!(PD, BoundaryDataRestriction(u; value = +1, regions = [reg_top]))
+
+    if periodic
+        assign_restriction!(PD, CoupledDofsRestriction(u, reg_left, reg_right))
+        assign_restriction!(PD, CoupledDofsRestriction(u, reg_front, reg_back))
+    end
+
+    ## solve
+    sol, SC = solve(PD, FES; return_config = true, kwargs...)
+    residual(SC) < 1.0e-10 || error("Residual is not zero!")
+
+    function exact_error!(out, u, qpinfo)
+        # exact solution here is u(x,y,z) = 2z - 1
+        val = qpinfo.x[3] * 2 - 1
+        out[1] = (val - u[1])^2
+        return nothing
+    end
+
+    plt = plot([grid(u), id(u)], sol; Plotter, width = 1200, height = 800, scene3d = :LScene, slice = :y => 0.5)
+
+    ErrorIntegrator = ItemIntegrator(exact_error!, [id(u)]; quadorder = 2)
+    L2error = sqrt(sum(evaluate(ErrorIntegrator, sol)))
+
+    @show L2error
+
+    return L2error, plt
+
+end
+
+generateplots = ExtendableFEM.default_generateplots(Example313_PeriodicPoisson, "example313.png") #hide
+function runtests()                                                                               #hide
+    error1, _ = main(periodic = true)                                                             #hide
+    @test error1 < 1.0e-12                                                                        #hide
+
+    error2, _ = main(periodic = false)                                                            #hide
+    @test error2 < 1.0e-12                                                                        #hide
+
+    return nothing                                                                                #hide
+end                                                                                               #hide
+
+end # module

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/boundarydata_restriction.jl

@@ -34,6 +34,14 @@ function assemble!(R::BoundaryDataRestriction, sol, SC; kwargs...)
 
     ## assemble boundary operator (-> computes fixed dofs and vals)
     assemble!(nothing, SC.b, sol, R.boundary_operator, SC; kwargs...)
+    apply_penalties!(
+        nothing, SC.b, sol, R.boundary_operator, SC;
+        assemble_matrix = false,
+        assemble_rhs = false,
+        assemble_sol = true, # only modify the solution with boundary data
+        kwargs...
+    )
+
     fixeddofs = fixed_dofs(R.boundary_operator)
     fixedvals = fixed_vals(R.boundary_operator)
     nvals = length(fixeddofs)