WIAS-PDELib · Da-Be-Ru · Nov 14, 2025 · Oct 30, 2025 · Oct 30, 2025 · Oct 30, 2025
diff --git 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
@@ -1,10 +1,11 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
+version = "1.5.0"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
-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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -46,13 +46,16 @@ 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
 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
@@ -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
@@ -27,13 +27,16 @@ 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]
     result[2] = result[1] * (qpinfo.xref[1] - 1 // 2)
     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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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]
-Original file line number
+Diff line change
@@ Expand Up @@
         return interpolate!(target, ON_CELLS, source; kwargs...)
     end
     """
     ````
     function nodevalues_subset!(
@@ Expand Down @@