From 6973f6c1042b5a03d569c6caf1e35ad082bb1487 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 15 May 2025 18:10:56 +0200 Subject: [PATCH 1/6] add missing dispatches to _mass_qr --- src/L2_projection.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 15cfdd8a8d..6c00045cd5 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -143,6 +143,12 @@ end function _mass_qr(::Lagrange{shape, 2}) where {shape <: RefSimplex} return QuadratureRule{shape}(4) end +function _mass_qr(::DiscontinuousLagrange{shape, order}) where {shape, order} + return _mass_qr(Lagrange{shape, order}()) +end +function _mass_qr(::Serendipity{shape, order}) where {shape, order} + return _mass_qr(Lagrange{shape, order}()) +end _mass_qr(ip::VectorizedInterpolation) = _mass_qr(ip.ip) function _assemble_L2_matrix(dh::DofHandler, qrs_lhs::Vector{<:QuadratureRule}) From 3cca5d25ce2b3bfbc870ffb720cc61af476b4e5b Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Fri, 16 May 2025 18:38:56 +0200 Subject: [PATCH 2/6] add test --- src/L2_projection.jl | 21 +++++--------- test/runtests.jl | 56 ++++++++++++++++++------------------ test/test_l2_projection.jl | 59 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 94 insertions(+), 42 deletions(-) diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 6c00045cd5..0db066942d 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -101,7 +101,7 @@ Add an interpolation `ip` on the cells in `set` to the `L2Projector` `proj`. projection equation, when calling [`project`](@ref). It should match the quadrature points used when creating the quadrature-point variables to project. * The *optional* `qr_lhs` sets the quadrature rule used to integrate the left-hand-side of the projection equation, - and defaults to a quadrature rule that integrates the mass-matrix exactly for the given interpolation `ip`. + and defaults to a quadrature rule that integrates the mass-matrix accurately for the given interpolation `ip`. """ function add!( @@ -137,19 +137,12 @@ function close!(proj::L2Projector) end # Quadrature sufficient for integrating a mass matrix -function _mass_qr(::Lagrange{shape, order}) where {shape <: AbstractRefShape, order} - return QuadratureRule{shape}(order + 1) -end -function _mass_qr(::Lagrange{shape, 2}) where {shape <: RefSimplex} - return QuadratureRule{shape}(4) -end -function _mass_qr(::DiscontinuousLagrange{shape, order}) where {shape, order} - return _mass_qr(Lagrange{shape, order}()) -end -function _mass_qr(::Serendipity{shape, order}) where {shape, order} - return _mass_qr(Lagrange{shape, order}()) -end -_mass_qr(ip::VectorizedInterpolation) = _mass_qr(ip.ip) +_mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip)) +_mass_qr_order(::Lagrange{<:AbstractRefShape, order}) where {order} = order + 1 +_mass_qr_order(::DiscontinuousLagrange{<:AbstractRefShape, order}) where {order} = order + 1 +_mass_qr_order(::Serendipity{<:AbstractRefShape, order}) where {order} = order + 1 +_mass_qr_order(::Lagrange{shape, 2}) where {shape <: RefSimplex} = 4 +_mass_qr_order(ip::VectorizedInterpolation) = _mass_qr_order(ip.ip) function _assemble_L2_matrix(dh::DofHandler, qrs_lhs::Vector{<:QuadratureRule}) M = Symmetric(allocate_matrix(dh)) diff --git a/test/runtests.jl b/test/runtests.jl index e27fd73f26..e9c19a8d68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,34 +16,34 @@ using HCubature: hcubature, hquadrature include("test_utils.jl") # Unit tests -include("test_collectionsofviews.jl") -include("test_refshapes.jl") -include("test_interpolations.jl") -include("test_cellvalues.jl") -include("test_facevalues.jl") -include("test_interfacevalues.jl") -include("test_quadrules.jl") -include("test_assemble.jl") -include("test_dofs.jl") -include("test_sparsity_patterns.jl") -include("test_constraints.jl") -include("test_grid_dofhandler_vtk.jl") -include("test_vtk_export.jl") -include("test_abstractgrid.jl") -include("test_grid_generators.jl") -include("test_grid_addboundaryset.jl") -include("test_mixeddofhandler.jl") +# include("test_collectionsofviews.jl") +# include("test_refshapes.jl") +# include("test_interpolations.jl") +# include("test_cellvalues.jl") +# include("test_facevalues.jl") +# include("test_interfacevalues.jl") +# include("test_quadrules.jl") +# include("test_assemble.jl") +# include("test_dofs.jl") +# include("test_sparsity_patterns.jl") +# include("test_constraints.jl") +# include("test_grid_dofhandler_vtk.jl") +# include("test_vtk_export.jl") +# include("test_abstractgrid.jl") +# include("test_grid_generators.jl") +# include("test_grid_addboundaryset.jl") +# include("test_mixeddofhandler.jl") include("test_l2_projection.jl") -include("test_pointevaluation.jl") -# include("test_notebooks.jl") -include("test_apply_rhs.jl") -include("test_apply_analytical.jl") -include("PoolAllocator.jl") -include("test_deprecations.jl") -include("blockarrays.jl") -include("test_assembler_extensions.jl") -include("test_continuity.jl") -include("test_examples.jl") +# include("test_pointevaluation.jl") +# # include("test_notebooks.jl") +# include("test_apply_rhs.jl") +# include("test_apply_analytical.jl") +# include("PoolAllocator.jl") +# include("test_deprecations.jl") +# include("blockarrays.jl") +# include("test_assembler_extensions.jl") +# include("test_continuity.jl") +# include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined # # See which is not defined if fails @@ -52,4 +52,4 @@ include("test_examples.jl") # end # Integration tests -include("integration/test_simple_scalar_convergence.jl") +# include("integration/test_simple_scalar_convergence.jl") diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index 035a916ff8..b73d570a7d 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -493,6 +493,64 @@ function test_l2proj_errorpaths() return end +function test_mass_qr() + function _ips(shape) + if shape <: RefCube + return ( + Lagrange{shape, 1}(), Lagrange{shape, 2}(), + Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), + DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), Serendipity{ + shape, 2, + }(), + ) + else + return ( + Lagrange{shape, 1}(), Lagrange{shape, 2}(), + Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), + DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), + ) + end + end + + function _mass_qr_test(ip::Interpolation{RefShape}) where {RefShape} + return QuadratureRule{RefShape}(Ferrite._mass_qr_order(ip) + 1) + end + + function _calculate_local_me(grid, ip::Interpolation, qr::QuadratureRule) + cv = CellValues(qr, ip) + cc = CellCache(grid) + Ferrite.reinit!(cc, 1) + Ferrite.reinit!(cv, cc) + + nbasef = getnbasefunctions(cv) + Me = zeros(nbasef, nbasef) + for q_point in 1:getnquadpoints(cv) + dV = getdetJdV(cv, q_point) + for i in 1:nbasef + Ni = shape_value(cv, q_point, i) + for j in 1:nbasef + Nj = shape_value(cv, q_point, j) + Me[i, j] += (Ni ⋅ Nj) * dV + end + end + end + return Me + end + + for (shape, refshape) in ((Quadrilateral, RefQuadrilateral), (Triangle, RefTriangle)) + grid = generate_grid(shape, (2, 2), Vec{2}((0.0, 0.0)), Vec{2}((1.0, 1.0))) + for ip in _ips(refshape) + me_1 = _calculate_local_me(grid, ip, Ferrite._mass_qr(ip)) + me_2 = _calculate_local_me(grid, ip, _mass_qr_test(ip)) + @test me_1 ≈ me_2 + if !(me_1 ≈ me_2) + println("The mass matrix calculated with order+1 is not identical to the mass matrix calculated with order for $ip") + end + end + end + return +end + @testset "Test L2-Projection" begin test_projection(1, RefQuadrilateral) test_projection(1, RefTriangle) @@ -505,4 +563,5 @@ end test_export(subset = true) test_show_l2() test_l2proj_errorpaths() + test_mass_qr() end From 9ad9e3ac70758ca2a485b0de895b86074ed881dc Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Fri, 16 May 2025 18:39:59 +0200 Subject: [PATCH 3/6] revert comment --- test/runtests.jl | 56 ++++++++++++++++++++++++------------------------ 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e9c19a8d68..e27fd73f26 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,34 +16,34 @@ using HCubature: hcubature, hquadrature include("test_utils.jl") # Unit tests -# include("test_collectionsofviews.jl") -# include("test_refshapes.jl") -# include("test_interpolations.jl") -# include("test_cellvalues.jl") -# include("test_facevalues.jl") -# include("test_interfacevalues.jl") -# include("test_quadrules.jl") -# include("test_assemble.jl") -# include("test_dofs.jl") -# include("test_sparsity_patterns.jl") -# include("test_constraints.jl") -# include("test_grid_dofhandler_vtk.jl") -# include("test_vtk_export.jl") -# include("test_abstractgrid.jl") -# include("test_grid_generators.jl") -# include("test_grid_addboundaryset.jl") -# include("test_mixeddofhandler.jl") +include("test_collectionsofviews.jl") +include("test_refshapes.jl") +include("test_interpolations.jl") +include("test_cellvalues.jl") +include("test_facevalues.jl") +include("test_interfacevalues.jl") +include("test_quadrules.jl") +include("test_assemble.jl") +include("test_dofs.jl") +include("test_sparsity_patterns.jl") +include("test_constraints.jl") +include("test_grid_dofhandler_vtk.jl") +include("test_vtk_export.jl") +include("test_abstractgrid.jl") +include("test_grid_generators.jl") +include("test_grid_addboundaryset.jl") +include("test_mixeddofhandler.jl") include("test_l2_projection.jl") -# include("test_pointevaluation.jl") -# # include("test_notebooks.jl") -# include("test_apply_rhs.jl") -# include("test_apply_analytical.jl") -# include("PoolAllocator.jl") -# include("test_deprecations.jl") -# include("blockarrays.jl") -# include("test_assembler_extensions.jl") -# include("test_continuity.jl") -# include("test_examples.jl") +include("test_pointevaluation.jl") +# include("test_notebooks.jl") +include("test_apply_rhs.jl") +include("test_apply_analytical.jl") +include("PoolAllocator.jl") +include("test_deprecations.jl") +include("blockarrays.jl") +include("test_assembler_extensions.jl") +include("test_continuity.jl") +include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined # # See which is not defined if fails @@ -52,4 +52,4 @@ include("test_l2_projection.jl") # end # Integration tests -# include("integration/test_simple_scalar_convergence.jl") +include("integration/test_simple_scalar_convergence.jl") From 9f319e74fa1b1525ca4efe92195ec01d8091f3fa Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Fri, 16 May 2025 18:40:53 +0200 Subject: [PATCH 4/6] fix test --- test/test_l2_projection.jl | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index b73d570a7d..1f57322d1c 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -495,21 +495,15 @@ end function test_mass_qr() function _ips(shape) + ips = ( + Lagrange{shape, 1}(), Lagrange{shape, 2}(), + Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), + DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), DiscontinuousLagrange{shape, 3}(), + ) if shape <: RefCube - return ( - Lagrange{shape, 1}(), Lagrange{shape, 2}(), - Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), - DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), Serendipity{ - shape, 2, - }(), - ) - else - return ( - Lagrange{shape, 1}(), Lagrange{shape, 2}(), - Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), - DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), - ) + ips = (ips..., Serendipity{shape, 2}()) end + return ips end function _mass_qr_test(ip::Interpolation{RefShape}) where {RefShape} From e40d45f2c0c4d27311ac634aff9d2fb5c9be6238 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 22 May 2025 13:45:28 +0200 Subject: [PATCH 5/6] fix for triangle and dl --- src/L2_projection.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 0db066942d..76284380c7 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -139,9 +139,9 @@ end # Quadrature sufficient for integrating a mass matrix _mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip)) _mass_qr_order(::Lagrange{<:AbstractRefShape, order}) where {order} = order + 1 -_mass_qr_order(::DiscontinuousLagrange{<:AbstractRefShape, order}) where {order} = order + 1 -_mass_qr_order(::Serendipity{<:AbstractRefShape, order}) where {order} = order + 1 -_mass_qr_order(::Lagrange{shape, 2}) where {shape <: RefSimplex} = 4 +_mass_qr_order(::DiscontinuousLagrange{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}()) +_mass_qr_order(::Serendipity{<:AbstractRefShape, order}) where {order} = _mass_qr_order(Lagrange{RefShape, order}()) +_mass_qr_order(::Lagrange{RefSimplex{refdim}, order}) where {order, refdim} = max(1, refdim * order) _mass_qr_order(ip::VectorizedInterpolation) = _mass_qr_order(ip.ip) function _assemble_L2_matrix(dh::DofHandler, qrs_lhs::Vector{<:QuadratureRule}) From 48b30cfb0c1d16841cd057bc8d34a7a02e3b8cd6 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 22 May 2025 14:04:49 +0200 Subject: [PATCH 6/6] fix serendipity --- src/L2_projection.jl | 2 +- test/test_l2_projection.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 76284380c7..511894803a 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -140,7 +140,7 @@ end _mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip)) _mass_qr_order(::Lagrange{<:AbstractRefShape, order}) where {order} = order + 1 _mass_qr_order(::DiscontinuousLagrange{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}()) -_mass_qr_order(::Serendipity{<:AbstractRefShape, order}) where {order} = _mass_qr_order(Lagrange{RefShape, order}()) +_mass_qr_order(::Serendipity{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}()) _mass_qr_order(::Lagrange{RefSimplex{refdim}, order}) where {order, refdim} = max(1, refdim * order) _mass_qr_order(ip::VectorizedInterpolation) = _mass_qr_order(ip.ip) diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index 1f57322d1c..ed286039d2 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -500,7 +500,7 @@ function test_mass_qr() Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(), DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), DiscontinuousLagrange{shape, 3}(), ) - if shape <: RefCube + if shape <: Ferrite.RefHypercube ips = (ips..., Serendipity{shape, 2}()) end return ips