Fix cellparents interpolation (#54)

* extracted incident cell of items into qpinfo in interpolators for use_cellparents in lazy_interpolate!

* correct itemcell fill for various assembly type delegations in interpolations

* adjusted itemcells computation for ON_EDGES assembly in MomentInterpolator

Author: Da-Be-Ru

CHANGELOG.md

@@ -1,10 +1,11 @@
 # CHANGES
 
-## v1.5.0 November 03, 2025
+## v1.5.0 November 21, 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
+  - updated `MomentInterpolator` to deliver incident cells to `QPInfos` to speed up `lazy_interpolate!` 
 
 ## v1.4.0 July 17, 2025
   - added new weighted reconstruction operator `WeightedReconstruct` (so far only for H1BR{2})

src/interpolators.jl

@@ -191,19 +191,29 @@ function MomentInterpolator(
     itemregions = xgrid[GridComponentRegions4AssemblyType(AT)]
     itemdofs = Dofmap4AssemblyType(FE, AT)
     has_normals = true
+    nitems::Int = num_sources(itemnodes)
     if AT <: ON_FACES
         itemnormals = xgrid[FaceNormals]
+        @views itemcells = xgrid[FaceCells][1, :]
     elseif AT <: ON_BFACES
         itemnormals = xgrid[FaceNormals][:, xgrid[BFaceFaces]]
+        itemcells = xgrid[BFaceCells]
     else
         has_normals = false
+        if AT <: ON_EDGES # e.g. H1P2 delegation for tetrahedron to edges
+            ec = xgrid[EdgeCells]
+            # extract first cell to which a given edge belongs for every edge
+            itemcells = view(ec.colentries, view(ec.colstart, 1:nitems))
+        else
+            # leave empty for delegations to cell dofs, fill up with items at call site
+            itemcells = []
+        end
     end
     EGs = xgrid[GridComponentUniqueGeometries4AssemblyType(AT)]
 
     @assert length(EGs) == 1 "MomentInterpolator currently only works with grids with a single element geometry"
     EG = EGs[1]
 
-    nitems::Int = num_sources(itemnodes)
     ncomponents::Int = get_ncomponents(FEType)
     edim::Int = dim_element(EG)
     order_FE = get_polynomialorder(FEType, EG)
@@ -403,6 +413,8 @@ function MomentInterpolator(
                 QP.xref = xref[qp]
                 eval_trafo!(QP.x, L2G, xref[qp])
 
+                QP.cell = itemcells[item]
+
                 exact_function!(result_f, QP)
                 if (bestapprox)
                     for m in 1:nmoments, k in 1:ncomponents
@@ -483,9 +495,8 @@ function MomentInterpolator(
         end
         QP.params = params === nothing ? [] : params
         QP.time = time
-        if isempty(items)
-            items = 1:nitems
-        end
+        isempty(items) && (items = 1:nitems)
+        isempty(itemcells) && (itemcells = 1:nitems)
         assembly_loop!(target, f_moments, items, exact_function!, QF, L2G, FEB, FEB_moments)
         return nothing
     end

src/lazy_interpolate.jl

@@ -87,43 +87,28 @@ function lazy_interpolate!(
     same_cells::Bool = xgrid == target.FES.xgrid
     CF::CellFinder{Tv, Ti} = CellFinder(xgrid)
 
-    if same_cells || use_cellparents == true
-        if same_cells
-            xCellParents = 1:num_cells(target.FES.xgrid)
-        else
-            xCellParents::Array{Ti, 1} = target.FES.xgrid[CellParents]
-        end
+    if same_cells || use_cellparents
+        xCellParents::Array{Ti, 1} = same_cells ? (1:num_cells(target.FES.xgrid)) : target.FES.xgrid[CellParents]
         function point_evaluation_parentgrid!(result, qpinfo)
-            x = qpinfo.x
-            cell = xCellParents[qpinfo.cell]
-            if xtrafo !== nothing
-                xtrafo(x_source, x)
-                cell = gFindLocal!(xref, CF, x_source; icellstart = cell, eps = eps, trybrute = !only_localsearch)
-            else
-                cell = gFindLocal!(xref, CF, x; icellstart = cell, eps = eps, trybrute = !only_localsearch)
-            end
-            evaluate_bary!(result, PE, xref, cell)
-            return nothing
+            x = xtrafo !== nothing ? xtrafo(x_source, qpinfo.x) : qpinfo.x
+            cell = gFindLocal!(xref, CF, x; icellstart = xCellParents[qpinfo.cell], eps = eps, trybrute = !only_localsearch)
+            return evaluate_bary!(result, PE, xref, cell)
         end
         fe_function = point_evaluation_parentgrid!
     else
         function point_evaluation_arbitrarygrids!(result, qpinfo)
-            x = qpinfo.x
-            if xtrafo !== nothing
-                xtrafo(x_source, x)
-                cell = gFindLocal!(xref, CF, x_source; icellstart = lastnonzerocell, eps = eps, trybrute = !only_localsearch)
-            else
-                cell = gFindLocal!(xref, CF, x; icellstart = lastnonzerocell, eps = eps, trybrute = !only_localsearch)
-            end
+            x = xtrafo !== nothing ? xtrafo(x_source, qpinfo.x) : qpinfo.x
+            cell = gFindLocal!(xref, CF, x; icellstart = lastnonzerocell, eps = eps, trybrute = !only_localsearch)
             if cell == 0
-                fill!(result, not_in_domain_value)
+                return fill!(result, not_in_domain_value)
             else
                 evaluate_bary!(result, PE, xref, cell)
                 lastnonzerocell = cell
+                return cell
             end
-            return nothing
         end
         fe_function = point_evaluation_arbitrarygrids!
     end
+
     return interpolate!(target, ON_CELLS, fe_function; resultdim = resultdim, items = items, kwargs...)
 end

test/test_interpolation_matrix.jl

@@ -214,7 +214,7 @@ function run_space_interpolation_matrix_tests()
         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"
+                @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