diff --git a/CHANGELOG.md b/CHANGELOG.md index c989cbbd..720021c6 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index 40d4a9eb..b153b361 100644 --- a/Project.toml +++ b/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 ", "Patrick Jaap "] [deps] diff --git a/docs/make.jl b/docs/make.jl index aac33327..87b51b5c 100644 --- a/docs/make.jl +++ b/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 diff --git a/docs/src/combinedofs.md b/docs/src/combinedofs.md index 5f977709..34b258b1 100644 --- a/docs/src/combinedofs.md +++ b/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 diff --git a/docs/src/examples_intro.md b/docs/src/examples_intro.md index de023c71..78dbc240 100644 --- a/docs/src/examples_intro.md +++ b/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) ``` diff --git a/examples/Example212_PeriodicBoundary2D.jl b/examples/Example212_PeriodicElasticity2D.jl similarity index 98% rename from examples/Example212_PeriodicBoundary2D.jl rename to examples/Example212_PeriodicElasticity2D.jl index 76e9307f..ecc97206 100644 --- a/examples/Example212_PeriodicBoundary2D.jl +++ b/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 diff --git a/examples/Example312_PeriodicBoundary3D.jl b/examples/Example312_PeriodicElasticity3D.jl similarity index 86% rename from examples/Example312_PeriodicBoundary3D.jl rename to examples/Example312_PeriodicElasticity3D.jl index e9809fee..c3f57f14 100644 --- a/examples/Example312_PeriodicBoundary3D.jl +++ b/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 diff --git a/examples/Example313_PeriodicPoisson.jl b/examples/Example313_PeriodicPoisson.jl new file mode 100644 index 00000000..3d663f2a --- /dev/null +++ b/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 diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 584931cd..30d428d4 100644 --- a/src/ExtendableFEM.jl +++ b/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, diff --git a/src/common_restrictions/boundarydata_restriction.jl b/src/common_restrictions/boundarydata_restriction.jl index 321a4e34..8777f52d 100644 --- a/src/common_restrictions/boundarydata_restriction.jl +++ b/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) diff --git a/src/common_restrictions/coupled_dofs_restriction.jl b/src/common_restrictions/coupled_dofs_restriction.jl index 79565eed..30ca61f5 100644 --- a/src/common_restrictions/coupled_dofs_restriction.jl +++ b/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 diff --git a/src/common_restrictions/linear_functional_restriction.jl b/src/common_restrictions/linear_functional_restriction.jl index 3ea37016..72895bce 100644 --- a/src/common_restrictions/linear_functional_restriction.jl +++ b/src/common_restrictions/linear_functional_restriction.jl @@ -72,7 +72,7 @@ function assemble!(R::LinearFunctionalRestriction{T, Tv}, sol, SC; kwargs...) wh assemble!(nothing, b, sol, R.linear_operator, SC; assemble_rhs = true, kwargs...) - R.parameters[:matrix] = reshape(b.entries, n, 1) + R.parameters[:matrix] = sparse(reshape(b.entries, n, 1)) R.parameters[:rhs] = Tv[R.value] R.parameters[:fixed_dofs] = [] diff --git a/src/solver_config.jl b/src/solver_config.jl index 1031ddb6..79bee97b 100644 --- a/src/solver_config.jl +++ b/src/solver_config.jl @@ -75,31 +75,32 @@ lastrhs(S::SolverConfiguration) = S.b # Default context information with help info. # default_solver_kwargs() = Dict{Symbol, Tuple{Any, String}}( - :target_residual => (1.0e-10, "stop if the absolute (nonlinear) residual is smaller than this number"), - :damping => (0, "amount of damping, value should be between in (0,1)"), :abstol => (1.0e-11, "abstol for linear solver (if iterative)"), - :reltol => (1.0e-11, "reltol for linear solver (if iterative)"), - :time => (0.0, "current time to be used in all time-dependent operators"), - :init => (nothing, "initial solution (also used to save the new solution)"), - :spy => (false, "show unicode spy plot of system matrix during solve"), - :symmetrize => (false, "make system matrix symmetric (replace by (A+A')/2)"), - :symmetrize_structure => (false, "make the system sparse matrix structurally symmetric (e.g. if [j,k] is also [k,j] must be set, all diagonal entries must be set)"), - :restrict_dofs => ([], "array of dofs for each unknown that should be solved (default: all dofs)"), :check_matrix => (false, "check matrix for symmetry and positive definiteness and largest/smallest eigenvalues"), - :verbosity => (0, "verbosity level"), - :show_config => (false, "show configuration at the beginning of solve"), - :show_matrix => (false, "show system matrix after assembly"), - :return_config => (false, "solver returns solver configuration (including A and b of last iteration)"), - :is_linear => ("auto", "linear problem (avoid reassembly of nonlinear operators to check residual)"), - :inactive => (Array{Unknown, 1}([]), "inactive unknowns (are made available in assembly, but not updated in solve)"), - :maxiterations => (10, "maximal number of nonlinear iterations/linear solves"), + :compress_restrictions => (true, "compress the column space of the given restriction matrices via QR to avoid overdetermined systems"), :constant_matrix => (false, "matrix is constant (skips reassembly and refactorization in solver)"), :constant_rhs => (false, "right-hand side is constant (skips reassembly)"), + :damping => (0, "amount of damping, value should be between in (0,1)"), + :inactive => (Array{Unknown, 1}([]), "inactive unknowns (are made available in assembly, but not updated in solve)"), + :init => (nothing, "initial solution (also used to save the new solution)"), + :initialized => (false, "linear system in solver configuration is already assembled (turns true after first solve)"), + :is_linear => ("auto", "linear problem (avoid reassembly of nonlinear operators to check residual)"), + :maxiterations => (10, "maximal number of nonlinear iterations/linear solves"), :method_linear => (UMFPACKFactorization(), "any solver or custom LinearSolveFunction compatible with LinearSolve.jl (default = UMFPACKFactorization())"), + :plot => (false, "plot all solved unknowns with a (very rough but fast) unicode plot"), :precon_linear => (nothing, "function that computes preconditioner for method_linear in case an iterative solver is chosen"), + :reltol => (1.0e-11, "reltol for linear solver (if iterative)"), + :restrict_dofs => ([], "array of dofs for each unknown that should be solved (default: all dofs)"), + :return_config => (false, "solver returns solver configuration (including A and b of last iteration)"), + :show_config => (false, "show configuration at the beginning of solve"), + :show_matrix => (false, "show system matrix after assembly"), + :spy => (false, "show unicode spy plot of system matrix during solve"), + :symmetrize => (false, "make system matrix symmetric (replace by (A+A')/2)"), + :symmetrize_structure => (false, "make the system sparse matrix structurally symmetric (e.g. if [j,k] is also [k,j] must be set, all diagonal entries must be set)"), + :target_residual => (1.0e-10, "stop if the absolute (nonlinear) residual is smaller than this number"), + :time => (0.0, "current time to be used in all time-dependent operators"), :timeroutputs => (:full, "configures show of timeroutputs (choose between :hide, :full, :compact)"), - :initialized => (false, "linear system in solver configuration is already assembled (turns true after first solve)"), - :plot => (false, "plot all solved unknowns with a (very rough but fast) unicode plot"), + :verbosity => (0, "verbosity level"), ) diff --git a/src/solvers.jl b/src/solvers.jl index b01cdde6..e06b5a71 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -271,6 +271,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, # block sizes for the block matrix block_sizes = [length(b_unrestricted), [ length(b) for b in restriction_rhss ]...] + total_size = sum(block_sizes) # total number of additional LM dofs nLMs = @views sum(block_sizes[2:end]) @@ -282,11 +283,46 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, rhs .-= B'sol_freedofs end - total_size = sum(block_sizes) Tv = eltype(b_unrestricted) ## create block matrix if linsolve_needs_matrix + if SC.parameters[:compress_restrictions] + # combine all restriction matrices into one: + combined_transposed_restriction_matrix = vcat(restriction_matrices'...) + combined_restriction_rhs = vcat(restriction_rhss...) + + # compute rank revealing QR decomposition (of the transposed matrix) + qr_result = qr(combined_transposed_restriction_matrix) + + # extract components (from docs: Q*R = combined_transposed_restriction_matrix[prow, pcol]) + (; Q, R, prow, pcol) = qr_result + + # we need the inverse column permutation + ipcol = invperm(pcol) + + # the new combined restriction block (add another transpose) + combined_restriction_matrix = R[:, ipcol]' + + # the new combined restriction rhs + combined_restriction_rhs = Q'combined_restriction_rhs[prow] + + # compress the column space + qr_rank = rank(qr_result) + @assert norm(combined_restriction_rhs[(qr_rank + 1):end]) ≤ 1.0e-12 * norm(combined_restriction_rhs) "the rhs of the restriction is not in the image" + + combined_restriction_matrix = combined_restriction_matrix[:, 1:qr_rank] + combined_restriction_rhs = combined_restriction_rhs[1:qr_rank] + + # replace by single entries + restriction_matrices = [combined_restriction_matrix] + restriction_rhss = [combined_restriction_rhs] + + # update sizes + block_sizes = [length(b_unrestricted), length(combined_restriction_rhs)] + total_size = sum(block_sizes) + end + A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes) A_block[Block(1, 1)] = A_unrestricted end @@ -298,7 +334,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, u_block = BlockVector(zeros(Tv, total_size), block_sizes) u_block[Block(1)] = u_unrestricted - for i in eachindex(PD.restrictions) + for i in eachindex(restriction_matrices) if linsolve_needs_matrix A_block[Block(1, i + 1)] = restriction_matrices[i] A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i]) diff --git a/test/runtests.jl b/test/runtests.jl index c03b0f00..2bcea4cb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,7 +36,7 @@ function run_examples() "Example207_AdvectionUpwindDG.jl", "Example210_LshapeAdaptivePoissonProblem.jl", "Example211_LshapeAdaptiveEQPoissonProblem.jl", - "Example212_PeriodicBoundary2D.jl", + "Example212_PeriodicElasticity2D.jl", "Example220_ReactionConvectionDiffusion.jl", "Example225_ObstacleProblem.jl", "Example226_Thermoforming.jl", @@ -56,7 +56,8 @@ function run_examples() "Example290_PoroElasticity.jl", "Example301_PoissonProblem.jl", "Example310_DivFreeBasis.jl", - "Example312_PeriodicBoundary3D.jl", + "Example312_PeriodicElasticity3D.jl", + "Example313_PeriodicPoisson.jl", ] return @testset "module examples" begin