From 1f3c5c1545b913d72a1f49ce83b35c997640bf55 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Thu, 30 Oct 2025 20:32:30 +0100 Subject: [PATCH 1/8] added function for computing jacobian of interpolation between finite element spaces --- Project.toml | 8 +- src/ExtendableFEMBase.jl | 4 + src/interpolations.jl | 56 ++++++ src/interpolators.jl | 14 ++ src/lazy_interpolate.jl | 2 +- src/point_evaluator.jl | 10 +- test/Project.toml | 9 +- test/runtests.jl | 4 + test/test_interpolation_matrix.jl | 282 ++++++++++++++++++++++++++++++ 9 files changed, 377 insertions(+), 12 deletions(-) create mode 100644 test/test_interpolation_matrix.jl diff --git a/Project.toml b/Project.toml index 29998fc8..d7110c16 100644 --- a/Project.toml +++ b/Project.toml @@ -1,10 +1,11 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" -authors = ["Christian Merdon ", "Patrick Jaap "] version = "1.4.0" +authors = ["Christian Merdon ", "Patrick Jaap "] [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/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index 6d3cef67..6f744825 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") @@ -114,6 +117,7 @@ export interpolate! # must be defined separately by each FEdefinition export nodevalues, continuify export nodevalues!, nodevalues_subset! export nodevalues_view +export compute_interpolation_jacobian export interpolator_matrix diff --git a/src/interpolations.jl b/src/interpolations.jl index f1582b1e..5a030ca5 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -91,6 +91,62 @@ function ExtendableGrids.interpolate!(target::FEVectorBlock, source; kwargs...) return interpolate!(target, ON_CELLS, source; kwargs...) end +""" +```` +function compute_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 +- This function can be used for computing prolongation or restriction operators if the `FESpace`s are defined on coarser/finer grids, respectively. + +""" +function compute_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 + """ ```` diff --git a/src/interpolators.jl b/src/interpolators.jl index 1a56e83f..d4a7c8c5 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 d85b35a6..9563e66b 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 f4845756..0715ddaa 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 2f754bcd..7c4879be 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 9b85cd5a..b5b793c4 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 00000000..f1f36d95 --- /dev/null +++ b/test/test_interpolation_matrix.jl @@ -0,0 +1,282 @@ +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, + ] + + 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, + ] + + 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, + ] + + 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, + ] + + 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, + ] + + 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 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_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, + ] + + + # V1: H1-conforming space, computes matrix into RT0 space + function RT0_interpolator(V1::FESpace, VRT0::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 + 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 + + # 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_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_interpolation_jacobian(target_FES, source_FES; use_cellparents = false) + RT0_matrix = RT0_interpolator(source_FES, target_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 From 115c0afca5691f2f4e57eef39b80fab18b360756 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Thu, 30 Oct 2025 20:34:48 +0100 Subject: [PATCH 2/8] fixed typos in some p1 elements --- src/fedefs/h1_p1.jl | 2 +- src/fedefs/l2_p1.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fedefs/h1_p1.jl b/src/fedefs/h1_p1.jl index e14625e8..192c09a3 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/l2_p1.jl b/src/fedefs/l2_p1.jl index 2fe07f2d..5c67e106 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} From 21dde9e22cbe7e8a4bf1f504d3743ef44e516b24 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Thu, 30 Oct 2025 20:36:42 +0100 Subject: [PATCH 3/8] added missing interpolation calls for Hcurl-conforming elements from @chmerdon --- src/fedefs/hcurl_n1.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/fedefs/hcurl_n1.jl b/src/fedefs/hcurl_n1.jl index b7e3b3a4..53dc7b1f 100644 --- a/src/fedefs/hcurl_n1.jl +++ b/src/fedefs/hcurl_n1.jl @@ -34,6 +34,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} @@ -65,6 +66,10 @@ 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 + get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...) + end # on faces dofs are only tangential fluxes From 2e92ba7e19dfb7e658ead3a82d5e3276f4f400cf Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Thu, 30 Oct 2025 20:37:16 +0100 Subject: [PATCH 4/8] updated docs and version bump --- CHANGELOG.md | 4 ++++ Project.toml | 2 +- docs/Project.toml | 12 ++++++------ docs/src/interpolations.md | 4 ++++ 4 files changed, 15 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ffdd99ee..091f3e96 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # CHANGES +## v1.5.0 October 30, 2025 + - added type parameter to `PointEvaluator` so `lazy_interpolate!` is now accessible for automatic differentiation + - added new `compute_interpolation_jacobian` function for computing Jacobians of interpolations between finite element spaces + ## 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 d7110c16..3780ef90 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" -version = "1.4.0" +version = "1.5.0" authors = ["Christian Merdon ", "Patrick Jaap "] [deps] diff --git a/docs/Project.toml b/docs/Project.toml index c512582b..7a359d3f 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/src/interpolations.md b/docs/src/interpolations.md index 2d867ca9..e638069d 100644 --- a/docs/src/interpolations.md +++ b/docs/src/interpolations.md @@ -37,6 +37,10 @@ Additionally, you can transfer finite element functions from one grid to another lazy_interpolate! ``` +```@docs +compute_interpolation_jacobian +``` + 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). From a2b8d6a914c493942871ab35c1cdd7f3fad3e84b Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Mon, 3 Nov 2025 13:21:59 +0100 Subject: [PATCH 5/8] fixed hcurln1 interpolation and rearranged matrix interpolation code --- src/ExtendableFEMBase.jl | 4 +- src/fedefs/hcurl_n1.jl | 7 +- src/interpolation_matrix_representations.jl | 103 ++++++++++++++++++++ src/interpolations.jl | 57 ----------- test/test_interpolation_matrix.jl | 62 +----------- test/test_interpolators.jl | 2 +- 6 files changed, 115 insertions(+), 120 deletions(-) create mode 100644 src/interpolation_matrix_representations.jl diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index 6f744825..ed179357 100644 --- a/src/ExtendableFEMBase.jl +++ b/src/ExtendableFEMBase.jl @@ -117,7 +117,6 @@ export interpolate! # must be defined separately by each FEdefinition export nodevalues, continuify export nodevalues!, nodevalues_subset! export nodevalues_view -export compute_interpolation_jacobian export interpolator_matrix @@ -175,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/hcurl_n1.jl b/src/fedefs/hcurl_n1.jl index 53dc7b1f..dbfebd01 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] @@ -57,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...) @@ -68,8 +70,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, end # set values of interior N1 functions such that P0 moments are preserved - get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...) - + return get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...) end # on faces dofs are only tangential fluxes diff --git a/src/interpolation_matrix_representations.jl b/src/interpolation_matrix_representations.jl new file mode 100644 index 00000000..0ffd28e1 --- /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 5a030ca5..f43c8ead 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -91,63 +91,6 @@ function ExtendableGrids.interpolate!(target::FEVectorBlock, source; kwargs...) return interpolate!(target, ON_CELLS, source; kwargs...) end -""" -```` -function compute_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 -- This function can be used for computing prolongation or restriction operators if the `FESpace`s are defined on coarser/finer grids, respectively. - -""" -function compute_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 nodevalues_subset!( diff --git a/test/test_interpolation_matrix.jl b/test/test_interpolation_matrix.jl index f1f36d95..d7775f35 100644 --- a/test/test_interpolation_matrix.jl +++ b/test/test_interpolation_matrix.jl @@ -10,13 +10,6 @@ function run_grid_interpolation_matrix_tests() H1Pk{1, 1, 5} => 5, ] - 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, - ] - TestCatalog2D = [ HCURLN0{2} => 0, HCURLN1{2} => 1, @@ -46,16 +39,6 @@ function run_grid_interpolation_matrix_tests() H1Pk{2, 2, 5} => 5, ] - 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, - ] - TestCatalog3D = [ HCURLN0{3} => 0, HDIVRT0{3} => 0, @@ -72,15 +55,6 @@ function run_grid_interpolation_matrix_tests() H1P3{3, 3} => 3, ] - 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 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) @@ -94,7 +68,7 @@ function run_grid_interpolation_matrix_tests() target_vector = FEVector(target_FES) interpolate!(target_vector[1], u; bonus_quadorder = order) - interpolation_matrix = compute_interpolation_jacobian(target_FES, source_FES; use_cellparents) + 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 @@ -177,34 +151,6 @@ function run_space_interpolation_matrix_tests() (H1P2B{3, 3}, H1BR{3}) => 1, ] - - # V1: H1-conforming space, computes matrix into RT0 space - function RT0_interpolator(V1::FESpace, VRT0::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 - 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 - # 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) @@ -218,7 +164,7 @@ function run_space_interpolation_matrix_tests() target_vector = FEVector(target_FES) interpolate!(target_vector[1], u; bonus_quadorder = order) - interpolation_matrix = compute_interpolation_jacobian(target_FES, source_FES; use_cellparents) + 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 @@ -250,8 +196,8 @@ function run_space_interpolation_matrix_tests() 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_interpolation_jacobian(target_FES, source_FES; use_cellparents = false) - RT0_matrix = RT0_interpolator(source_FES, target_FES) + 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 diff --git a/test/test_interpolators.jl b/test/test_interpolators.jl index c8ef9b11..7cc63940 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 From 6251a891b57f0bf70190971f53a0f39b90f8c687 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Mon, 3 Nov 2025 13:23:23 +0100 Subject: [PATCH 6/8] updated docs + included changelog entry --- CHANGELOG.md | 6 ++++-- docs/make.jl | 1 + docs/src/changelog.md | 6 ++++++ docs/src/interpolations.md | 6 +++++- 4 files changed, 16 insertions(+), 3 deletions(-) create mode 100644 docs/src/changelog.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 091f3e96..31497f9e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,10 @@ # CHANGES -## v1.5.0 October 30, 2025 +## v1.5.0 November 03, 2025 - added type parameter to `PointEvaluator` so `lazy_interpolate!` is now accessible for automatic differentiation - - added new `compute_interpolation_jacobian` function for computing Jacobians of interpolations between finite element spaces + - 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}) diff --git a/docs/make.jl b/docs/make.jl index a41f7638..15e7cc82 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 00000000..789e3fc6 --- /dev/null +++ b/docs/src/changelog.md @@ -0,0 +1,6 @@ +```@eval +using Markdown +Markdown.parse(""" +$(read("../../CHANGELOG.md",String)) +""") +``` \ No newline at end of file diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md index e638069d..5b681836 100644 --- a/docs/src/interpolations.md +++ b/docs/src/interpolations.md @@ -38,7 +38,11 @@ lazy_interpolate! ``` ```@docs -compute_interpolation_jacobian +compute_lazy_interpolation_jacobian +``` + +```@docs +H1Pk_to_HDIVRT0_interpolator ``` The following function continuously interpolates finite element function into a H1Pk space by From 3a93140ac00ec8a2f6b5d32636e75ddd29f36100 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Mon, 3 Nov 2025 13:27:25 +0100 Subject: [PATCH 7/8] added newline to changelog.md --- docs/src/changelog.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 789e3fc6..d5dc3f51 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -3,4 +3,4 @@ using Markdown Markdown.parse(""" $(read("../../CHANGELOG.md",String)) """) -``` \ No newline at end of file +``` From 94081918bea6364d553efacb413046cc3b5ac126 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Fri, 14 Nov 2025 11:27:42 +0100 Subject: [PATCH 8/8] modified Ex290 for ringsector --- .../Example290_InterpolationBetweenMeshes.jl | 24 ++++++++++++------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/examples/Example290_InterpolationBetweenMeshes.jl b/examples/Example290_InterpolationBetweenMeshes.jl index bb3a3e48..72e19608 100644 --- a/examples/Example290_InterpolationBetweenMeshes.jl +++ b/examples/Example290_InterpolationBetweenMeshes.jl @@ -15,7 +15,7 @@ module Example290_InterpolationBetweenMeshes using ExtendableFEMBase using ExtendableGrids -using GridVisualize +#using GridVisualize ## function to interpolate function u!(result, qpinfo) @@ -25,10 +25,13 @@ function u!(result, qpinfo) end ## everything is wrapped in a main function -function main(; ν = 1.0e-3, nrefs = 4, Plotter = nothing) +## crashes with nrefs=0 +#xgrid = ringsector(0.5:0.1:1., linspace(0.,2π,10);eltype = Triangle2D); xgrid = uniform_refine(xgrid, 0); @show xgrid[CellNodes];source_space = FESpace{H1Pk{2,2,1}}(xgrid); @show source_space.dofgrid[CellFaces] ; target_space = FESpace{H1Pk{2,2,1}}(uniform_refine(xgrid,1; store_parents = true)); source_fevector = FEVector(source_space); interpolate!(source_fevector[1], (result,qpinfo) -> (result .= qpinfo.x)); target_fevector = FEVector(target_space); +function main(; ν = 1.0e-3, nrefs = 0, Plotter = nothing) ## generate two grids - xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), nrefs) + xgrid1 = uniform_refine(ringsector(0.5:0.1:1., linspace(0.,2π,10); eltype = Triangle2D), nrefs) + #xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), nrefs) xgrid2 = uniform_refine(xgrid1, 3; store_parents = true) @show xgrid1 xgrid2 @@ -44,19 +47,22 @@ function main(; ν = 1.0e-3, nrefs = 4, Plotter = nothing) FEFunction2 = FEVector(FES2) ## interpolate function onto first grid + # breaks for ringsector at nrefs = 0, but not at nrefs ≥ 1 @time interpolate!(FEFunction1[1], u!) @time interpolate!(FEFunction2[1], u!) ## interpolate onto other grid - @time lazy_interpolate!(FEFunction2[1], FEFunction1) - @time lazy_interpolate!(FEFunction2[1], FEFunction1; use_cellparents = true) + #@time lazy_interpolate!(FEFunction2[1], FEFunction1) + + # freezes up with ringsector + #@time lazy_interpolate!(FEFunction2[1], FEFunction1; use_cellparents = true) ## plot - p = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400)) - scalarplot!(p[1, 1], xgrid1, view(nodevalues(FEFunction1[1]), 1, :), levels = 11, title = "u_h ($FEType1, coarse grid)") - scalarplot!(p[1, 2], xgrid2, view(nodevalues(FEFunction2[1]), 1, :), levels = 11, title = "u_h ($FEType2, fine grid)") + # p = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400)) + # scalarplot!(p[1, 1], xgrid1, view(nodevalues(FEFunction1[1]), 1, :), levels = 11, title = "u_h ($FEType1, coarse grid)") + # scalarplot!(p[1, 2], xgrid2, view(nodevalues(FEFunction2[1]), 1, :), levels = 11, title = "u_h ($FEType2, fine grid)") - return p + #return p end function generateplots(dir = pwd(); Plotter = nothing, kwargs...)