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

Author: Da-Be-Ru

src/interpolators.jl

@@ -112,6 +112,7 @@ function NodalInterpolator(
                     continue
                 end
                 cell = xNodeCells[1, j]
+                #@show cell
                 QP.item = cell
                 QP.cell = cell
                 QP.region = xCellRegions[cell]
@@ -186,24 +187,37 @@ function MomentInterpolator(
         kwargs...
     ) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
 
+    #@info "MomentInterpolator $(AT)"
+
     itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)]
     itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)]
     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]
+            itemcells = []
+            for item in 1:nitems
+                push!(itemcells, ec[1, item])
+            end
+        else
+            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 +417,9 @@ function MomentInterpolator(
                 QP.xref = xref[qp]
                 eval_trafo!(QP.x, L2G, xref[qp])
 
+                QP.cell = itemcells[item]
+                #@show QP.item QP.cell
+
                 exact_function!(result_f, QP)
                 if (bestapprox)
                     for m in 1:nmoments, k in 1:ncomponents
@@ -483,9 +500,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 = items)
         assembly_loop!(target, f_moments, items, exact_function!, QF, L2G, FEB, FEB_moments)
         return nothing
     end

src/lazy_interpolate.jl

@@ -87,43 +87,30 @@ 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]
+        #@show xCellParents
         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
+            #@show qpinfo.x qpinfo.cell
+            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