From 836c351b7051c37a14971d3fcade941503c89df4 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Fri, 14 Nov 2025 11:26:11 +0100 Subject: [PATCH 1/4] extracted incident cell of items into qpinfo in interpolators for use_cellparents in lazy_interpolate! --- src/interpolators.jl | 24 ++++++++++++++++---- src/lazy_interpolate.jl | 37 ++++++++++--------------------- test/test_interpolation_matrix.jl | 2 +- 3 files changed, 33 insertions(+), 30 deletions(-) diff --git a/src/interpolators.jl b/src/interpolators.jl index d4a7c8c..d4d13d6 100644 --- a/src/interpolators.jl +++ b/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 diff --git a/src/lazy_interpolate.jl b/src/lazy_interpolate.jl index 9563e66..cb1c56d 100644 --- a/src/lazy_interpolate.jl +++ b/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 diff --git a/test/test_interpolation_matrix.jl b/test/test_interpolation_matrix.jl index d7775f3..79a7e97 100644 --- a/test/test_interpolation_matrix.jl +++ b/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 From 2017672569b8f3c4711c9eac1c12184aaeb475c8 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Fri, 21 Nov 2025 12:01:02 +0100 Subject: [PATCH 2/4] correct itemcell fill for various assembly type delegations in interpolations --- src/interpolators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolators.jl b/src/interpolators.jl index d4d13d6..0befcd0 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -501,7 +501,7 @@ function MomentInterpolator( QP.params = params === nothing ? [] : params QP.time = time isempty(items) && (items = 1:nitems) - isempty(itemcells) && (itemcells = items) + isempty(itemcells) && (itemcells = 1:nitems) assembly_loop!(target, f_moments, items, exact_function!, QF, L2G, FEB, FEB_moments) return nothing end From e77fffe736fda9285e62e87b13e4c2bbaa8d5dfb Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Fri, 21 Nov 2025 14:43:50 +0100 Subject: [PATCH 3/4] removed debug printouts and updated changelog --- CHANGELOG.md | 3 ++- src/interpolators.jl | 4 ---- src/lazy_interpolate.jl | 2 -- 3 files changed, 2 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 31497f9..0d3a5f4 100644 --- a/CHANGELOG.md +++ b/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}) diff --git a/src/interpolators.jl b/src/interpolators.jl index 0befcd0..15404c3 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -112,7 +112,6 @@ function NodalInterpolator( continue end cell = xNodeCells[1, j] - #@show cell QP.item = cell QP.cell = cell QP.region = xCellRegions[cell] @@ -187,8 +186,6 @@ function MomentInterpolator( kwargs... ) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} - #@info "MomentInterpolator $(AT)" - itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)] itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)] itemregions = xgrid[GridComponentRegions4AssemblyType(AT)] @@ -418,7 +415,6 @@ function MomentInterpolator( eval_trafo!(QP.x, L2G, xref[qp]) QP.cell = itemcells[item] - #@show QP.item QP.cell exact_function!(result_f, QP) if (bestapprox) diff --git a/src/lazy_interpolate.jl b/src/lazy_interpolate.jl index cb1c56d..24b22fc 100644 --- a/src/lazy_interpolate.jl +++ b/src/lazy_interpolate.jl @@ -89,10 +89,8 @@ function lazy_interpolate!( 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 = 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 From b7b9645ef3337b5a281b26e477b6c8ba9708b693 Mon Sep 17 00:00:00 2001 From: Daniel Runge Date: Sat, 22 Nov 2025 14:28:31 +0100 Subject: [PATCH 4/4] adjusted itemcells computation for ON_EDGES assembly in MomentInterpolator --- src/interpolators.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/interpolators.jl b/src/interpolators.jl index 15404c3..4c5174f 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -202,11 +202,10 @@ function MomentInterpolator( 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 + # 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