Fix remaining issues

Author: termi-official

src/mesh/meshes.jl

@@ -40,6 +40,7 @@ function add_subdomain!(dh::DofHandler{<:Any, <:SimpleMesh}, name::String, appro
     cells = mesh.grid.cells
     haskey(mesh.volumetric_subdomains, name) || error("Volumetric Subdomain $name not found on mesh. Available subdomains: $(keys(mesh.volumetric_subdomains))")
     for (celltype, cellset) in mesh.volumetric_subdomains[name].data
+        # @info name, length(cellset)
         sdh = SubDofHandler(dh, OrderedSet{Int}([idx.idx for idx in cellset]))
         for ad in approxmations
             add!(sdh, ad.sym, getinterpolation(ad.ipc, cells[first(sdh.cellset)]))

src/modeling/electrophysiology/ecg.jl

@@ -154,10 +154,10 @@ stiffness matrix with applied homogeneous Dirichlet boundary condition at the fi
 """
 struct PoissonECGReconstructionCache{DiffusionOperatorType1, DiffusionOperatorType2, TransferOperatorType, SolutionVectorType, SolverCacheType, PHType, CHType}
     torso_op::DiffusionOperatorType1  # Operator on the torso mesh for ∇κ∇
-    heart_op::DiffusionOperatorType2  # Operator on the heart mesh for ∇κᵢ∇
+    source_op::DiffusionOperatorType2  # Operator on the heart mesh for ∇κᵢ∇
     transfer_op::TransferOperatorType # Transfer from heart to torso mesh
     ϕₑ::SolutionVectorType            # Solution vector buffer
-    κ∇φₘ_h::SolutionVectorType        # Source term buffer on heart
+    φₘ_t::SolutionVectorType          # Solution vector buffer on torso
     κ∇φₘ_t::SolutionVectorType        # Source term buffer on torso
     inner_solver::SolverCacheType     # Linear solver
     ph::PHType                        # PointEvalHandler on the torso
@@ -237,16 +237,16 @@ function PoissonECGReconstructionCache(
         subdomains_to = get_subdofhandler_indices_on_subdomains(torso_dh, torso_heart_domain)
     )
 
-    heart_op = setup_assembled_operator(
+    source_op = setup_assembled_operator(
         BilinearDiffusionIntegrator(
             heart_diffusion_tensor_field,
             qrc,
             extracellular_potential_symbol,
         ),
         system_matrix_type,
-        heart_dh,
+        torso_dh,
     )
-    update_operator!(heart_op, 0.) # Trigger assembly
+    update_operator!(source_op, 0.) # Trigger assembly
 
     torso_op = setup_assembled_operator(
         BilinearDiffusionIntegrator(
@@ -268,7 +268,7 @@ function PoissonECGReconstructionCache(
     PoissonECGReconstructionCache(
         heart_fun,
         torso_fun,
-        heart_op,
+        source_op,
         torso_op,
         transfer_op,
         ph;
@@ -280,7 +280,7 @@ end
 function PoissonECGReconstructionCache(
     heart_fun::AffineODEFunction,
     torso_fun::AffineSteadyStateFunction,
-    heart_op::AssembledBilinearOperator,
+    source_op::AssembledBilinearOperator,
     torso_op::AssembledBilinearOperator,
     transfer_op::AbstractTransferOperator,
     ph::PointEvalHandler;
@@ -292,7 +292,7 @@ function PoissonECGReconstructionCache(
     grid = get_grid(torso_dh)
     length(torso_dh.field_names) == 1 || @warn "Multiple fields detected. Setup might be broken..."
 
-    κ∇φₘh = create_system_vector(solution_vector_type, heart_fun) # RHS buffer before transfer
+    φₘt   = create_system_vector(solution_vector_type, torso_fun) # RHS buffer for source term
     κ∇φₘt = create_system_vector(solution_vector_type, torso_fun) # RHS buffer after transfer
     κ∇φₘt .= 0.0
     ϕₑ    = create_system_vector(solution_vector_type, torso_fun) # Solution vector
@@ -302,14 +302,14 @@ function PoissonECGReconstructionCache(
     )
     lincache = init(linprob, linear_solver)
 
-    return PoissonECGReconstructionCache(torso_op, heart_op, transfer_op, ϕₑ, κ∇φₘh, κ∇φₘt, lincache, ph, torso_ch)
+    return PoissonECGReconstructionCache(torso_op, source_op, transfer_op, ϕₑ, φₘt, κ∇φₘt, lincache, ph, torso_ch)
 end
 
 function update_ecg!(cache::PoissonECGReconstructionCache, φₘ::AbstractVector)
-    # Compute κᵢ∇φₘ on the heart
-    mul!(cache.κ∇φₘ_h, cache.heart_op, φₘ)
-    # Transfer κᵢ∇φₘ to the torso
-    transfer!(cache.κ∇φₘ_t, cache.transfer_op, cache.κ∇φₘ_h)
+    # Transfer φₘ to the torso
+    transfer!(cache.φₘ_t, cache.transfer_op, φₘ)
+    # Compute κᵢ∇φₘ on the torso
+    mul!(cache.κ∇φₘ_t, cache.source_op, cache.φₘ_t)
     cache.κ∇φₘ_t[isnan.(cache.κ∇φₘ_t)] .= 0.0 # FIXME
     # "Move to right hand side
     cache.κ∇φₘ_t .*= -1.0
@@ -348,9 +348,9 @@ V(t)=\\int \\nabla Z(\\boldsymbol{x}) \\cdot \\boldsymbol{\\kappa}_\\mathrm{i} \
 ```
 """
 struct Geselowitz1989ECGLeadCache{TZ <: AbstractMatrix, DiffusionOperatorType, TransferOperatorType, SolutionVectorType <: AbstractVector, ElectrodesVecType}
-    heart_op::DiffusionOperatorType   # Operator on the heart mesh for ∇κᵢ∇
+    source_op::DiffusionOperatorType   # Operator on the heart mesh for ∇κᵢ∇
     transfer_op::TransferOperatorType # Transfer from heart to torso mesh
-    κ∇φₘ_h::SolutionVectorType        # Source term buffer on heart
+    φₘ_t::SolutionVectorType          # Potential field on torso
     κ∇φₘ_t::SolutionVectorType        # Source term buffer on torso
     Z::TZ                             # Lead field
     electrode_positions::ElectrodesVecType
@@ -419,7 +419,7 @@ end
 function Geselowitz1989ECGLeadCache(
     heart_fun::AffineODEFunction,
     torso_grid::AbstractGrid,
-    heart_diffusion_tensor_field, # κᵢ - diffusion tensor description for heart on heart grid
+    heart_diffusion_tensor_field, # κᵢ - diffusion tensor description for heart on heart grid alone
     full_diffusion_tensor_field,  # κ - diffusion tensor description for heart and torso on torso grid
     electrode_positions::AbstractVector{Vector{VertexIndex}};
     ipc                  = LagrangeCollection{1}(),
@@ -454,15 +454,14 @@ function Geselowitz1989ECGLeadCache(
         torso_grid
     )
 
-    heart_grid = get_grid(heart_fun.dh)
     sourcefun = semidiscretize(
         source_model,
         FiniteElementDiscretization(
             Dict(tmpsym => ipc),
             Dirichlet[],
-            subdomain_names(heart_grid),
+            subdomain_names(torso_grid),
         ),
-        heart_grid,
+        torso_grid,
     )
 
     ϕₘ_op = setup_assembled_operator(
@@ -472,7 +471,7 @@ function Geselowitz1989ECGLeadCache(
             tmpsym,
         ),
         system_matrix_type,
-        sourcefun.dh,
+        lead_field_fun.dh,
     )
     update_operator!(ϕₘ_op, 0.) # Trigger assembly
 
@@ -525,14 +524,13 @@ function Geselowitz1989ECGLeadCache(
     lead_dh  = lead_op.dh
     length(lead_dh.field_names) == 1 || @warn "Multiple fields detected. Setup might be broken..."
     nelectrodes = length(electrode_positions)
-    ∇Njκ∇φₘ_h   = create_system_vector(solution_vector_type, heart_fun) # Solution vector
-    ∇Njκ∇φₘ_t   = create_system_vector(solution_vector_type, lead_fun)  # Solution vector
-    Z           = zeros(eltype(∇Njκ∇φₘ_t), nelectrodes, length(∇Njκ∇φₘ_t))
+    φₘ_t        = create_system_vector(solution_vector_type, lead_fun) # Solution vector
+    ∇φₘ_t       = create_system_vector(solution_vector_type, lead_fun)  # Solution vector
+    Z           = zeros(eltype(∇φₘ_t), nelectrodes, length(∇φₘ_t))
     ϕₑ          = zeros(nelectrodes)
 
-    lead_rhs = zeros(eltype(∇Njκ∇φₘ_t), nelectrodes, length(∇Njκ∇φₘ_t))
+    lead_rhs = zeros(eltype(∇φₘ_t), nelectrodes, length(∇φₘ_t))
 
-    @info size(lead_op.A), size(lead_rhs[1,:])
     leadprob = LinearSolve.LinearProblem(
         lead_op.A, copy(lead_rhs[1,:])
     )
@@ -544,17 +542,12 @@ function Geselowitz1989ECGLeadCache(
         for j in 2:length(electrode_set)
             _add_electrode!(current_rhs, lead_dh, electrode_set[j], -1.0/(length(electrode_set)-1), lead_field_sym)
         end
-        # leadprob = LinearSolve.LinearProblem(
-        #     lead_op.A, copy(current_rhs);# u0=Z[i,:]
-        # )
-        # lincache = init(leadprob, linear_solver)
-        # LinearSolve.solve!(lincache)
-        LinearSolve.set_b(lincache, copy(current_rhs))
+        lincache.b .= current_rhs
         LinearSolve.solve!(lincache)
         Z[i,:] .= lincache.u
     end
 
-    return Geselowitz1989ECGLeadCache(source_op, transfer_op, ∇Njκ∇φₘ_h, ∇Njκ∇φₘ_t, Z, electrode_positions)
+    return Geselowitz1989ECGLeadCache(source_op, transfer_op, φₘ_t, ∇φₘ_t, Z, electrode_positions)
 end
 
 function _add_electrode!(f::AbstractVector{T}, dh::DofHandler, electrode::VertexIndex, weight, lead_field_sym::Symbol) where {T<:Number}
@@ -565,10 +558,10 @@ function _add_electrode!(f::AbstractVector{T}, dh::DofHandler, electrode::Vertex
 end
 
 function update_ecg!(cache::Geselowitz1989ECGLeadCache, φₘ::AbstractVector)
-    # Compute κᵢ∇φₘ on the heart
-    mul!(cache.κ∇φₘ_h, cache.heart_op, φₘ)
     # Transfer κᵢ∇φₘ to the torso
-    transfer!(cache.κ∇φₘ_t, cache.transfer_op, cache.κ∇φₘ_h)
+    transfer!(cache.φₘ_t, cache.transfer_op, φₘ)
+    # Compute κᵢ∇φₘ on the heart
+    mul!(cache.κ∇φₘ_t, cache.source_op, cache.φₘ_t)
     cache.κ∇φₘ_t[isnan.(cache.κ∇φₘ_t)] .= 0.0 # FIXME
     return nothing
 end

test/integration/test_ecg.jl

@@ -1,9 +1,12 @@
+using Thunderbolt, Test, SparseArrays
+import Thunderbolt: to_mesh, OrderedSet
+
 @testset "ECG" begin
-    @testset "Blocks with $geo" for geo in (Hexahedron, Tetrahedron)
-        size = 4.
+    @testset "Blocks with $geo" for geo in (Tetrahedron,Hexahedron)
+        size = 2.
         signal_strength = 0.04
-        nel_heart = (18, 18, 16)
-        nel_torso = (4, 4, 4) .* Int(size) .* 2
+        nel_heart = (6, 6, 6) .* 1
+        nel_torso = (8, 8, 8) .* Int(size)
         ground_vertex = Vec(0.0, 0.0, 0.0)
         electrodes = [
             ground_vertex,
@@ -17,13 +20,16 @@
         electrode_pairs = [[i,1] for i in 2:length(electrodes)]
 
         heart_grid = generate_mesh(geo, nel_heart)
-        # Ferrite.transform_coordinates!(heart_grid, x->Vec{3}(sign.(x) .* x.^2))
+        Ferrite.transform_coordinates!(heart_grid, x->Vec{3}(sign.(x) .* x.^2))
 
-        κ  = SymmetricTensor{2,3,Float64}((1.0, 0, 0, 1.0, 0, 1.0))
-        κᵢ = SymmetricTensor{2,3,Float64}((1.0, 0, 0, 1.0, 0, 1.0))
+        κ  = ConstantCoefficient(SymmetricTensor{2,3,Float64}((1.0, 0, 0, 1.0, 0, 1.0)))
+        κᵢ = AnalyticalCoefficient(
+            (x,t) -> norm(x,Inf) ≤ 1.0 ? SymmetricTensor{2,3,Float64}((1.0, 0, 0, 1.0, 0, 1.0)) : SymmetricTensor{2,3,Float64}((0.0, 0, 0, 0.0, 0, 0.0)),
+            CartesianCoordinateSystem{3}()
+        )
 
         heart_model = TransientDiffusionModel(
-            ConstantCoefficient(κ),
+            κᵢ,
             NoStimulationProtocol(), # Poisoning to detecte if we accidentally touch these
             :φₘ
         )
@@ -38,7 +44,7 @@
 
         op = Thunderbolt.setup_assembled_operator(
             Thunderbolt.BilinearDiffusionIntegrator(
-                ConstantCoefficient(κ),
+                κ,
                 QuadratureRuleCollection(2),
                 :φₘ,
             ),
@@ -48,24 +54,38 @@
         Thunderbolt.update_operator!(op, 0.0) # trigger assembly
 
         torso_grid_ = generate_grid(geo, nel_torso, Vec((-size,-size,-size)), Vec((size,size,size)))
-        addcellset!(torso_grid_, "heart", x->norm(x,1) ≤ 1.0)
-        addcellset!(torso_grid_, "surrounding-tissue", x->norm(x,1) ≥ 1.0)
+        addcellset!(torso_grid_, "heart", x->norm(x,Inf) ≤ 1.0)
+        # addcellset!(torso_grid_, "surrounding-tissue", x->norm(x,Inf) ≥ 1.0)
+        torso_grid_.cellsets["surrounding-tissue"] = OrderedSet([i for i in 1:getncells(torso_grid_) if i ∉ torso_grid_.cellsets["heart"]])
         torso_grid = to_mesh(torso_grid_)
-        torso_fun = semidiscretize(
-            heart_model,
-            FiniteElementDiscretization(
-                Dict(:φₘ => LagrangeCollection{1}()),
-                Dirichlet[],
-                ["heart", "surrounding-tissue"]
-            ),
-            torso_grid
-        )
         u = zeros(Thunderbolt.solution_size(heart_fun))
         plonsey_ecg = Thunderbolt.Plonsey1964ECGGaussCache(op, u)
-        poisson_ecg = Thunderbolt.PoissonECGReconstructionCache(heart_fun, torso_grid, ConstantCoefficient(κᵢ), ConstantCoefficient(κ), electrodes; ground = OrderedSet([Thunderbolt.get_closest_vertex(ground_vertex, torso_grid)]), torso_heart_domain=["heart"])
+        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}
+        )
 
         geselowitz_electrodes = [[electrodes[1], electrodes[i]] for i in 2:length(electrodes)]
-        geselowitz_ecg = Thunderbolt.Geselowitz1989ECGLeadCache(heart_fun, torso_grid, ConstantCoefficient(κᵢ), ConstantCoefficient(κ), geselowitz_electrodes; ground = OrderedSet([Thunderbolt.get_closest_vertex(ground_vertex, torso_grid)]), torso_heart_domain=["heart"])
+        geselowitz_ecg = Thunderbolt.Geselowitz1989ECGLeadCache(
+            heart_fun,
+            torso_grid,
+            κᵢ, κ,
+            geselowitz_electrodes; 
+            ground = OrderedSet([Thunderbolt.get_closest_vertex(ground_vertex, torso_grid)]),
+            torso_heart_domain=["heart"],
+            ipc                  = LagrangeCollection{1}(),
+            qrc                  = QuadratureRuleCollection(3),
+            linear_solver        = Thunderbolt.LinearSolve.UMFPACKFactorization(),
+            system_matrix_type   = SparseMatrixCSC{Float64,Int64}
+        )
 
         @testset "Equilibrium" begin
             u .= 0.0
@@ -115,8 +135,8 @@
             end
         end
 
-        @testset "Planar wave dim=$dim" for dim in 1:3
-            Ferrite.apply_analytical!(u, heart_fun.dh, :φₘ, x->x[dim]^3 )
+        @testset "Planar wave dim=$dim" for dim in 1:1# 1:3
+            Ferrite.apply_analytical!(u, heart_fun.dh, :φₘ, x->x[dim]^3)
 
             @testset "Plonsey1964 xᵢ³" begin
                 Thunderbolt.update_ecg!(plonsey_ecg, u)
@@ -135,7 +155,7 @@
                 @test ecg_vals[1] ≈ 0.0 atol=1e-12 # Ground
                 for dim2 in 1:3
                     if dim2 == dim
-                        @test ecg_vals[2*dim2+1]-ecg_vals[2*dim2] < 0.0
+                        @test ecg_vals[2*dim2+1]-ecg_vals[2*dim2] ≈ -2*0.37 atol=1e-2
                     else
                         @test ecg_vals[2*dim2+1]-ecg_vals[2*dim2] ≈ 0.0 atol=1e-4
                     end
@@ -147,7 +167,7 @@
                 ecg_vals = Thunderbolt.evaluate_ecg(geselowitz_ecg)
                 for dim2 in 1:3
                     if dim2 == dim
-                        @test ecg_vals[2*dim2]-ecg_vals[2*dim2-1] < 0.0
+                        @test ecg_vals[2*dim2]-ecg_vals[2*dim2-1] ≈ -2*0.37 atol=1e-2
                     else
                         @test ecg_vals[2*dim2]-ecg_vals[2*dim2-1] ≈ 0.0 atol=1e-4
                     end
@@ -175,7 +195,7 @@
                     ecg_vals = Thunderbolt.evaluate_ecg(poisson_ecg) 
 
                     if i == dim
-                        @test ecg_vals[2*i+1]-ecg_vals[2*i] > 0.0
+                        @test ecg_vals[2*i+1]-ecg_vals[2*i] ≈ 2*0.37 atol=1e-2
                     else
                         @test ecg_vals[2*i+1]-ecg_vals[2*i] ≈ 0.0 atol=1e-4
                     end
@@ -189,7 +209,7 @@
                     ecg_vals = Thunderbolt.evaluate_ecg(geselowitz_ecg) 
 
                     if i == dim
-                        @test ecg_vals[2*i]-ecg_vals[2*i-1] > 0.0
+                        @test ecg_vals[2*i]-ecg_vals[2*i-1] ≈ 2*0.37 atol=1e-2
                     else
                         @test ecg_vals[2*i]-ecg_vals[2*i-1] ≈ 0.0 atol=1e-4
                     end