Bug: ringsector breaks interpolation call

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

Project.toml

@@ -1,10 +1,11 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
+version = "1.5.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
-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"

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"

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[

docs/src/changelog.md

@@ -0,0 +1,6 @@
+```@eval
+using Markdown
+Markdown.parse("""
+$(read("../../CHANGELOG.md",String))
+""")
+```

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).
 

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...)

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
 

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}

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

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}

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

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!(

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