added function for computing jacobian of interpolation between finite element spaces

Author: Da-Be-Ru

Project.toml

@@ -1,10 +1,11 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 version = "1.4.0"
+authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
@@ -13,6 +14,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 SpecialPolynomials = "a25cea48-d430-424a-8ee7-0d3ad3742e9e"
 
 [weakdeps]
@@ -23,6 +26,7 @@ ExtendableFEMBaseUnicodePlotsExt = ["UnicodePlots"]
 
 [compat]
 DiffResults = "1"
+DifferentiationInterface = "0.7.10"
 DocStringExtensions = "0.8,0.9"
 ExtendableGrids = "1.13.0"
 ExtendableSparse = "1.5.1"
@@ -31,6 +35,8 @@ LinearAlgebra = "1.9"
 Polynomials = "2.0.21, 3, 4"
 Printf = "1.9"
 SparseArrays = "1.9"
+SparseConnectivityTracer = "1.1.2"
+SparseMatrixColorings = "0.4.22"
 SpecialPolynomials = "0.4.9, 0.5"
 UnicodePlots = "3.6"
 julia = "1.9"

src/ExtendableFEMBase.jl

@@ -46,13 +46,16 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     AbstractExtendableSparseMatrixCSC, ExtendableSparseMatrixCSC, MTExtendableSparseMatrixCSC,
     rawupdateindex!
+using DifferentiationInterface: AutoForwardDiff, AutoSparse, jacobian
 using ForwardDiff: ForwardDiff, DiffResults
 using LinearAlgebra: LinearAlgebra, convert, det, diagm, dot, eigen, ldiv!, lu,
     mul!, norm, transpose
 using Polynomials: Polynomials, Polynomial, coeffs
 using Printf: Printf, @printf
 using SparseArrays: SparseArrays, AbstractSparseArray, AbstractSparseMatrix,
     SparseMatrixCSC, nzrange, rowvals
+using SparseConnectivityTracer: TracerSparsityDetector
+using SparseMatrixColorings: GreedyColoringAlgorithm
 using SpecialPolynomials: SpecialPolynomials, ShiftedLegendre, basis
 
 include("functionoperators.jl")
@@ -114,6 +117,7 @@ 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
 

src/interpolations.jl

@@ -91,6 +91,62 @@ 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
+
 
 """
 ````

src/interpolators.jl

@@ -65,6 +65,9 @@ function NodalInterpolator(
         xCellDofs = FES[CellDofs]
         ncells = num_cells(grid)
         function evaluate_broken!(target, exact_function!, items; time = 0, params = [], kwargs...)
+            if !(eltype(target) <: T)
+                result = zeros(eltype(target), ncomponents)
+            end
             QP.time = time
             QP.params = params === nothing ? [] : params
             if isempty(items)
@@ -96,6 +99,9 @@ function NodalInterpolator(
         nnodes = num_nodes(grid)
         xNodeCells = atranspose(grid[CellNodes])
         function evaluate!(target, exact_function!, items; time = 0, params = [], kwargs...)
+            if !(eltype(target) <: T)
+                result = zeros(eltype(target), ncomponents)
+            end
             QP.time = time
             QP.params = params === nothing ? [] : params
             if isempty(items)
@@ -370,6 +376,10 @@ function MomentInterpolator(
     interiordofs = zeros(Int, length(idofs))
 
     function assembly_loop!(target, f_moments, items, exact_function!, QF, L2G, FEB, FEB_moments)
+        if !(eltype(target) <: Tv)
+            result_f = zeros(eltype(target), ncomponents)
+            f_moments = zeros(eltype(target), nmoments)
+        end
         weights, xref = QF.w, QF.xref
         nweights = length(weights)
         for item::Int in items
@@ -582,6 +592,10 @@ function FunctionalInterpolator(
     FEB = FEEvaluator(FE, operator, QF; AT = AT, T = Tv, L2G = L2G)
 
     function assembly_loop!(target, f_fluxes, items, exact_function!, QF, L2G, FEB)
+        if !(eltype(target) <: Tv)
+            result_f = zeros(eltype(target), ncomponents)
+            f_fluxes = zeros(eltype(target), nfluxes)
+        end
         weights, xref = QF.w, QF.xref
         nweights = length(weights)
         for item::Int in items

src/lazy_interpolate.jl

@@ -79,7 +79,7 @@ function lazy_interpolate!(
     if xdim_source != xdim_target
         @assert xtrafo !== nothing "grids have different coordinate dimensions, need xtrafo!"
     end
-    PE = PointEvaluator(postprocess, operators, source)
+    PE = PointEvaluator(postprocess, operators, source; TCoeff = T1)
     xref = zeros(Tv, xdim_source)
     x_source = zeros(Tv, xdim_source)
     cell::Int = start_cell

src/point_evaluator.jl

@@ -1,4 +1,4 @@
-mutable struct PointEvaluator{Tv <: Real, UT, KFT <: Function}
+mutable struct PointEvaluator{Tv <: Real, TCoeff <: Real, UT, KFT <: Function}
     u_args::Array{UT, 1}
     ops_args::Array{DataType, 1}
     kernel::KFT
@@ -61,11 +61,11 @@ $(_myprint(default_peval_kwargs()))
 After construction, call `initialize!` to prepare the evaluator for a given solution, then use `evaluate!` or `evaluate_bary!` to perform point evaluations.
 
 """
-function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, kwargs...)
+function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, TCoeff = Float64, kwargs...)
     parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_peval_kwargs())
     _update_params!(parameters, kwargs)
     @assert length(u_args) == length(ops_args)
-    PE = PointEvaluator{Tv, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, 1, nothing, nothing, nothing, zeros(Tv, 2), parameters)
+    PE = PointEvaluator{Tv, TCoeff, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, 1, nothing, nothing, nothing, zeros(Tv, 2), parameters)
     if sol !== nothing
         initialize!(PE, sol)
     end
@@ -121,7 +121,7 @@ $(_myprint(default_peval_kwargs()))
 # Notes
 - This function must be called before using `evaluate!` or `evaluate_bary!` with the `PointEvaluator`.
 """
-function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {T, UT}
+function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...) where {T, TCoeff, UT}
     _update_params!(O.parameters, kwargs)
     if UT <: Integer
         ind_args = O.u_args
@@ -159,7 +159,7 @@ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {
     op_lengths_args = [size(O.BE_args[1][j].cvals, 1) for j in 1:nargs]
     op_offsets_args = [0]
     append!(op_offsets_args, cumsum(op_lengths_args))
-    input_args = zeros(T, op_offsets_args[end])
+    input_args = zeros(TCoeff, op_offsets_args[end])
 
     FEATs_args = [EffAT4AssemblyType(get_AT(FES_args[j]), AT) for j in 1:nargs]
     itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j in 1:nargs]

test/Project.toml

@@ -1,13 +1,12 @@
 [deps]
-ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
+ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]

test/runtests.jl

@@ -1,6 +1,7 @@
 using Test
 using ExtendableGrids
 using ExtendableFEMBase
+using ExtendableSparse
 using ExplicitImports
 using ExampleJuggler
 using SparseArrays
@@ -28,6 +29,7 @@ end
 
 include("test_quadrature.jl")
 include("test_interpolators.jl")
+include("test_interpolation_matrix.jl")
 include("test_operators.jl")
 include("test_febasis.jl")
 include("test_segmentintegrator.jl")
@@ -150,6 +152,8 @@ function run_all_tests()
         run_operator_tests()
         run_quadrature_tests()
         run_interpolator_tests()
+        run_grid_interpolation_matrix_tests()
+        run_space_interpolation_matrix_tests()
         run_segmentintegrator_tests()
         run_pointevaluator_tests()
     end