diff --git a/examples/Example212_PeriodicBoundary2D.jl b/examples/Example212_PeriodicBoundary2D.jl new file mode 100644 index 00000000..e853ad79 --- /dev/null +++ b/examples/Example212_PeriodicBoundary2D.jl @@ -0,0 +1,180 @@ +#= + +# 212 : Periodic Boundary 2D +([source code](@__SOURCE_URL__)) + +This is a simple demonstration and validation of the generic periodic boundary operator. + +We construct an unstructured periodic 2D grid and solve a simple linear elastic problem +with periodic boundary coupling along the x-axis. + +![](example212.png) +=# + +module Example212_PeriodicBoundary2D + +using ExtendableFEM +using ExtendableGrids +using SimplexGridFactory +using GridVisualize +using Triangulate +using UnicodePlots +using StaticArrays +using LinearAlgebra +using Test #hide +using SparseArrays + +## enumerate the boundary regions +const reg_left = 4 +const reg_right = 2 +const reg_dirichlet = 1 +const reg_default = 3 + +## 2D reduction of the material used in Example 312 +## in Voigt notation +function material_tensor() + c11 = 396.0 + c12 = 137.0 + c44 = 116.0 + + return @SArray [ + c11 c12 0 + c12 c11 0 + 0 0 c44 + ] +end + +## generate the kernels for the linear problem +## ๐‚: Hooke tensor, ๐‘“: body force +function make_kernels(๐‚, ๐‘“) + + ## linear stress-strain mapping + bilinear_kernel!(ฯƒ, ฮตv, qpinfo) = mul!(ฯƒ, ๐‚, ฮตv) + + ## plain body force + linear_kernel!(result, qpinfo) = (result .= ๐‘“) + + return bilinear_kernel!, linear_kernel! +end + + +""" + create a 2D grid with Dirichlet boundary region at the bottom center +""" +function create_grid(; h, height, width) + builder = SimplexGridBuilder(; Generator = Triangulate) + + ## bottom points + b1 = point!(builder, 0, 0) + b2 = point!(builder, 0.45 * width, 0) + b3 = point!(builder, 0.5 * width, 0) + b4 = point!(builder, 0.55 * width, 0) + b5 = point!(builder, width, 0) + + ## top points + t1 = point!(builder, 0, height) + t2 = point!(builder, width / 2, height) + t3 = point!(builder, width, height) + + ## default faces + facetregion!(builder, reg_default) + facet!(builder, b1, b2) + facet!(builder, b4, b5) + facet!(builder, t1, t2) + facet!(builder, t2, t3) + + ## left face + facetregion!(builder, reg_left) + facet!(builder, b1, t1) + + ## right face + facetregion!(builder, reg_right) + facet!(builder, b5, t3) + + ## Dirichlet face + facetregion!(builder, reg_dirichlet) + facet!(builder, b3, b4) + facet!(builder, b2, b3) + + ## divider + facetregion!(builder, reg_default) + facet!(builder, t2, b3) + + + cellregion!(builder, 1) + maxvolume!(builder, h) + regionpoint!(builder, width / 3, height / 2) + + cellregion!(builder, 2) + # much finer grid on the right half to make periodic coupling non-trivial + maxvolume!(builder, 0.1 * h) + regionpoint!(builder, 2width / 3, height / 2) + + return simplexgrid(builder) +end + +function main(; + order = 1, + periodic = true, + Plotter = nothing, + force = 10.0, + h = 1.0e-2, + width = 6.0, + height = 1.0, + kwargs... + ) + xgrid = create_grid(; h, width, height) + + ## create finite element space and solution vector + if order == 1 + FES = FESpace{H1P1{2}}(xgrid) + elseif order == 2 + FES = FESpace{H1P2{2, 2}}(xgrid) + end + + ## problem description + PD = ProblemDescription() + u = Unknown("u"; name = "displacement") + assign_unknown!(PD, u) + + ๐‚ = material_tensor() + ๐‘“ = force * [0, 1] + + bilinear_kernel!, linear_kernel! = make_kernels(๐‚, ๐‘“) + assign_operator!(PD, BilinearOperator(bilinear_kernel!, [ฮตV(u, 1.0)]; kwargs...)) + assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...)) + + assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet])) + + if periodic + function give_opposite!(y, x) + y .= x + y[1] = width - x[1] + return nothing + end + + coupling_matrix = get_periodic_coupling_matrix(FES, xgrid, reg_left, reg_right, give_opposite!) + display(coupling_matrix) + assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...)) + end + + sol = solve(PD, FES) + + plt = GridVisualizer(; Plotter, size = (1300, 800)) + + magnification = 1 + displaced_grid = deepcopy(xgrid) + displace_mesh!(displaced_grid, sol[1], magnify = magnification) + gridplot!(plt, displaced_grid, linewidth = 1, title = "displaced mesh, $(magnification)x magnified", scene3d = :LScene) + + return sol, plt +end + +generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide +function runtests() #hide + sol, plt = main() #hide + @test abs(maximum(view(sol[1])) - 1.3447465095618172) < 1.0e-3 #hide + return nothing #hide +end #hide + +end # module diff --git a/examples/Example312_PeriodicBoundary3D.jl b/examples/Example312_PeriodicBoundary3D.jl new file mode 100644 index 00000000..002836e8 --- /dev/null +++ b/examples/Example312_PeriodicBoundary3D.jl @@ -0,0 +1,196 @@ +#= + +# 312 : Periodic Boundary 3D +([source code](@__SOURCE_URL__)) + +This is a simple demonstration and validation of the generic periodic boundary operator. + +We construct an unstructured periodic 3D grid and solve a simple linear elastic problem +with periodic coupling along the x-axis. + +![](example312.png) +=# + +module Example312_PeriodicBoundary3D + +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_dirichlet = 3 +const reg_default = 4 + + +# define the Hooke tensor for AlN (from DOI 10.1063/1.1368156) +function material_tensor() + c11 = 396.0 + c12 = 137.0 + c13 = 108.0 + c33 = 373.0 + c44 = 116.0 + + return @SArray [ + c11 c12 c13 0 0 0 + c12 c11 c13 0 0 0 + c13 c13 c33 0 0 0 + 0 0 0 c44 0 0 + 0 0 0 0 c44 0 + 0 0 0 0 0 c44 + ] +end + +## generate the kernels for the linear problem +## ๐‚: Hooke tensor, ๐‘“: body force +function make_kernels(๐‚, ๐‘“) + + ## linear stress-strain mapping + bilinear_kernel!(ฯƒ, ฮตv, qpinfo) = mul!(ฯƒ, ๐‚, ฮตv) + + ## body force + linear_kernel!(result, qpinfo) = (result .= ๐‘“) + + return bilinear_kernel!, linear_kernel! +end + + +""" + create 3D grid with Dirichlet boundary region at the bottom center +""" +function create_grid(; h, height, width, depth) + builder = SimplexGridBuilder(; Generator = TetGen) + + ## bottom points + b01 = point!(builder, 0, 0, 0) + b02 = point!(builder, 0.45 * width, 0, 0) + b03 = point!(builder, 0.55 * width, 0, 0) + b04 = point!(builder, width, 0, 0) + + b11 = point!(builder, 0, depth, 0) + b12 = point!(builder, 0.45 * width, depth, 0) + b13 = point!(builder, 0.55 * width, depth, 0) + b14 = point!(builder, width, depth, 0) + + ## top points + t01 = point!(builder, 0, 0, height) + t02 = point!(builder, width, 0, height) + + t11 = point!(builder, 0, depth, height) + t12 = point!(builder, width, depth, height) + + ## center points + c01 = point!(builder, 0.5 * width, 0, 0) + c02 = point!(builder, 0.5 * width, 0, height) + c11 = point!(builder, 0.5 * width, depth, 0) + c12 = point!(builder, 0.5 * width, depth, height) + + ## default faces + facetregion!(builder, reg_default) + facet!(builder, b01, b02, b12, b11) + facet!(builder, b03, b04, b14, b13) + facet!(builder, [t01, c02, c12, t11]) + facet!(builder, [c02, t02, t12, c12]) + facet!(builder, [b01, b02, c01, c02, t01]) + facet!(builder, [c01, b03, b04, t02, c02]) + facet!(builder, [b11, b12, c11, c12, t11]) + facet!(builder, [c11, b13, b14, t12, c12]) + facet!(builder, c01, c02, c12, c11) + + ## left face + facetregion!(builder, reg_left) + facet!(builder, b01, t01, t11, b11) + + ## right face + facetregion!(builder, reg_right) + facet!(builder, b04, t02, t12, b14) + + ## Dirichlet face + facetregion!(builder, reg_dirichlet) + facet!(builder, [b02, c01, c11, b12]) + facet!(builder, [c01, b03, b13, c11]) + + cellregion!(builder, 1) + maxvolume!(builder, h) + regionpoint!(builder, width / 3, depth / 2, height / 2) + cellregion!(builder, 2) + ## finer grid on the right half to make the periodic coupling non-trivial + maxvolume!(builder, 0.3 * h) + regionpoint!(builder, 2 * width / 3, depth / 2, height / 2) + + return simplexgrid(builder) +end + +function main(; + order = 1, + periodic = true, + Plotter = nothing, + force = 1.0, + h = 1.0e-4, + width = 6.0, + height = 0.2, + depth = 1, + kwargs... + ) + + xgrid = create_grid(; h, width, height, depth) + + ## create finite element space and solution vector + if order == 1 + FES = FESpace{H1P1{3}}(xgrid) + elseif order == 2 + FES = FESpace{H1P2{3, 3}}(xgrid) + end + + ## problem description + PD = ProblemDescription() + u = Unknown("u"; name = "displacement") + assign_unknown!(PD, u) + + ๐‚ = material_tensor() + ๐‘“ = force * [0, 0, 1] + + bilinear_kernel!, linear_kernel! = make_kernels(๐‚, ๐‘“) + assign_operator!(PD, BilinearOperator(bilinear_kernel!, [ฮตV(u, 1.0)]; kwargs...)) + assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...)) + + assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet])) + + if periodic + function give_opposite!(y, x) + y .= x + y[1] = width - x[1] + return nothing + end + coupling_matrix = get_periodic_coupling_matrix(FES, xgrid, reg_left, reg_right, give_opposite!) + display(coupling_matrix) + assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...)) + end + + ## solve + sol = solve(PD, FES; kwargs...) + + plt = GridVisualizer(; Plotter, size = (1200, 800)) + magnification = 1 + displaced_grid = deepcopy(xgrid) + displace_mesh!(displaced_grid, sol[1], magnify = magnification) + gridplot!(plt, displaced_grid, linewidth = 1, title = "displaced mesh, $(magnification)x magnified", scene3d = :LScene) + + return sol, plt + +end + +generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide +function runtests() #hide + sol, plt = main() #hide + @test abs(maximum(view(sol[1])) - 1.8038107003663026) < 1.0e-3 #hide + return nothing #hide +end #hide + +end # module diff --git a/src/common_operators/combinedofs.jl b/src/common_operators/combinedofs.jl index 034a4fac..1e9c3035 100644 --- a/src/common_operators/combinedofs.jl +++ b/src/common_operators/combinedofs.jl @@ -84,46 +84,83 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1} if CD.parameters[:verbosity] > 0 @info ".... coupling $(length(coupling_matrix.nzval)) dofs" end - function assemble(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T} + function assemble!(A::AbstractSparseArray{T}, b::AbstractVector{T}, assemble_matrix::Bool, assemble_rhs::Bool, kwargs...) where {T} if assemble_matrix + + # transpose the matrix once for efficient row access + transposed_coupling_matrix = sparse(transpose(coupling_matrix)) + + # go through each coupled dof and update the FE adjacency info + # from the constrained dofs here + for dof_i in 1:size(coupling_matrix, 2) # this col-view is efficient coupling_i = @views coupling_matrix[:, dof_i] + # do nothing if dof_k is not coupled to any constrained dof + if nnz(coupling_i) == 0 + continue + end + + # write the FE adjacency of the constrained dofs into this row + targetrow = dof_i + offsetX + + # extract the constrained dofs and the weights + coupled_dofs_i, weights_i = findnz(coupling_i) + + # parse through all cols and update the entries + for dof_j in 1:size(coupling_matrix, 2) + # this col-view is efficient + coupling_j = @views coupling_matrix[:, dof_j] + + # if both dof_i and dof_j are coupled to a constrained dof, then + # the FE adjacency A_ij is not updated: this is covered by the linear combinations + # expressed in the rows of the constrained dofs_on_boundary + # Hence, check that dof_j is not coupled to anything + if nnz(coupling_j) == 0 + targetcol = dof_j + offsetY + for (dof_k, weight_ik) in zip(coupled_dofs_i, weights_i) + sourcerow = dof_k + offsetX + sourcecol = targetcol + val = A[sourcerow, sourcecol] + _addnz(A, targetrow, targetcol, val, weight_ik) + end + end + end + end + + # replace the geometric coupling rows based + # on the original coupling matrix + for dof_i in 1:size(transposed_coupling_matrix, 2) + + coupling_i = transposed_coupling_matrix[:, dof_i] # do nothing if no coupling for dof_i if nnz(coupling_i) == 0 continue end # get the coupled dofs of dof_i and the corresponding weights - coupled_dofs, weights = findnz(coupling_i) + coupled_dofs_i, weights_i = findnz(coupling_i) - # transfer all assembly information of dof_i to the - # coupled dofs: add the corresponding matrix row - sourcerow = dof_i + offsetY - for sourcecol in 1:size(A, 2) - val = A[sourcerow, sourcecol] - for dof_j in coupled_dofs - targetrow = dof_j + offsetX - targetcol = sourcecol - _addnz(A, targetrow, targetcol, val, 1) - end - # eliminate the sourcerow - A[sourcerow, sourcecol] = 0 - end + sourcerow = dof_i + offsetX + # eliminate the sourcerow + for col in 1:size(A, 2) + A[sourcerow, col] = 0 + end # replace sourcerow with coupling linear combination _addnz(A, sourcerow, sourcerow, -1.0, 1) - for (dof_j, weight) in zip(coupled_dofs, weights) - # set negative weight for dofแตข - โˆ‘โฑผ wโฑผdofโฑผ = 0 - _addnz(A, sourcerow, dof_j + offsetY, weight, 1) + for (dof_j, weight_ij) in zip(coupled_dofs_i, weights_i) + # weights for โˆ‘โฑผ wโฑผdofโฑผ - dofแตข = 0 + _addnz(A, sourcerow, dof_j + offsetY, weight_ij, 1) end - flush!(A) end + flush!(A) end if assemble_rhs + for dof_i in 1:size(coupling_matrix, 2) # this col-view is efficient coupling_i = @views coupling_matrix[:, dof_i] @@ -132,22 +169,35 @@ function build_assembler!(CD::CombineDofs{UT, CT}, FE::Array{<:FEVectorBlock, 1} continue end + # get the coupled dofs of dof_i and the corresponding weights coupled_dofs, weights = findnz(coupling_i) - sourcerow = dof_i + offsetY - for dof_j in coupled_dofs - targetrow = dof_j + offsetX - b[targetrow] += b[sourcerow] + # transfer all assembly information to dof_i + targetrow = dof_i + offsetY + for (dof_j, weight) in zip(coupled_dofs, weights) + sourcerow = dof_j + offsetY + b[targetrow] += weight * b[sourcerow] + end + end + + + # now set the rows of the constrained dofs to zero to enforce the linear combination + for dof_i in 1:size(transposed_coupling_matrix, 2) + coupling_i = transposed_coupling_matrix[:, dof_i] + # do nothing if no coupling for dof_i + if nnz(coupling_i) == 0 + continue end - b[sourcerow] = 0.0 + b[dof_i + offsetX] = 0.0 + end end return nothing end - CD.assembler = assemble + CD.assembler = assemble! CD.FESX = FESX CD.FESY = FESY end diff --git a/test/Project.toml b/test/Project.toml index d4613150..c2d0d078 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -13,8 +13,10 @@ OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" diff --git a/test/runtests.jl b/test/runtests.jl index 6f0710c8..98e42924 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,6 +35,7 @@ function run_examples() "Example207_AdvectionUpwindDG.jl", "Example210_LshapeAdaptivePoissonProblem.jl", "Example211_LshapeAdaptiveEQPoissonProblem.jl", + "Example212_PeriodicBoundary2D.jl", "Example220_ReactionConvectionDiffusion.jl", "Example225_ObstacleProblem.jl", "Example226_Thermoforming.jl", @@ -54,6 +55,7 @@ function run_examples() "Example290_PoroElasticity.jl", "Example301_PoissonProblem.jl", "Example310_DivFreeBasis.jl", + "Example312_PeriodicBoundary3D.jl", ] return @testset "module examples" begin @@ -87,7 +89,10 @@ function run_all_tests() run_dgblf_tests() run_nonlinear_operator_tests() run_itemintegrator_tests() - return run_dt_tests() + run_dt_tests() + run_test_helper_functions() + + return nothing end run_all_tests() diff --git a/test/test_helper_functions.jl b/test/test_helper_functions.jl index c314ecb6..2b3bcaf7 100644 --- a/test/test_helper_functions.jl +++ b/test/test_helper_functions.jl @@ -1,102 +1,107 @@ -# check if coords are opposite to each other -function test_coords(X, Y, xgrid, is_opposite) - coords = xgrid[Coordinates] - for (x, y) in zip(X, Y) - if @views !is_opposite(coords[:, x], coords[:, y]) - @show coords[:, x] coords[:, y] - return false +function run_test_helper_functions() + + # check if coords are opposite to each other + function test_coords(X, Y, xgrid, is_opposite) + coords = xgrid[Coordinates] + for (x, y) in zip(X, Y) + if @views !is_opposite(coords[:, x], coords[:, y]) + @show coords[:, x] coords[:, y] + return false + end end + return true end - return true -end -@testset "get_periodic_coupling_info" begin + @testset "get_periodic_coupling_info" begin - xgrid = simplexgrid(0:100, 0:100) - # couple left <-> right - let - is_opposite(x, y) = abs(x[2] - y[2]) < 1.0e-12 - FES = FESpace{H1Pk{1, 2, 1}}(xgrid) - X, Y, F = get_periodic_coupling_info(FES, xgrid, 4, 2, is_opposite) - @test test_coords(X, Y, xgrid, is_opposite) - end - # couple top <-> bottom - let - is_opposite(x, y) = abs(x[1] - y[1]) < 1.0e-12 - FES = FESpace{H1Pk{1, 2, 1}}(xgrid) - X, Y, F = get_periodic_coupling_info(FES, xgrid, 1, 3, is_opposite) - @test test_coords(X, Y, xgrid, is_opposite) + xgrid = simplexgrid(0:100, 0:100) + # couple left <-> right + let + is_opposite(x, y) = abs(x[2] - y[2]) < 1.0e-12 + FES = FESpace{H1Pk{1, 2, 1}}(xgrid) + X, Y, F = get_periodic_coupling_info(FES, xgrid, 4, 2, is_opposite) + @test test_coords(X, Y, xgrid, is_opposite) + end + # couple top <-> bottom + let + is_opposite(x, y) = abs(x[1] - y[1]) < 1.0e-12 + FES = FESpace{H1Pk{1, 2, 1}}(xgrid) + X, Y, F = get_periodic_coupling_info(FES, xgrid, 1, 3, is_opposite) + @test test_coords(X, Y, xgrid, is_opposite) + end end -end -# check if matrix is a proper coupling matrix -# for testing, we assume that the FES on the "from" boundary is -# contained in the FES in the "to" boundary -function test_matrix(matrix; structured_grid = true) + # check if matrix is a proper coupling matrix + # for testing, we assume that the FES on the "from" boundary is + # contained in the FES in the "to" boundary + function test_matrix(matrix; structured_grid = true) + + if structured_grid + # only 1.0 entries! + if findfirst(โ‰‰(1.0, atol = 1.0e-8), matrix.nzval) !== nothing + @show "found entries โ‰‰ 1.0 " + return false + end - if structured_grid - # only 1.0 entries! - if findfirst(โ‰‰(1.0, atol = 1.0e-8), matrix.nzval) !== nothing - @show "found entries โ‰‰ 1.0 " - return false + # at most one entry in each col + for i in 1:size(matrix, 2) + if sum(matrix[:, i]) > 1.0 + 1.0e-8 + @show sum(matrix[:, i]) i + return false + end + end end - # at most one entry in each col - for i in 1:size(matrix, 2) - if sum(matrix[:, i]) > 1.0 + 1.0e-8 - @show sum(matrix[:, i]) i + # row sum is 0.0 or 1.0 + for i in 1:size(matrix, 1) + row_sum = sum(matrix[i, :]) + if !(row_sum == 0.0 || row_sum โ‰ˆ 1.0) + @show row_sum i return false end end - end - # row sum is 0.0 or 1.0 - for i in 1:size(matrix, 1) - row_sum = sum(matrix[i, :]) - if !(row_sum == 0.0 || row_sum โ‰ˆ 1.0) - @show row_sum i - return false - end + return true end - return true -end + @testset "get_periodic_coupling_matrix" begin -@testset "get_periodic_coupling_matrix" begin + # combine left/right at x=0/1 + function give_opposite!(y, x) + y .= x + y[1] = 1 - x[1] + return nothing + end - # combine left/right at x=0/1 - function give_opposite!(y, x) - y .= x - y[1] = 1 - x[1] - return nothing - end + let # 3D P1 + xgrid = simplexgrid(0:0.1:1.0, 0:0.1:1.0, 0:0.1:1.0) + FES = FESpace{H1P1{1}}(xgrid) + A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) + @test test_matrix(A) + end - let # 3D P1 - xgrid = simplexgrid(0:0.1:1.0, 0:0.1:1.0, 0:0.1:1.0) - FES = FESpace{H1P1{1}}(xgrid) - A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) - @test test_matrix(A) - end + let # 3D P2 with 2 components + xgrid = simplexgrid(0:0.5:1.0, 0:0.5:1.0, 0:0.5:1.0) + FES = FESpace{H1P2{2, 3}}(xgrid) + A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) + @test test_matrix(A) + end - let # 3D P2 with 2 components - xgrid = simplexgrid(0:0.5:1.0, 0:0.5:1.0, 0:0.5:1.0) - FES = FESpace{H1P2{2, 3}}(xgrid) - A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) - @test test_matrix(A) - end + let # 2D unstructured + b = SimplexGridBuilder(Generator = Triangulate) + grid1 = simplexgrid(0:1.0, 0:1.0) + grid2 = simplexgrid(0:1.0, 0:0.5:1.0) + bregions!(b, grid1, 1 => 1, 3 => 3, 4 => 4) + bregions!(b, grid2, 2 => 2) + xgrid = simplexgrid(b) - let # 2D unstructured - b = SimplexGridBuilder(Generator = Triangulate) - grid1 = simplexgrid(0:1.0, 0:1.0) - grid2 = simplexgrid(0:1.0, 0:0.5:1.0) - bregions!(b, grid1, 1 => 1, 3 => 3, 4 => 4) - bregions!(b, grid2, 2 => 2) - xgrid = simplexgrid(b) - - FES = FESpace{H1P1{1}}(xgrid) - A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) - @test test_matrix(A; structured_grid = false) + FES = FESpace{H1P1{1}}(xgrid) + A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8) + @test test_matrix(A; structured_grid = false) + end end + + return nothing end