diff --git a/CHANGELOG.md b/CHANGELOG.md index ffdd99e..31497f9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # CHANGES +## v1.5.0 November 03, 2025 + - added type parameter to `PointEvaluator` so `lazy_interpolate!` is now accessible for automatic differentiation + - added new `compute_lazy_interpolation_jacobian` function for computing Jacobians of interpolations between finite element spaces + - added `H1Pk_to_HDIVRT0_interpolator` function for (faster) computation of interpolation matrix from `H1Pk`-conforming spaces into `HDIVRT0`-conforming spaces on the same grid + - fixed interpolation of `HCURLN1` finite element + ## v1.4.0 July 17, 2025 - added new weighted reconstruction operator `WeightedReconstruct` (so far only for H1BR{2}) - fixed H1CR{2} reconstruction diff --git a/Project.toml b/Project.toml index 29998fc..3780ef9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,10 +1,11 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" +version = "1.5.0" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "1.4.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" @@ -13,6 +14,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" +SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" SpecialPolynomials = "a25cea48-d430-424a-8ee7-0d3ad3742e9e" [weakdeps] @@ -23,6 +26,7 @@ ExtendableFEMBaseUnicodePlotsExt = ["UnicodePlots"] [compat] DiffResults = "1" +DifferentiationInterface = "0.7.10" DocStringExtensions = "0.8,0.9" ExtendableGrids = "1.13.0" ExtendableSparse = "1.5.1" @@ -31,6 +35,8 @@ LinearAlgebra = "1.9" Polynomials = "2.0.21, 3, 4" Printf = "1.9" SparseArrays = "1.9" +SparseConnectivityTracer = "1.1.2" +SparseMatrixColorings = "0.4.22" SpecialPolynomials = "0.4.9, 0.5" UnicodePlots = "3.6" julia = "1.9" diff --git a/docs/Project.toml b/docs/Project.toml index c512582..7a359d3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,17 +1,17 @@ [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" -Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" PlutoStaticHTML = "359b1769-a58e-495b-9770-312e911026ad" PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" -DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" [compat] Literate = ">=0.2.7" diff --git a/docs/make.jl b/docs/make.jl index a41f763..15e7cc8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -50,6 +50,7 @@ function make_all(; with_examples::Bool = true) doctest = true, pages = [ "Home" => "index.md", + "Changelog" => "changelog.md", "Index" => "package_index.md", "List of Finite Elements" => "fems.md", "Base Structures" => Any[ diff --git a/docs/src/changelog.md b/docs/src/changelog.md new file mode 100644 index 0000000..d5dc3f5 --- /dev/null +++ b/docs/src/changelog.md @@ -0,0 +1,6 @@ +```@eval +using Markdown +Markdown.parse(""" +$(read("../../CHANGELOG.md",String)) +""") +``` diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md index 2d867ca..5b68183 100644 --- a/docs/src/interpolations.md +++ b/docs/src/interpolations.md @@ -37,6 +37,14 @@ Additionally, you can transfer finite element functions from one grid to another lazy_interpolate! ``` +```@docs +compute_lazy_interpolation_jacobian +``` + +```@docs +H1Pk_to_HDIVRT0_interpolator +``` + The following function continuously interpolates finite element function into a H1Pk space by point evaluations at the Lagrange nodes of the H1Pk element (averaged over all neighbours). diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index 6d3cef6..ed17935 100644 --- a/src/ExtendableFEMBase.jl +++ b/src/ExtendableFEMBase.jl @@ -46,6 +46,7 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry, using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!, AbstractExtendableSparseMatrixCSC, ExtendableSparseMatrixCSC, MTExtendableSparseMatrixCSC, rawupdateindex! +using DifferentiationInterface: AutoForwardDiff, AutoSparse, jacobian using ForwardDiff: ForwardDiff, DiffResults using LinearAlgebra: LinearAlgebra, convert, det, diagm, dot, eigen, ldiv!, lu, mul!, norm, transpose @@ -53,6 +54,8 @@ using Polynomials: Polynomials, Polynomial, coeffs using Printf: Printf, @printf using SparseArrays: SparseArrays, AbstractSparseArray, AbstractSparseMatrix, SparseMatrixCSC, nzrange, rowvals +using SparseConnectivityTracer: TracerSparsityDetector +using SparseMatrixColorings: GreedyColoringAlgorithm using SpecialPolynomials: SpecialPolynomials, ShiftedLegendre, basis include("functionoperators.jl") @@ -171,6 +174,9 @@ export PointEvaluator, evaluate!, evaluate_bary!, eval_func, eval_func_bary include("lazy_interpolate.jl") export lazy_interpolate! +include("interpolation_matrix_representations.jl") +export compute_lazy_interpolation_jacobian +export H1Pk_to_HDIVRT0_interpolator # ExtendableFEMBaseUnicodePlotsExt extension diff --git a/src/fedefs/h1_p1.jl b/src/fedefs/h1_p1.jl index e14625e..192c09a 100644 --- a/src/fedefs/h1_p1.jl +++ b/src/fedefs/h1_p1.jl @@ -41,7 +41,7 @@ end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT} # delegate edge nodes to node interpolation subitems = slice(FE.dofgrid[EdgeNodes], items) - return interpolaste!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...) + return interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...) end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT} diff --git a/src/fedefs/hcurl_n1.jl b/src/fedefs/hcurl_n1.jl index b7e3b3a..dbfebd0 100644 --- a/src/fedefs/hcurl_n1.jl +++ b/src/fedefs/hcurl_n1.jl @@ -27,6 +27,8 @@ get_dofmap_pattern(FEType::Type{<:HCURLN1{2}}, ::Union{Type{FaceDofs}, Type{BFac isdefined(FEType::Type{<:HCURLN1}, ::Type{<:Triangle2D}) = true +interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HCURLN1{2}}, ::Type{<:Triangle2D}) = 6 + function N1_tangentflux_eval_2d!(result, f, qpinfo) result[1] = -f[1] * qpinfo.normal[2] # rotated normal = tangent result[1] += f[2] * qpinfo.normal[1] @@ -34,6 +36,7 @@ function N1_tangentflux_eval_2d!(result, f, qpinfo) return nothing end init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HCURLN1{2}, APT} = FunctionalInterpolator(N1_tangentflux_eval_2d!, FES, ON_FACES; bonus_quadorder = 1) +init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: HCURLN1{2}, APT} = MomentInterpolator(FES, ON_CELLS) function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HCURLN1, APT} @@ -56,7 +59,7 @@ end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], kwargs...) where {Tv, Ti, FEType <: HCURLN1, APT} edim = get_ncomponents(FEType) - return if edim == 2 + if edim == 2 # delegate cell faces to face interpolation subitems = slice(FE.dofgrid[CellFaces], items) interpolate!(Target, FE, ON_FACES, data; items = subitems, kwargs...) @@ -65,6 +68,9 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, subitems = slice(FE.dofgrid[CellEdges], items) interpolate!(Target, FE, ON_EDGES, data; items = subitems, kwargs...) end + + # set values of interior N1 functions such that P0 moments are preserved + return get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...) end # on faces dofs are only tangential fluxes diff --git a/src/fedefs/l2_p1.jl b/src/fedefs/l2_p1.jl index 2fe07f2..5c67e10 100644 --- a/src/fedefs/l2_p1.jl +++ b/src/fedefs/l2_p1.jl @@ -43,7 +43,7 @@ end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: L2P1, APT} # delegate edge nodes to node interpolation subitems = slice(FE.dofgrid[EdgeNodes], items) - return interpolaste!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...) + return interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...) end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: L2P1, APT} diff --git a/src/interpolation_matrix_representations.jl b/src/interpolation_matrix_representations.jl new file mode 100644 index 0000000..0ffd28e --- /dev/null +++ b/src/interpolation_matrix_representations.jl @@ -0,0 +1,103 @@ +""" +```` +function compute_lazy_interpolation_jacobian( + target_space::FESpace, + source_space::FESpace; + use_cellparents::Bool = false, + kwargs... +) +```` + +Compute the Jacobian of the [`lazy_interpolate!`](@ref) call with respect to the +`source_space` degrees of freedom, i.e. for functions ``v = \\sum_j \\alpha_j \\, \\varphi_j`` of the +`source_space` and the interpolation operator ``I(v) = \\sum_k L_k(v)\\,\\phi_k = \\sum_k L_k\\left(\\sum_j \\alpha_j \\varphi_j\\right) \\, \\phi_k`` +into the `target_space`, this function computes the jacobian ``\\left[\\frac{\\partial L_k}{\\partial \\alpha_j}\\right]_{k,\\,j}`` +and returns its sparse matrix representation. + +# Arguments +- `target_space::FESpace`: Finite element space into which the interpolation ``I(v)`` is directed. +- `source_space::FESpace`: Finite element space from which ``v`` is taken. + +# Keyword Arguments +- `use_cellparents`: Use parent cell information if available (can speed up the calculation if the `target_space` is defined on a subgrid of `source_space`). +- `kwargs...`: Additional keyword arguments passed to lower-level `lazy_interpolate!` call. + +# Notes +- Since [`lazy_interpolate!`](@ref) is based on evaluating functions from the `source_space` +using a [`ExtendableFEMBase.PointEvaluator`](@ref), this should be used carefully on finer +grids as this is not the most efficient method, but will work out of the box for any +`source` and `target` spaces. +- This function can be used for computing prolongation or restriction operators if the `FESpace`s are defined on coarser/finer grids, respectively. + +""" +function compute_lazy_interpolation_jacobian(target_space::FESpace, source_space::FESpace; use_cellparents::Bool = false, kwargs...) + # DifferentiationInterface.jacobian needs a function of signature + # AbstractVector -> AbstractVector + function do_interpolation(source_vector::AbstractVector; use_cellparents = use_cellparents) + T = valtype(source_vector) + target_vector = FEVector{T}(target_space)[1] + source_FE_Vector = FEVector{T}(source_space) + source_FE_Vector.entries .= source_vector + + lazy_interpolate!(target_vector, source_FE_Vector, [(1, Identity)]; use_cellparents, kwargs...) + + return target_vector.entries + end + + n = ndofs(source_space) + + dense_forward_backend = AutoForwardDiff() + sparse_forward_backend = AutoSparse( + dense_forward_backend; + sparsity_detector = TracerSparsityDetector(), + coloring_algorithm = GreedyColoringAlgorithm(), + ) + + M = jacobian(x -> do_interpolation(x; use_cellparents), sparse_forward_backend, ones(n)) + + return M +end + +""" +```` +function H1Pk_to_HDIVRT0_interpolator( + VRT0::FESpace, + V1::FESpace +) +```` + +Computes the interpolation matrix from the [`H1Pk`](@ref)-conforming source space `V1` +into the [`HDIVRT0`](@ref)-conforming target space `VRT0` *defined on the same grid* and +returns its sparse matrix representation. + +# Arguments +- `VRT0::FESpace`: [`HDIVRT0`](@ref) target space into which the interpolation is directed. +- `V1::FESpace`: [`H1Pk`](@ref) source space from which the interpolant is taken. + +""" +function H1Pk_to_HDIVRT0_interpolator(VRT0::FESpace, V1::FESpace) + FEType = get_FEType(V1) + facedofs_V1::VariableTargetAdjacency{Int32} = V1[FaceDofs] + dim = get_ncomponents(V1) + EGF = dim == 2 ? Edge1D : Triangle2D + ndofs = get_ndofs(ON_FACES, FEType, EGF) + order = get_polynomialorder(FEType, EGF) + face_basis = get_basis(ON_FACES, FEType, EGF) + result = zeros(ndofs, dim) + ref_integrate!(result, EGF, order, face_basis) + + F0 = FEMatrix(V1, VRT0) + FE::ExtendableSparseMatrix{Float64, Int64} = F0.entries + xgrid = V1.xgrid + @assert xgrid == VRT0.xgrid + xFaceVolumes = xgrid[FaceVolumes] + xFaceNormals = xgrid[FaceNormals] + nfaces = num_sources(xFaceNormals) + for face in 1:nfaces + for dof1 in 1:ndofs + FE[facedofs_V1[dof1, face], face] = dot(view(result, dof1, :), view(xFaceNormals, :, face)) * xFaceVolumes[face] + end + end + flush!(FE) + return F0 +end diff --git a/src/interpolations.jl b/src/interpolations.jl index f1582b1..f43c8ea 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -91,7 +91,6 @@ function ExtendableGrids.interpolate!(target::FEVectorBlock, source; kwargs...) return interpolate!(target, ON_CELLS, source; kwargs...) end - """ ```` function nodevalues_subset!( diff --git a/src/interpolators.jl b/src/interpolators.jl index 1a56e83..d4a7c8c 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -65,6 +65,9 @@ function NodalInterpolator( xCellDofs = FES[CellDofs] ncells = num_cells(grid) function evaluate_broken!(target, exact_function!, items; time = 0, params = [], kwargs...) + if !(eltype(target) <: T) + result = zeros(eltype(target), ncomponents) + end QP.time = time QP.params = params === nothing ? [] : params if isempty(items) @@ -96,6 +99,9 @@ function NodalInterpolator( nnodes = num_nodes(grid) xNodeCells = atranspose(grid[CellNodes]) function evaluate!(target, exact_function!, items; time = 0, params = [], kwargs...) + if !(eltype(target) <: T) + result = zeros(eltype(target), ncomponents) + end QP.time = time QP.params = params === nothing ? [] : params if isempty(items) @@ -370,6 +376,10 @@ function MomentInterpolator( interiordofs = zeros(Int, length(idofs)) function assembly_loop!(target, f_moments, items, exact_function!, QF, L2G, FEB, FEB_moments) + if !(eltype(target) <: Tv) + result_f = zeros(eltype(target), ncomponents) + f_moments = zeros(eltype(target), nmoments) + end weights, xref = QF.w, QF.xref nweights = length(weights) for item::Int in items @@ -582,6 +592,10 @@ function FunctionalInterpolator( FEB = FEEvaluator(FE, operator, QF; AT = AT, T = Tv, L2G = L2G) function assembly_loop!(target, f_fluxes, items, exact_function!, QF, L2G, FEB) + if !(eltype(target) <: Tv) + result_f = zeros(eltype(target), ncomponents) + f_fluxes = zeros(eltype(target), nfluxes) + end weights, xref = QF.w, QF.xref nweights = length(weights) for item::Int in items diff --git a/src/lazy_interpolate.jl b/src/lazy_interpolate.jl index d85b35a..9563e66 100644 --- a/src/lazy_interpolate.jl +++ b/src/lazy_interpolate.jl @@ -79,7 +79,7 @@ function lazy_interpolate!( if xdim_source != xdim_target @assert xtrafo !== nothing "grids have different coordinate dimensions, need xtrafo!" end - PE = PointEvaluator(postprocess, operators, source) + PE = PointEvaluator(postprocess, operators, source; TCoeff = T1) xref = zeros(Tv, xdim_source) x_source = zeros(Tv, xdim_source) cell::Int = start_cell diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl index f484575..0715dda 100644 --- a/src/point_evaluator.jl +++ b/src/point_evaluator.jl @@ -1,4 +1,4 @@ -mutable struct PointEvaluator{Tv <: Real, UT, KFT <: Function} +mutable struct PointEvaluator{Tv <: Real, TCoeff <: Real, UT, KFT <: Function} u_args::Array{UT, 1} ops_args::Array{DataType, 1} kernel::KFT @@ -61,11 +61,11 @@ $(_myprint(default_peval_kwargs())) After construction, call `initialize!` to prepare the evaluator for a given solution, then use `evaluate!` or `evaluate_bary!` to perform point evaluations. """ -function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, kwargs...) +function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, TCoeff = Float64, kwargs...) parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_peval_kwargs()) _update_params!(parameters, kwargs) @assert length(u_args) == length(ops_args) - PE = PointEvaluator{Tv, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, 1, nothing, nothing, nothing, zeros(Tv, 2), parameters) + PE = PointEvaluator{Tv, TCoeff, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, 1, nothing, nothing, nothing, zeros(Tv, 2), parameters) if sol !== nothing initialize!(PE, sol) end @@ -121,7 +121,7 @@ $(_myprint(default_peval_kwargs())) # Notes - This function must be called before using `evaluate!` or `evaluate_bary!` with the `PointEvaluator`. """ -function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {T, UT} +function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...) where {T, TCoeff, UT} _update_params!(O.parameters, kwargs) if UT <: Integer ind_args = O.u_args @@ -159,7 +159,7 @@ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where { op_lengths_args = [size(O.BE_args[1][j].cvals, 1) for j in 1:nargs] op_offsets_args = [0] append!(op_offsets_args, cumsum(op_lengths_args)) - input_args = zeros(T, op_offsets_args[end]) + input_args = zeros(TCoeff, op_offsets_args[end]) FEATs_args = [EffAT4AssemblyType(get_AT(FES_args[j]), AT) for j in 1:nargs] itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j in 1:nargs] diff --git a/test/Project.toml b/test/Project.toml index 2f754bc..7c4879b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,13 +1,12 @@ [deps] -ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c" -ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" -ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" +ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] diff --git a/test/runtests.jl b/test/runtests.jl index 9b85cd5..b5b793c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using Test using ExtendableGrids using ExtendableFEMBase +using ExtendableSparse using ExplicitImports using ExampleJuggler using SparseArrays @@ -28,6 +29,7 @@ end include("test_quadrature.jl") include("test_interpolators.jl") +include("test_interpolation_matrix.jl") include("test_operators.jl") include("test_febasis.jl") include("test_segmentintegrator.jl") @@ -150,6 +152,8 @@ function run_all_tests() run_operator_tests() run_quadrature_tests() run_interpolator_tests() + run_grid_interpolation_matrix_tests() + run_space_interpolation_matrix_tests() run_segmentintegrator_tests() run_pointevaluator_tests() end diff --git a/test/test_interpolation_matrix.jl b/test/test_interpolation_matrix.jl new file mode 100644 index 0000000..d7775f3 --- /dev/null +++ b/test/test_interpolation_matrix.jl @@ -0,0 +1,228 @@ +function run_grid_interpolation_matrix_tests() + # list of FETypes that should be tested + TestCatalog1D = [ + L2P0{1} => 0, + H1P1{1} => 1, + H1P2{1, 1} => 2, + H1P3{1, 1} => 3, + H1Pk{1, 1, 3} => 3, + H1Pk{1, 1, 4} => 4, + H1Pk{1, 1, 5} => 5, + ] + + TestCatalog2D = [ + HCURLN0{2} => 0, + HCURLN1{2} => 1, + HDIVRT0{2} => 0, + HDIVRTk{2, 0} => 0, + HDIVBDM1{2} => 1, + HDIVRT1{2} => 1, + HDIVRTk{2, 1} => 1, + HDIVBDM2{2} => 2, + HDIVRTk{2, 2} => 2, + HDIVRTk{2, 3} => 3, + HDIVRTk{2, 4} => 4, + L2P0{2} => 0, + L2P1{2} => 1, + H1P1{2} => 1, + H1Q1{2} => 1, + H1CR{2} => 1, + H1MINI{2, 2} => 1, + H1P1TEB{2} => 1, + H1BR{2} => 1, + H1P2{2, 2} => 2, + H1P2B{2, 2} => 2, + H1Q2{2, 2} => 2, + H1P3{2, 2} => 3, + H1Pk{2, 2, 3} => 3, + H1Pk{2, 2, 4} => 4, + H1Pk{2, 2, 5} => 5, + ] + + TestCatalog3D = [ + HCURLN0{3} => 0, + HDIVRT0{3} => 0, + HDIVBDM1{3} => 1, + HDIVRT1{3} => 1, + L2P0{3} => 0, + H1P1{3} => 1, + H1Q1{3} => 1, + H1CR{3} => 1, + H1MINI{3, 3} => 1, + H1P1TEB{3} => 1, + H1BR{3} => 1, + H1P2{3, 3} => 2, + H1P3{3, 3} => 3, + ] + + # test interpolation of same space between refined grids + function test_grid_matrix_computation(xgrid, FEType, order; broken::Bool = false, use_cellparents::Bool = false) + u, ~ = exact_function(Val(dim_grid(xgrid)), order) + + source_FES = FESpace{FEType}(xgrid; broken) + target_FES = FESpace{FEType}(uniform_refine(xgrid); broken) + + source_vector = FEVector(source_FES) + interpolate!(source_vector[1], u; bonus_quadorder = order) + + target_vector = FEVector(target_FES) + interpolate!(target_vector[1], u; bonus_quadorder = order) + + interpolation_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents) + matrix_interpolated_entries = interpolation_matrix * source_vector.entries + + return @test norm(target_vector.entries - matrix_interpolated_entries) < tolerance + end + + @testset "Grid Interpolation Matrix Tests" begin + println("\n") + println("============================") + println("Testing Grid Interpolation Matrices in 1D") + println("============================") + xgrid = testgrid(Edge1D) + for (element, order) in TestCatalog1D + @info "Element: ($(element), $(order)) \n" + test_grid_matrix_computation(xgrid, element, order; broken = false) + test_grid_matrix_computation(xgrid, element, order; broken = true) + end + + println("\n") + println("============================") + println("Testing Grid Interpolation Matrices in 2D") + println("============================") + for EG in [Triangle2D, Parallelogram2D] + xgrid = uniform_refine(reference_domain(EG), 1) + for (element, order) in TestCatalog2D, broken in (false, true) + @info "Element: ($(EG), $(element), $(order), broken = $(broken)) \n" + if ExtendableFEMBase.isdefined(element, EG, broken) + test_grid_matrix_computation(xgrid, element, order; broken) + else + @warn "($(element),$(order)) (broken = $(broken)) not defined on $(EG) (skipping test case)" + end + end + end + + println("\n") + println("============================") + println("Testing Grid Interpolation Matrices in 3D") + println("============================") + for EG in [Tetrahedron3D, Parallelepiped3D] + xgrid = uniform_refine(reference_domain(EG), 1) + for (element, order) in TestCatalog3D, broken in (false, true) + @info "Element: ($(EG), $(element), $(order), broken = $(broken)) \n" + if ExtendableFEMBase.isdefined(element, EG, broken) + test_grid_matrix_computation(xgrid, element, order; broken) + else + @warn "($(element),$(order)) (broken = $(broken)) not defined on $(EG) (skipping test case)" + end + end + end + end + + return println("") +end + +function run_space_interpolation_matrix_tests() + # list of space pairs that should be interpolated into one another + + PairTestCatalog1D = [ + (L2P0{1}, H1P1{1}) => 0, + (H1P1{1}, H1P2{1, 1}) => 1, + (H1P2{1, 1}, H1P3{1, 1}) => 2, + (H1Pk{1, 1, 3}, H1Pk{1, 1, 5}) => 3, + ] + + PairTestCatalog2D = [ + (H1P1{2}, HDIVRT0{2}) => 0, + (H1P2{2, 2}, HDIVRT0{2}) => 0, + (H1P2{2, 2}, HDIVRT1{2}) => 1, + (H1P2{2, 2}, HDIVRTk{2, 2}) => 2, + (HDIVRT1{2}, HDIVBDM1{2}) => 1, + (L2P0{2}, L2P1{2}) => 0, + (H1P2B{2, 2}, H1BR{2}) => 1, + ] + + PairTestCatalog3D = [ + (H1P1{3}, HDIVRT0{3}) => 0, + (H1P2{3, 3}, HDIVRT0{3}) => 0, + (H1P2{3, 3}, HDIVRT1{3}) => 1, + (HDIVRT1{3}, HDIVBDM1{3}) => 1, + (L2P0{3}, L2P1{3}) => 0, + (H1P2B{3, 3}, H1BR{3}) => 1, + ] + + # test interpolation for different elements on same grid + function test_space_matrix_computation(xgrid, source_FEType, target_FEType, order; broken::Bool = false, use_cellparents::Bool = false) + u, ~ = exact_function(Val(dim_grid(xgrid)), order) + + source_FES = FESpace{source_FEType}(xgrid) + target_FES = FESpace{target_FEType}(xgrid; broken) + + source_vector = FEVector(source_FES) + interpolate!(source_vector[1], u; bonus_quadorder = order) + + target_vector = FEVector(target_FES) + interpolate!(target_vector[1], u; bonus_quadorder = order) + + interpolation_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents) + matrix_interpolated_entries = interpolation_matrix * source_vector.entries + + return @test norm(target_vector.entries - matrix_interpolated_entries) < tolerance + end + + @testset "Space Interpolation Matrix Tests" begin + println("\n") + println("============================") + println("Testing Space Interpolation Matrices in 1D") + println("============================") + xgrid = testgrid(Edge1D) + for ((source_element, target_element), order) in PairTestCatalog1D + @info "Element pair: ($(source_element), $(target_element)), order: $(order) \n" + test_space_matrix_computation(xgrid, source_element, target_element, order; broken = false) + test_space_matrix_computation(xgrid, source_element, target_element, order; broken = true) + end + + println("\n") + println("============================") + println("Testing Space Interpolation Matrices in 2D") + println("============================") + for EG in [Triangle2D, Parallelogram2D] + xgrid = uniform_refine(reference_domain(EG), 1) + for ((source_element, target_element), order) in PairTestCatalog2D, broken in (false, true) + @info "Element pair: ($(EG), $(source_element), $(target_element)), order: $(order), broken = $(broken) \n" + if ExtendableFEMBase.isdefined(target_element, EG, broken) && ExtendableFEMBase.isdefined(source_element, EG, broken) + test_space_matrix_computation(xgrid, source_element, target_element, order; broken) + + if (source_element, target_element) == (H1P1{2}, HDIVRT0{2}) && !broken + source_FES = FESpace{source_element}(xgrid; broken) + target_FES = FESpace{target_element}(xgrid; broken) + autodiff_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents = false) + RT0_matrix = H1Pk_to_HDIVRT0_interpolator(target_FES, source_FES) + + @test norm(autodiff_matrix' - RT0_matrix[1]) < tolerance + end + else + @warn "($(target_element),$(order)) (broken = $(broken)) not defined on $(EG) (skipping test case)" + end + end + end + + println("\n") + println("============================") + println("Testing Space Interpolation Matrices in 3D") + println("============================") + for EG in [Tetrahedron3D, Parallelepiped3D] + xgrid = uniform_refine(reference_domain(EG), 1) + for ((source_element, target_element), order) in PairTestCatalog3D, broken in (false, true) + @info "Element pair: ($(EG), $(source_element), $(target_element)), order: $(order), broken =$(broken) \n" + if ExtendableFEMBase.isdefined(target_element, EG, broken) && ExtendableFEMBase.isdefined(source_element, EG, broken) + test_space_matrix_computation(xgrid, source_element, target_element, order; broken) + else + @warn "($(target_element),$(order)) (broken = $(broken)) not defined on $(EG) (skipping test case)" + end + end + end + end + + return println("") +end diff --git a/test/test_interpolators.jl b/test/test_interpolators.jl index c8ef9b1..7cc6394 100644 --- a/test/test_interpolators.jl +++ b/test/test_interpolators.jl @@ -109,7 +109,7 @@ function run_interpolator_tests() show(devnull, Solution) # compute error - error = compute_error(Solution[1], u, order) + error = compute_error(Solution[1], u, order + 1) println("FEType = $FEType $(broken ? "broken" : "") $AT | ndofs = $(FES.ndofs) | order = $order | error = $(norm(error, Inf))") return @test norm(error) < tolerance end