diff --git a/src/exports.jl b/src/exports.jl index 09a0ab1e5e..6ab0507a7a 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -17,6 +17,7 @@ export DiscontinuousLagrange, Serendipity, getnbasefunctions, + getrefshape, # Quadrature QuadratureRule, diff --git a/src/interpolations.jl b/src/interpolations.jl index efe0f8b195..8aa03dd533 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1268,6 +1268,7 @@ end ####################################### # Taken from https://defelement.com/elements/bubble-enriched-lagrange.html getnbasefunctions(::BubbleEnrichedLagrange{RefTriangle,1}) = 4 +adjust_dofs_during_distribution(::BubbleEnrichedLagrange{RefTriangle,1}) = false vertexdof_indices(::BubbleEnrichedLagrange{RefTriangle,1}) = ((1,), (2,), (3,)) facedof_indices(::BubbleEnrichedLagrange{RefTriangle,1}) = ((1,2), (2,3), (3,1)) diff --git a/test/integration/test_simple_scalar_convergence.jl b/test/integration/test_simple_scalar_convergence.jl new file mode 100644 index 0000000000..d6387fc0d3 --- /dev/null +++ b/test/integration/test_simple_scalar_convergence.jl @@ -0,0 +1,216 @@ +using Ferrite, Test +import Ferrite: getdim, default_interpolation + +module ConvergenceTestHelper + +using Ferrite, SparseArrays, ForwardDiff, Test +import LinearAlgebra: diag + +get_geometry(::Ferrite.Interpolation{RefLine}) = Line +get_geometry(::Ferrite.Interpolation{RefQuadrilateral}) = Quadrilateral +get_geometry(::Ferrite.Interpolation{RefTriangle}) = Triangle +get_geometry(::Ferrite.Interpolation{RefPrism}) = Wedge +get_geometry(::Ferrite.Interpolation{RefHexahedron}) = Hexahedron +get_geometry(::Ferrite.Interpolation{RefTetrahedron}) = Tetrahedron +get_geometry(::Ferrite.Interpolation{RefPyramid}) = Pyramid + +get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = 2*order+1 +get_quadrature_order(::BubbleEnrichedLagrange{shape, order}) where {shape, order} = 2*order + +get_num_elements(::Ferrite.Interpolation{shape, 1}) where {shape} = 21 +get_num_elements(::Ferrite.Interpolation{shape, 2}) where {shape} = 7 +get_num_elements(::Ferrite.Interpolation{RefHexahedron, 1}) = 11 +get_num_elements(::Ferrite.Interpolation{RefHexahedron, 2}) = 4 +get_num_elements(::Ferrite.Interpolation{shape, 3}) where {shape} = 8 +get_num_elements(::Ferrite.Interpolation{shape, 4}) where {shape} = 5 +get_num_elements(::Ferrite.Interpolation{shape, 5}) where {shape} = 3 + +analytical_solution(x) = prod(cos, x*π/2) +analytical_rhs(x) = -Tensors.laplace(analytical_solution,x) + +# Standard assembly copy pasta for Poisson problem +function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues, coords) + n_basefuncs = getnbasefunctions(cellvalues) + ## Reset to 0 + fill!(Ke, 0) + fill!(fe, 0) + ## Loop over quadrature points + for q_point in 1:getnquadpoints(cellvalues) + ## Get the quadrature weight + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + ## Loop over test shape functions + for i in 1:n_basefuncs + δu = shape_value(cellvalues, q_point, i) + ∇δu = shape_gradient(cellvalues, q_point, i) + ## Add contribution to fe + fe[i] += analytical_rhs(x) * δu * dΩ + ## Loop over trial shape functions + for j in 1:n_basefuncs + ∇u = shape_gradient(cellvalues, q_point, j) + ## Add contribution to Ke + Ke[i, j] += (∇δu ⋅ ∇u) * dΩ + end + end + end + return Ke, fe +end + +# Standard assembly copy pasta for Poisson problem +function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler) + ## Allocate the element stiffness matrix and element force vector + n_basefuncs = getnbasefunctions(cellvalues) + Ke = zeros(n_basefuncs, n_basefuncs) + fe = zeros(n_basefuncs) + ## Allocate global force vector f + f = zeros(ndofs(dh)) + ## Create an assembler + assembler = start_assemble(K, f) + ## Loop over all cels + for cell in CellIterator(dh) + ## Reinitialize cellvalues for this cell + reinit!(cellvalues, cell) + coords = getcoordinates(cell) + ## Compute element contribution + assemble_element!(Ke, fe, cellvalues, coords) + ## Assemble Ke and fe into K and f + assemble!(assembler, celldofs(cell), Ke, fe) + end + return K, f +end + +# Compute norms +function check_and_compute_convergence_norms(dh, u, cellvalues, testatol) + L2norm = 0.0 + ∇L2norm = 0.0 + L∞norm = 0.0 + for cell in CellIterator(dh) + reinit!(cellvalues, cell) + n_basefuncs = getnbasefunctions(cellvalues) + coords = getcoordinates(cell) + uₑ = u[celldofs(cell)] + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + uₐₙₐ = prod(cos, x*π/2) + uₐₚₚᵣₒₓ = function_value(cellvalues, q_point, uₑ) + L∞norm = max(L∞norm, norm(uₐₙₐ-uₐₚₚᵣₒₓ)) + L2norm += norm(uₐₙₐ-uₐₚₚᵣₒₓ)^2*dΩ + + ∇uₐₙₐ = gradient(x-> prod(cos, x*π/2), x) + ∇uₐₚₚᵣₒₓ = function_gradient(cellvalues, q_point, uₑ) + ∇L2norm += norm(∇uₐₙₐ-∇uₐₚₚᵣₒₓ)^2*dΩ + + # Pointwise convergence + @test uₐₙₐ ≈ uₐₚₚᵣₒₓ atol=testatol + end + end + √(L2norm), √(∇L2norm), L∞norm +end + +# Assemble and solve +function solve(dh, ch, cellvalues) + K, f = assemble_global(cellvalues, create_sparsity_pattern(dh), dh); + apply!(K, f, ch) + u = K \ f; +end + +function setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + # Construct Ferrite stuff + dh = DofHandler(grid) + add!(dh, :u, interpolation) + close!(dh); + + ch = ConstraintHandler(dh); + ∂Ω = union( + values(grid.facesets)... + ); + dbc = Dirichlet(:u, ∂Ω, (x, t) -> analytical_solution(x)) + add!(ch, dbc); + close!(ch); + + cellvalues = CellValues(qr, interpolation, interpolation_geo); + + return dh, ch, cellvalues +end + +end # module ConvergenceTestHelper + +# These test only for convergence within margins +@testset "convergence analysis" begin + @testset "$interpolation" for interpolation in ( + Lagrange{RefTriangle, 3}(), + Lagrange{RefTriangle, 4}(), + Lagrange{RefTriangle, 5}(), + Lagrange{RefHexahedron, 1}(), + Lagrange{RefTetrahedron, 1}(), + Lagrange{RefPrism, 1}(), + Lagrange{RefPyramid, 1}(), + # + Serendipity{RefQuadrilateral, 2}(), + Serendipity{RefHexahedron, 2}(), + # + BubbleEnrichedLagrange{RefTriangle, 1}(), + # + CrouzeixRaviart{RefTriangle, 1}(), + ) + # Generate a grid ... + geometry = ConvergenceTestHelper.get_geometry(interpolation) + interpolation_geo = default_interpolation(geometry) + N = ConvergenceTestHelper.get_num_elements(interpolation) + grid = generate_grid(geometry, ntuple(x->N, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 1e-2) + end +end + +# These test also for correct convergence rates +@testset "convergence rate" begin + @testset "$interpolation" for interpolation in ( + Lagrange{RefLine, 1}(), + Lagrange{RefLine, 2}(), + Lagrange{RefQuadrilateral, 1}(), + Lagrange{RefQuadrilateral, 2}(), + Lagrange{RefQuadrilateral, 3}(), + Lagrange{RefTriangle, 1}(), + Lagrange{RefTriangle, 2}(), + Lagrange{RefHexahedron, 2}(), + Lagrange{RefTetrahedron, 2}(), + Lagrange{RefPrism, 2}(), + ) + # Generate a grid ... + geometry = ConvergenceTestHelper.get_geometry(interpolation) + interpolation_geo = default_interpolation(geometry) + # "Coarse case" + N₁ = ConvergenceTestHelper.get_num_elements(interpolation) + grid = generate_grid(geometry, ntuple(x->N₁, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + L2₁, H1₁, _ = ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 1e-2) + + # "Fine case" + N₂ = 2*N₁ + grid = generate_grid(geometry, ntuple(x->N₂, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + L2₂, H1₂, _ = ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 5e-3) + + @test -(log(L2₂)-log(L2₁))/(log(N₂)-log(N₁)) ≈ Ferrite.getorder(interpolation)+1 atol=0.1 + @test -(log(H1₂)-log(H1₁))/(log(N₂)-log(N₁)) ≈ Ferrite.getorder(interpolation) atol=0.1 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d11ed6d582..fe3c1692c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,8 @@ if HAS_EXTENSIONS && MODULE_CAN_BE_TYPE_PARAMETER end include("test_utils.jl") + +# Unit tests include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") @@ -38,3 +40,6 @@ include("test_deprecations.jl") HAS_EXTENSIONS && include("blockarrays.jl") include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined + +# Integration tests +include("integration/test_simple_scalar_convergence.jl")