restructuring of Interpolators (#30)

* completely replaced old point_evaluation!, point_evaluation_broken! by NodalInterpolator
and ensure_moments! by MomentInterpolator and many custom interpolatioon
in non-H1 finite elements by a common new FunctionalInterpolator structure;
now every finite element interpolation is a composition of these three
interpolating structures
* better structured interpolation tests

---------

Co-authored-by: Patrick Jaap <patrick.jaap@wias-berlin.de>

Author: chmerdon

CHANGELOG.md

@@ -1,6 +1,11 @@
 # CHANGES
 
-## v1.0.0
+## v1.1.0 May 15, 2025
+  - overhaul of interpolators: each finite element
+    interpolation is composed of the three new
+    structures NodalInterpolator, MomentInterpolator, FunctionalInterpolator
+  
+## v1.0.0 April 7, 2025
   - doc improvements
   - several bugfixes regarding finite elements defined on subgrids
   - integrate now reacts on regions parameter and has proper docstrings

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-authors = ["Christian Merdon <merdon@wias-berlin.de>"]
-version = "1.0.0"
+authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
+version = "1.1.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
@@ -24,7 +24,7 @@ ExtendableFEMBaseUnicodePlotsExt = ["UnicodePlots"]
 [compat]
 DiffResults = "1"
 DocStringExtensions = "0.8,0.9"
-ExtendableGrids = "1.9.2"
+ExtendableGrids = "1.13.0"
 ExtendableSparse = "1.5.1"
 ForwardDiff = "0.10.35, 1"
 LinearAlgebra = "1.9"

examples/Example281_DiscontinuousPlot.jl

@@ -61,9 +61,11 @@ function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing)
     nodevals4nodes2 = nodevalues(FEFunction[1], Identity; abs = abs, regions = [2], nodes = subnodes2)
 
 	## plot
-	p = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 500))
-    gridplot!(p[1,1], xgrid)
-	scalarplot!(p[1, 2], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1,:), view(nodevals4nodes2,:)], cellwise = false, levels = 11, title = "u")
+    if Plotter !== nothing
+	    p = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 500))
+        gridplot!(p[1,1], xgrid)
+	    scalarplot!(p[1, 2], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1,:), view(nodevals4nodes2,:)], cellwise = false, levels = 11, title = "u")
+    end
 
     return p
 end

src/ExtendableFEMBase.jl

@@ -19,10 +19,11 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
     CellFinder, CellGeometries, CellNodes, CellParents,
     CellRegions, CellVolumes, Coordinates, Edge1D,
     EdgeCells, EdgeGeometries, EdgeNodes, EdgeTangents,
-    EdgeVolumes, ElementGeometries, ExtendableGrid,
+    EdgeParents, EdgeVolumes, ElementGeometries, ExtendableGrid,
     FaceCells, FaceEdgeSigns, FaceEdges, FaceGeometries,
     FaceNodes, FaceNormals, FaceParents, FaceRegions,
     FaceVolumes, GridComponentGeometries4AssemblyType,
+    NodeCells,
     GridComponentNodes4AssemblyType,
     GridComponentRegions4AssemblyType,
     GridComponentUniqueGeometries4AssemblyType,
@@ -107,6 +108,8 @@ export get_polynomialorder, get_ndofs, get_ndofs_all
 export get_ncomponents, get_edim
 export get_basis, get_coefficients, get_basissubset
 
+export NodalInterpolator
+
 export interpolate! # must be defined separately by each FEdefinition
 export nodevalues, continuify
 export nodevalues!, nodevalues_subset!

src/fedefs/h1_bubble.jl

@@ -35,25 +35,10 @@ isdefined(FEType::Type{<:H1BUBBLE}, ::Type{<:Tetrahedron3D}) = true
 
 interior_dofs_offset(::Type{ON_CELLS}, ::Type{H1BUBBLE{ncomponents}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents} = 0
 
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: H1BUBBLE, APT} = MomentInterpolator(FES, ON_CELLS)
+
 function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: H1BUBBLE, APT}
-    xCellVolumes = FE.dofgrid[CellVolumes]
-    ncells = num_sources(FE.dofgrid[CellNodes])
-    if items == []
-        items = 1:ncells
-    else
-        items = filter(!iszero, items)
-    end
-    ncomponents = get_ncomponents(FEType)
-    integrals4cell = zeros(T, ncomponents, ncells)
-    integrate!(integrals4cell, FE.dofgrid, ON_CELLS, exact_function!; items = items, kwargs...)
-    for cell in items
-        if cell != 0
-            for c in 1:ncomponents
-                Target[(cell - 1) * ncomponents + c] = integrals4cell[c, cell] / xCellVolumes[cell]
-            end
-        end
-    end
-    return
+    get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items, kwargs...)
 end
 
 function get_basis(::Type{<:AssemblyType}, ::Type{H1BUBBLE{ncomponents}}, ::Type{<:AbstractElementGeometry1D}) where {ncomponents}

src/fedefs/h1_mini.jl

@@ -42,10 +42,11 @@ get_ref_cellmoments(::Type{<:H1MINI}, ::Type{<:Quadrilateral2D}) = [1 // 4, 1 //
 
 interior_dofs_offset(::Type{ON_CELLS}, ::Type{H1MINI{ncomponents, edim}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents, edim} = num_nodes(EG)
 
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}) where {Tv, Ti, FEType <: H1MINI, APT} = NodalInterpolator(FES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: H1MINI, APT} = MomentInterpolator(FES, ON_CELLS)
+
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
-    nnodes = size(FE.dofgrid[Coordinates], 2)
-    ncells = num_sources(FE.dofgrid[CellNodes])
-    return point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = nnodes + ncells, kwargs...)
+    return get_interpolator(FE, AT_NODES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1MINI, APT}
@@ -66,7 +67,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
     interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
     # fix cell bubble value by preserving integral mean
-    return ensure_moments!(Target, FE, ON_CELLS, exact_function!; items = items, kwargs...)
+    return get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function nodevalues!(Target::AbstractArray{<:Real, 2}, Source::AbstractArray{<:Real, 1}, FE::FESpace{<:H1MINI})

src/fedefs/h1_p1.jl

@@ -32,31 +32,32 @@ isdefined(FEType::Type{<:H1P1}, ::Type{<:AbstractElementGeometry1D}) = true
 isdefined(FEType::Type{<:H1P1}, ::Type{<:Triangle2D}) = true
 isdefined(FEType::Type{<:H1P1}, ::Type{<:Tetrahedron3D}) = true
 
-function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
-    nnodes = size(FE.dofgrid[Coordinates], 2)
-    return point_evaluation!(Target, FE, AT_NODES, exact_function; items = items, component_offset = nnodes, kwargs...)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}) where {Tv, Ti, FEType <: H1P1, APT} = NodalInterpolator(FES)
+
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
+    return get_interpolator(FE, AT_NODES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
-function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
+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 interpolate!(Target, FE, AT_NODES, exact_function; items = subitems, kwargs...)
+    return interpolaste!(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}
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
     # delegate face nodes to node interpolation
     subitems = slice(FE.dofgrid[FaceNodes], items)
-    return interpolate!(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_CELLS}, exact_function; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P1, APT}
     return if FE.broken == true
-        # broken interpolation
-        point_evaluation_broken!(Target, FE, ON_CELLS, exact_function; items = items, time = time)
+        # broken nodal interpolation piecewise on cells
+        interpolate!(Target, FE, AT_NODES, exact_function!; items = items, kwargs...)
     else
         # delegate cell nodes to node interpolation
         subitems = slice(FE.dofgrid[CellNodes], items)
-        interpolate!(Target, FE, AT_NODES, exact_function; items = subitems, kwargs...)
+        interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
     end
 end
 

src/fedefs/h1_p2.jl

@@ -20,10 +20,10 @@ end
 get_ncomponents(FEType::Type{<:H1P2}) = FEType.parameters[1]
 get_edim(FEType::Type{<:H1P2}) = FEType.parameters[2]
 
-get_ndofs(::Type{<:AssemblyType}, FEType::Type{<:H1P2}, EG::Type{<:AbstractElementGeometry0D}) = FEType.parameters[1]
-get_ndofs(::Union{Type{<:ON_EDGES}, Type{<:ON_BEDGES}}, FEType::Type{<:H1P2}, EG::Type{<:AbstractElementGeometry1D}) = 3 * FEType.parameters[1]
-get_ndofs(::Union{Type{<:ON_FACES}, Type{<:ON_BFACES}}, FEType::Type{<:H1P2}, EG::Type{<:Union{AbstractElementGeometry1D, Triangle2D, Tetrahedron3D}}) = Int((FEType.parameters[2]) * (FEType.parameters[2] + 1) / 2 * FEType.parameters[1])
-get_ndofs(::Type{<:ON_CELLS}, FEType::Type{<:H1P2}, EG::Type{<:Union{AbstractElementGeometry1D, Triangle2D, Tetrahedron3D}}) = Int((FEType.parameters[2] + 1) * (FEType.parameters[2] + 2) / 2 * FEType.parameters[1])
+get_ndofs(::Type{<:AssemblyType}, FEType::Type{H1P2{n, e}}, EG::Type{<:AbstractElementGeometry0D}) where {n, e} = n
+get_ndofs(::Type{<:AssemblyType}, FEType::Type{H1P2{n, e}}, EG::Type{<:AbstractElementGeometry1D}) where {n, e} = 3 * n
+get_ndofs(::Type{<:AssemblyType}, FEType::Type{H1P2{n, e}}, EG::Type{<:Triangle2D}) where {n, e} = 6 * n
+get_ndofs(::Type{<:AssemblyType}, FEType::Type{H1P2{n, e}}, EG::Type{<:Tetrahedron3D}) where {n, e} = 10 * n
 
 get_polynomialorder(::Type{<:H1P2}, ::Type{<:Edge1D}) = 2;
 get_polynomialorder(::Type{<:H1P2}, ::Type{<:Triangle2D}) = 2;
@@ -43,18 +43,13 @@ isdefined(FEType::Type{<:H1P2}, ::Type{<:Tetrahedron3D}) = true
 
 interior_dofs_offset(::Type{<:AssemblyType}, ::Type{H1P2{ncomponents, edim}}, ::Type{Edge1D}) where {ncomponents, edim} = 2
 
-function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2, APT}
-    edim = get_edim(FEType)
-    nnodes = size(FE.dofgrid[Coordinates], 2)
-    if edim == 1
-        nedges = num_sources(FE.dofgrid[CellNodes])
-    elseif edim == 2
-        nedges = num_sources(FE.dofgrid[FaceNodes])
-    elseif edim == 3
-        nedges = num_sources(FE.dofgrid[EdgeNodes])
-    end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}) where {Tv, Ti, FEType <: H1P2, APT} = NodalInterpolator(FES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: H1P2, APT} = MomentInterpolator(FES, ON_FACES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}) where {Tv, Ti, FEType <: H1P2, APT} = MomentInterpolator(FES, ON_EDGES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: H1P2, APT} = MomentInterpolator(FES, ON_CELLS)
 
-    return point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = nnodes + nedges, kwargs...)
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2, APT}
+    return get_interpolator(FE, AT_NODES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2, APT}
@@ -65,7 +60,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
         interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
         # perform edge mean interpolation
-        ensure_moments!(Target, FE, ON_EDGES, exact_function!; items = items, kwargs...)
+        get_interpolator(FE, ON_EDGES).evaluate!(Target, exact_function!, items)
     end
 end
 
@@ -77,7 +72,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
         interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
         # perform face mean interpolation
-        ensure_moments!(Target, FE, ON_FACES, exact_function!; items = items, kwargs...)
+        get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
     elseif edim == 3
         # delegate face edges to edge interpolation
         subitems = slice(FE.dofgrid[FaceEdges], items)
@@ -106,7 +101,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
         interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
         # preserve cell integral
-        ensure_moments!(Target, FE, ON_CELLS, exact_function!; items = items, kwargs...)
+        get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...)
     end
 end
 

src/fedefs/h1_p2b.jl

@@ -35,17 +35,12 @@ interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{H1P2B{ncomponents, edim}}, ::Typ
 
 get_ref_cellmoments(::Type{<:H1P2B}, ::Type{<:Triangle2D}) = [0 // 1, 0 // 1, 0 // 1, 1 // 3, 1 // 3, 1 // 3, 1 // 1] # integrals of 1D basis functions over reference cell (divided by volume)
 
-function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2B, APT}
-    edim = get_edim(FEType)
-    nnodes = size(FE.dofgrid[Coordinates], 2)
-    offset = nnodes + num_sources(FE.dofgrid[CellNodes])
-    if edim == 2
-        offset += num_sources(FE.dofgrid[FaceNodes])
-    elseif edim == 3
-        offset += num_sources(FE.dofgrid[EdgeNodes])
-    end
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}) where {Tv, Ti, FEType <: H1P2B, APT} = NodalInterpolator(FES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: H1P2B, APT} = MomentInterpolator(FES, ON_FACES)
+init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: H1P2B, APT} = MomentInterpolator(FES, ON_CELLS)
 
-    return point_evaluation!(Target, FE, AT_NODES, exact_function!; items = items, component_offset = offset, kwargs...)
+function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{AT_NODES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2B, APT}
+    return get_interpolator(FE, AT_NODES).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_EDGES}, exact_function!; items = [], kwargs...) where {Tv, Ti, FEType <: H1P2B, APT}
@@ -56,7 +51,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
         interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
         # perform edge mean interpolation
-        ensure_moments!(Target, FE, ON_EDGES, exact_function!; items = items, kwargs...)
+        get_interpolator(FE, ON_EDGES).evaluate!(Target, exact_function!, items)
     end
 end
 
@@ -68,7 +63,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
         interpolate!(Target, FE, AT_NODES, exact_function!; items = subitems, kwargs...)
 
         # perform face mean interpolation
-        ensure_moments!(Target, FE, ON_FACES, exact_function!; items = items, kwargs...)
+        get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...)
     elseif edim == 3
         # delegate face edges to edge interpolation
         subitems = slice(FE.dofgrid[FaceEdges], items)
@@ -100,7 +95,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
     end
 
     # fix cell bubble value by preserving integral mean
-    return ensure_moments!(Target, FE, ON_CELLS, exact_function!; items = items, kwargs...)
+    get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...)
 end
 
 function get_basis(AT::Union{Type{<:ON_FACES}, Type{<:ON_BFACES}}, ::Type{H1P2B{ncomponents, edim}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents, edim}