fixed hcurln1 interpolation and rearranged matrix interpolation code

Author: Da-Be-Ru

src/ExtendableFEMBase.jl

@@ -117,7 +117,6 @@ export interpolate! # must be defined separately by each FEdefinition
 export nodevalues, continuify
 export nodevalues!, nodevalues_subset!
 export nodevalues_view
-export compute_interpolation_jacobian
 
 export interpolator_matrix
 
@@ -175,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/hcurl_n1.jl

@@ -27,6 +27,8 @@ 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]
@@ -57,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...)
@@ -68,8 +70,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT},
     end
 
     # set values of interior N1 functions such that P0 moments are preserved
-    get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...)
-
+    return get_interpolator(FE, ON_CELLS).evaluate!(Target, data, items; kwargs...)
 end
 
 # on faces dofs are only tangential fluxes

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,63 +91,6 @@ function ExtendableGrids.interpolate!(target::FEVectorBlock, source; kwargs...)
     return interpolate!(target, ON_CELLS, source; kwargs...)
 end
 
-"""
-````
-function compute_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
-- This function can be used for computing prolongation or restriction operators if the `FESpace`s are defined on coarser/finer grids, respectively.
-
-"""
-function compute_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 nodevalues_subset!(

test/test_interpolation_matrix.jl

@@ -10,13 +10,6 @@ function run_grid_interpolation_matrix_tests()
         H1Pk{1, 1, 5} => 5,
     ]
 
-    PairTestCatalog1D = [
-        (L2P0{1}, H1P1{1}) => 0,
-        (H1P1{1}, H1P2{1, 1}) => 1,
-        (H1P2{1, 1}, H1P3{1, 1}) => 2,
-        (H1Pk{1, 1, 3}, H1Pk{1, 1, 5}) => 3,
-    ]
-
     TestCatalog2D = [
         HCURLN0{2} => 0,
         HCURLN1{2} => 1,
@@ -46,16 +39,6 @@ function run_grid_interpolation_matrix_tests()
         H1Pk{2, 2, 5} => 5,
     ]
 
-    PairTestCatalog2D = [
-        (H1P1{2}, HDIVRT0{2}) => 0,
-        (H1P2{2, 2}, HDIVRT0{2}) => 0,
-        (H1P2{2, 2}, HDIVRT1{2}) => 1,
-        (H1P2{2, 2}, HDIVRTk{2, 2}) => 2,
-        (HDIVRT1{2}, HDIVBDM1{2}) => 1,
-        (L2P0{2}, L2P1{2}) => 0,
-        (H1P2B{2, 2}, H1BR{2}) => 1,
-    ]
-
     TestCatalog3D = [
         HCURLN0{3} => 0,
         HDIVRT0{3} => 0,
@@ -72,15 +55,6 @@ function run_grid_interpolation_matrix_tests()
         H1P3{3, 3} => 3,
     ]
 
-    PairTestCatalog3D = [
-        (H1P1{3}, HDIVRT0{3}) => 0,
-        (H1P2{3, 3}, HDIVRT0{3}) => 0,
-        (H1P2{3, 3}, HDIVRT1{3}) => 1,
-        (HDIVRT1{3}, HDIVBDM1{3}) => 1,
-        (L2P0{3}, L2P1{3}) => 0,
-        (H1P2B{3, 3}, H1BR{3}) => 1,
-    ]
-
     # test interpolation of same space between refined grids
     function test_grid_matrix_computation(xgrid, FEType, order; broken::Bool = false, use_cellparents::Bool = false)
         u, ~ = exact_function(Val(dim_grid(xgrid)), order)
@@ -94,7 +68,7 @@ function run_grid_interpolation_matrix_tests()
         target_vector = FEVector(target_FES)
         interpolate!(target_vector[1], u; bonus_quadorder = order)
 
-        interpolation_matrix = compute_interpolation_jacobian(target_FES, source_FES; use_cellparents)
+        interpolation_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents)
         matrix_interpolated_entries = interpolation_matrix * source_vector.entries
 
         return @test norm(target_vector.entries - matrix_interpolated_entries) < tolerance
@@ -177,34 +151,6 @@ function run_space_interpolation_matrix_tests()
         (H1P2B{3, 3}, H1BR{3}) => 1,
     ]
 
-
-    # V1: H1-conforming space, computes matrix into RT0 space
-    function RT0_interpolator(V1::FESpace, VRT0::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
-        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
-
     # test interpolation for different elements on same grid
     function test_space_matrix_computation(xgrid, source_FEType, target_FEType, order; broken::Bool = false, use_cellparents::Bool = false)
         u, ~ = exact_function(Val(dim_grid(xgrid)), order)
@@ -218,7 +164,7 @@ function run_space_interpolation_matrix_tests()
         target_vector = FEVector(target_FES)
         interpolate!(target_vector[1], u; bonus_quadorder = order)
 
-        interpolation_matrix = compute_interpolation_jacobian(target_FES, source_FES; use_cellparents)
+        interpolation_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents)
         matrix_interpolated_entries = interpolation_matrix * source_vector.entries
 
         return @test norm(target_vector.entries - matrix_interpolated_entries) < tolerance
@@ -250,8 +196,8 @@ function run_space_interpolation_matrix_tests()
                     if (source_element, target_element) == (H1P1{2}, HDIVRT0{2}) && !broken
                         source_FES = FESpace{source_element}(xgrid; broken)
                         target_FES = FESpace{target_element}(xgrid; broken)
-                        autodiff_matrix = compute_interpolation_jacobian(target_FES, source_FES; use_cellparents = false)
-                        RT0_matrix = RT0_interpolator(source_FES, target_FES)
+                        autodiff_matrix = compute_lazy_interpolation_jacobian(target_FES, source_FES; use_cellparents = false)
+                        RT0_matrix = H1Pk_to_HDIVRT0_interpolator(target_FES, source_FES)
 
                         @test norm(autodiff_matrix' - RT0_matrix[1]) < tolerance
                     end

test/test_interpolators.jl

@@ -109,7 +109,7 @@ function run_interpolator_tests()
         show(devnull, Solution)
 
         # compute error
-        error = compute_error(Solution[1], u, order)
+        error = compute_error(Solution[1], u, order + 1)
         println("FEType = $FEType $(broken ? "broken" : "") $AT | ndofs = $(FES.ndofs) | order = $order | error = $(norm(error, Inf))")
         return @test norm(error) < tolerance
     end