FIMH2021 verification

docs/src/literate-tutorials/ep04_geselowitz-ecg.jl

@@ -118,7 +118,7 @@ integrator = init(problem, timestepper, dt=dt₀, verbose=true)
 
 # Now that the time integrator is ready we setup the ECG problem.
 torso_mesh_κᵢ = ConstantCoefficient(1.0)
-torso_mesh_κ  = ConstantCoefficient(1.0)
+torso_mesh_κ  = ConstantCoefficient(2.0)
 # !!! todo
 #     Show how to transfer `diffusion_tensor_field` onto the torso mesh.
 geselowitz_ecg = Thunderbolt.Geselowitz1989ECGLeadCache(

src/ferrite-addons/transfer_operators.jl

@@ -80,7 +80,7 @@ struct NodalIntergridInterpolation{PH <: PointEvalHandler, DH1 <: AbstractDofHan
         grid_to   = Ferrite.get_grid(dh_to)
         grid_from = Ferrite.get_grid(dh_from)
         nodes = Vector{Ferrite.get_coordinate_type(grid_to)}(undef, length(dofset))
-        for sdh in dh_to.subdofhandlers[subdomains_from]
+        for sdh in dh_to.subdofhandlers[subdomains_to]
             # Skip subdofhandler if field is not present
             field_name_to ∈ Ferrite.getfieldnames(sdh) || continue
             # Grab the reference coordinates of the field to interpolate
@@ -94,7 +94,7 @@ struct NodalIntergridInterpolation{PH <: PointEvalHandler, DH1 <: AbstractDofHan
             _compute_dof_nodes_barrier!(nodes, sdh, Ferrite.dof_range(sdh, field_name_to), gip, dof_to_node_map, ref_coords)
         end
 
-        ph = PointEvalHandler(Ferrite.get_grid(dh_from), nodes; warn=false)
+        ph = PointEvalHandler(Ferrite.get_grid(dh_from), nodes, [sdh.cellset for sdh in dh_from.subdofhandlers[subdomains_from]]; warn=false)
 
         n_missing = sum(x -> x === nothing, ph.cells)
         n_missing == 0 || @warn "Constructing the interpolation for $field_name_from to $field_name_to failed. $n_missing (out of $(length(ph.cells))) points not found."

src/modeling/electrophysiology/ecg.jl

@@ -1,29 +1,30 @@
+# TODO we should use this function to study how we can pull the fluxes with the operator infrastructure directly
 function compute_quadrature_fluxes!(fluxdata, dh, u, field_name, integrator)
     grid = get_grid(dh)
     sdim = getspatialdim(grid)
     for sdh in dh.subdofhandlers
-        ip          = Ferrite.getfieldinterpolation(sdh, field_name)
-        firstcell   = getcells(grid, first(sdh.cellset))
-        ip_geo      = Ferrite.geometric_interpolation(typeof(firstcell))^sdim
-        element_qr  = getquadraturerule(integrator.qrc, firstcell)
-        cv = CellValues(element_qr, ip, ip_geo)
-        _compute_quadrature_fluxes_on_subdomain!(fluxdata,sdh,cv,u,integrator)
+        element_cache  = setup_element_cache(integrator, sdh)
+        # Function barrier
+        _compute_quadrature_fluxes_on_subdomain!(fluxdata,sdh,u,element_cache)
     end
 end
 
-function _compute_quadrature_fluxes_on_subdomain!(κ∇u,sdh,cv,u,integrator::BilinearDiffusionIntegrator)
-    n_basefuncs = getnbasefunctions(cv)
+function _compute_quadrature_fluxes_on_subdomain!(κ∇u,sdh,u,element_cache::BilinearDiffusionElementCache)
+    @unpack cellvalues, Dcache = element_cache
+    n_basefuncs = getnbasefunctions(cellvalues)
+
+    uₑ = zeros(n_basefuncs)
     for cell ∈ CellIterator(sdh)
         κ∇ucell = get_data_for_index(κ∇u, cellid(cell))
 
-        reinit!(cv, cell)
-        uₑ = @view u[celldofs(cell)]
+        reinit!(cellvalues, cell)
+        uₑ .= u[celldofs(cell)]
 
-        for qp in QuadratureIterator(cv)
-            D_loc = evaluate_coefficient(integrator.D, cell, qp, time)
+        for qp in QuadratureIterator(cellvalues)
+            D_loc = evaluate_coefficient(Dcache, cell, qp, time)
             # dΩ = getdetJdV(cellvalues, qp)
             for i in 1:n_basefuncs
-                ∇Nᵢ = shape_gradient(cv, qp, i)
+                ∇Nᵢ = shape_gradient(cellvalues, qp, i)
                 κ∇ucell[qp.i] += D_loc ⋅ ∇Nᵢ * uₑ[i]
             end
         end
@@ -52,12 +53,12 @@ struct Plonsey1964ECGGaussCache{BufferType, OperatorType}
     op::OperatorType
 end
 
-function Plonsey1964ECGGaussCache(op::AssembledBilinearOperator, φₘ::AbstractVector{T}) where T
+function Plonsey1964ECGGaussCache(op::AssembledBilinearOperator, φₘ::AbstractVector{T}, subdomains::AbstractVector{String} = [""]) where T
     @unpack dh, integrator = op
     @assert length(dh.field_names) == 1 "Multiple fields detected. Problem setup might be broken..."
     grid = get_grid(dh)
     sdim = Ferrite.getspatialdim(grid)
-    κ∇φₘ = construct_qvector(Vector{Vec{sdim,T}}, Vector{Int64}, grid, integrator.qrc)
+    κ∇φₘ = construct_qvector(Vector{Vec{sdim,T}}, Vector{Int64}, grid, integrator.qrc, subdomains)
     compute_quadrature_fluxes!(κ∇φₘ, dh, φₘ, dh.field_names[1], integrator)
     Plonsey1964ECGGaussCache(κ∇φₘ, op)
 end

test/integration/test_ecg.jl

@@ -1,10 +1,11 @@
 using Thunderbolt, Test, SparseArrays
 import Thunderbolt: to_mesh, OrderedSet
+using FerriteGmsh
 
 @testset "ECG" begin
     @testset "Blocks with $geo" for geo in (Tetrahedron,Hexahedron)
         size = 2.
-        signal_strength = 0.04
+        signal_strength = 0.02
         nel_heart = (6, 6, 6) .* 1
         nel_torso = (8, 8, 8) .* Int(size)
         ground_vertex = Vec(0.0, 0.0, 0.0)
@@ -48,7 +49,7 @@ import Thunderbolt: to_mesh, OrderedSet
                 QuadratureRuleCollection(2),
                 :φₘ,
             ),
-            Thunderbolt.SparseMatrixCSC,
+            SparseMatrixCSC{Float64,Int64},
             heart_fun.dh,
         )
         Thunderbolt.update_operator!(op, 0.0) # trigger assembly
@@ -274,4 +275,107 @@ import Thunderbolt: to_mesh, OrderedSet
             end
         end
     end
+
+    @testset "FIMH2021" begin
+        heart_elsize = 0.1
+        torso_elsize = 5.0
+        heart_radius = 1.0
+        torso_radius = 50.0
+        λ = 0.2
+
+        gmsh.initialize()
+        gmsh.model.add("sis-50")
+        gmsh.model.occ.add_disk(0.0, 0.0, 0.0, torso_radius, torso_radius, 1)
+        heart_handle = gmsh.model.occ.add_disk(0.0, 0.0, 0.0, heart_radius, heart_radius, 2)
+        gmsh.model.occ.cut([(2,1)], [(2,2)], 3, true, false)
+        handle2 = gmsh.model.occ.fragment([(2,3)], [(2,2)])
+        gmsh.model.occ.synchronize()
+        gmsh.model.mesh.renumberNodes()
+        gmsh.model.mesh.renumberElements()
+        # Set size of heart and torso
+        gmsh.model.mesh.setSize([(0,2)], heart_elsize);
+        gmsh.model.mesh.setSize([(0,3)], torso_elsize);
+        gmsh.model.mesh.generate(2)
+        gmsh.model.addPhysicalGroup(2, [2], -1, "heart")
+        gmsh.model.addPhysicalGroup(2, [3], -1, "torso")
+        nodes = tonodes()
+        elements, gmsh_elementidx = toelements(2)
+        cellsets = tocellsets(2, gmsh_elementidx)
+        gmsh.finalize()
+        grid = Grid(elements, nodes, cellsets=Dict(["heart" => cellsets["heart"], "torso" => cellsets["torso"]]))
+        torso_grid = heart_grid = to_mesh(grid)
+
+        κ  = AnalyticalCoefficient(
+            (x,t) -> norm(x,2) ≤ heart_radius ? SymmetricTensor{2,2,Float64,3}((2λ, 0, 2λ)) : SymmetricTensor{2,2,Float64,3}((λ, 0.0, λ)),
+            CartesianCoordinateSystem{2}()
+        )
+        κᵢ = AnalyticalCoefficient(
+            (x,t) -> norm(x,2) ≤ heart_radius ? SymmetricTensor{2,2,Float64,3}((λ, 0, λ)) : SymmetricTensor{2,2,Float64,3}((0.0, 0.0, 0.0)),
+            CartesianCoordinateSystem{2}()
+        )
+
+        heart_model = TransientDiffusionModel(
+            κᵢ,
+            NoStimulationProtocol(), # Poisoning to detecte if we accidentally touch these
+            :φₘ
+        )
+        heart_fun = semidiscretize(
+            heart_model,
+            FiniteElementDiscretization(
+                Dict(:φₘ => LagrangeCollection{1}()),
+                Dirichlet[],
+                ["heart"]
+            ),
+            heart_grid
+        )
+
+        ground_vertex = Vec(torso_radius, 0.0)
+        electrodes = [
+            ground_vertex,
+            Vec(0., 10.),
+        ]
+        electrode_pairs = [[2,1]]
+
+        op = Thunderbolt.setup_assembled_operator(
+            Thunderbolt.BilinearDiffusionIntegrator(
+                κ,
+                QuadratureRuleCollection(2),
+                :φₘ,
+            ),
+            SparseMatrixCSC{Float64,Int64},
+            heart_fun.dh,
+        )
+        Thunderbolt.update_operator!(op, 0.0) # trigger assembly
+
+        u = zeros(Thunderbolt.solution_size(heart_fun))
+        apply_analytical!(u, heart_fun.dh, :φₘ, x->max(0.0,norm(x-Vec((0.0,-1.0)))), getcellset(heart_grid, "heart"))
+
+        plonsey_ecg = Thunderbolt.Plonsey1964ECGGaussCache(op, u, ["heart"])
+        Thunderbolt.update_ecg!(plonsey_ecg, u)
+        plonsey_vals = Thunderbolt.evaluate_ecg(plonsey_ecg, electrodes[2], λ)
+
+        poisson_ecg = Thunderbolt.PoissonECGReconstructionCache(
+            heart_fun,
+            torso_grid,
+            κᵢ, κ,
+            electrodes;
+            ground               = OrderedSet([Thunderbolt.get_closest_vertex(ground_vertex, torso_grid)]),
+            torso_heart_domain   = ["heart"],
+            ipc                  = LagrangeCollection{1}(),
+            qrc                  = QuadratureRuleCollection(2),
+            linear_solver        = Thunderbolt.LinearSolve.UMFPACKFactorization(),
+            system_matrix_type   = SparseMatrixCSC{Float64,Int64}
+        )
+        Thunderbolt.update_ecg!(poisson_ecg, u)
+        poisson_vals = Thunderbolt.evaluate_ecg(poisson_ecg)
+        # VTKGridFile("FIMH2021-Debug", grid) do vtk
+        #     Ferrite.write_cellset(vtk, grid)
+        #     Ferrite.write_solution(vtk, poisson_ecg.torso_op.dh, poisson_ecg.φₘ_t, "1")
+        #     Ferrite.write_solution(vtk, poisson_ecg.torso_op.dh, poisson_ecg.κ∇φₘ_t, "2")
+        #     Ferrite.write_solution(vtk, poisson_ecg.torso_op.dh, poisson_ecg.ϕₑ, "3")
+        # end
+
+        # TODO investigate where the 2π comes from
+        @test plonsey_vals*2π ≈ poisson_vals[1] - poisson_vals[2] rtol = 1e-1
+    end
 end

test/test_transfer.jl

@@ -64,10 +64,14 @@
         end
     end
 
-    target_grid_nonmatching = generate_grid(Triangle, (40, 44), Vec((-2.0,-2.0)), Vec((2.0,2.0)))
+    # target_grid_nonmatching = generate_grid(Triangle, (40, 44), Vec((-2.0,-2.0)), Vec((2.0,2.0)))
+    target_grid_nonmatching = generate_grid(Triangle, (4, 4), Vec((-2.0,-2.0)), Vec((2.0,2.0)))
     addcellset!(target_grid_nonmatching, "hole", x->norm(x) ≤ 1.0)
-    addcellset!(target_grid_nonmatching, "remaining", x->norm(x) ≥ 1.0)
+    target_grid_nonmatching.cellsets["remaining"] = OrderedSet([i for i in 1:getncells(target_grid_nonmatching) if i ∉ target_grid_nonmatching.cellsets["hole"]])
     target_mesh_nonmatching = to_mesh(target_grid_nonmatching)
+    # VTKGridFile("FIMH2021-Debug", target_grid_nonmatching) do vtk
+    #     Ferrite.write_cellset(vtk, target_grid_nonmatching)
+    # end
 
     @testset "Nodal Intergrid Interpolation on Non-Matching Grids" begin
         source_dh = DofHandler(source_mesh)
@@ -124,4 +128,67 @@
             end
         end
     end
+
+    @testset "Nodal Intergrid Interpolation between Subdomains" failfast=true begin
+        source_dh = DofHandler(target_mesh_nonmatching)
+        source_sdh_hole = SubDofHandler(source_dh, getcellset(target_mesh_nonmatching, "hole"))
+        add!(source_sdh_hole, :v, Lagrange{RefTriangle,2}())
+        add!(source_sdh_hole, :z, Lagrange{RefTriangle,1}())
+        # source_sdh_remaining = SubDofHandler(source_dh, getcellset(target_mesh_nonmatching, "remaining"))
+        # add!(source_sdh_remaining, :v, Lagrange{RefTriangle,2}())
+        # add!(source_sdh_remaining, :w, Lagrange{RefTriangle,1}())
+        close!(source_dh)
+
+        source_u = ones(ndofs(source_dh))
+        apply_analytical!(source_u, source_dh, :v, x->-norm(x))
+        apply_analytical!(source_u, source_dh, :z, x-> norm(x))
+
+        target_dh = DofHandler(target_mesh_nonmatching)
+        target_sdh_hole = SubDofHandler(target_dh, getcellset(target_mesh_nonmatching, "hole"))
+        add!(target_sdh_hole, :v, Lagrange{RefTriangle,2}())
+        add!(target_sdh_hole, :w, Lagrange{RefTriangle,1}())
+        # target_sdh_remaining = SubDofHandler(target_dh, getcellset(target_mesh_nonmatching, "remaining"))
+        # add!(target_sdh_remaining, :v, Lagrange{RefTriangle,2}())
+        # add!(target_sdh_remaining, :w, Lagrange{RefTriangle,1}())
+        close!(target_dh)
+
+        v_range = dof_range(target_dh.subdofhandlers[1], :v)
+        w_range = dof_range(target_dh.subdofhandlers[1], :w)
+
+        target_sdhids = Thunderbolt.get_subdofhandler_indices_on_subdomains(target_dh, ["hole"])
+
+        op = Thunderbolt.NodalIntergridInterpolation(source_dh, target_dh, :v, :v; subdomains_to=target_sdhids)
+        target_u = [NaN for i in 1:ndofs(target_dh)]
+        Thunderbolt.transfer!(target_u, op, source_u)
+        cvv = CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefTriangle,2}())
+        for cc in CellIterator(target_dh.subdofhandlers[1])
+            reinit!(cvv, cc)
+            dofs_v = @view celldofs(cc)[v_range]
+            dofs_w = @view celldofs(cc)[w_range]
+            for qp in QuadratureIterator(cvv)
+                x = Thunderbolt.spatial_coordinate(Lagrange{RefTriangle,1}(), qp.ξ, getcoordinates(cc))
+                @test function_value(cvv, qp, target_u[dofs_v]) ≈ -norm(x) atol=3e-1
+                @test all(isnan.(target_u[dofs_w]))
+            end
+        end
+        VTKGridFile("FIMH2021-Debug", target_grid_nonmatching) do vtk
+            Ferrite.write_cellset(vtk, target_grid_nonmatching)
+            Ferrite.write_solution(vtk, target_dh, target_u, "tgt")
+            Ferrite.write_solution(vtk, source_dh, source_u, "src")
+        end
+        op = Thunderbolt.NodalIntergridInterpolation(source_dh, target_dh, :z, :w; subdomains_to=target_sdhids)
+        target_u = [NaN for i in 1:ndofs(target_dh)]
+        Thunderbolt.transfer!(target_u, op, source_u)
+        cvw = CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefTriangle,1}())
+        for cc in CellIterator(target_dh.subdofhandlers[1])
+            reinit!(cvw, cc)
+            dofs_v = @view celldofs(cc)[v_range]
+            dofs_w = @view celldofs(cc)[w_range]
+            for qp in QuadratureIterator(cvw)
+                x = Thunderbolt.spatial_coordinate(Lagrange{RefTriangle,1}(), qp.ξ, getcoordinates(cc))
+                @test all(isnan.(target_u[dofs_v]))
+                @test function_value(cvw, qp, target_u[dofs_w]) ≈ norm(x) atol=3e-1
+            end
+        end
+    end
 end