From cefb1a951ce6d5d3327d0ba348bd79f026e8e5ee Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 1 Apr 2025 21:51:52 +0200 Subject: [PATCH 1/3] Baseline for FIMH2021 reproducer --- src/modeling/electrophysiology/ecg.jl | 4 +- test/integration/test_ecg.jl | 98 ++++++++++++++++++++++++++- test/test_transfer.jl | 71 ++++++++++++++++++- 3 files changed, 168 insertions(+), 5 deletions(-) diff --git a/src/modeling/electrophysiology/ecg.jl b/src/modeling/electrophysiology/ecg.jl index 935994c25..eabd7590b 100644 --- a/src/modeling/electrophysiology/ecg.jl +++ b/src/modeling/electrophysiology/ecg.jl @@ -52,12 +52,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 diff --git a/test/integration/test_ecg.jl b/test/integration/test_ecg.jl index a0ab58ba2..0184b2441 100644 --- a/test/integration/test_ecg.jl +++ b/test/integration/test_ecg.jl @@ -1,5 +1,6 @@ using Thunderbolt, Test, SparseArrays import Thunderbolt: to_mesh, OrderedSet +using FerriteGmsh @testset "ECG" begin @testset "Blocks with $geo" for geo in (Tetrahedron,Hexahedron) @@ -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,99 @@ import Thunderbolt: to_mesh, OrderedSet end end end + + @testset "FIMH2021" begin + heart_elsize = 5.0#0.1 + torso_elsize = 20.0#5.0 + λ = 0.2 / (0.01*140.0) + + gmsh.initialize() + gmsh.model.add("sis-50") + gmsh.model.occ.add_disk(0.0, 0.0, 0.0, 50.0, 50.0, 1) + heart_handle = gmsh.model.occ.add_disk(0.0, 0.0, 0.0, 2.0, 2.0, 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) + # VTKGridFile("FIMH2021-Debug", grid) do vtk + # Ferrite.write_cellset(vtk, grid) + # end + + κ = ConstantCoefficient(SymmetricTensor{2,2,Float64,3}((λ, 0, λ))) + κᵢ = AnalyticalCoefficient( + (x,t) -> norm(x,2) ≤ 2.0 ? 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(50.0, 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) + 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) + + @test plonsey_vals ≈ poisson_vals + end end diff --git a/test/test_transfer.jl b/test/test_transfer.jl index 484d4d49e..fda2e854c 100644 --- a/test/test_transfer.jl +++ b/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 From 7cdcdbf7ac2c808a8dcbf1c9d1c3f44363ec7bc5 Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 1 Apr 2025 22:56:41 +0200 Subject: [PATCH 2/3] Almost --- src/ferrite-addons/transfer_operators.jl | 4 +-- src/modeling/electrophysiology/ecg.jl | 27 +++++++++--------- test/integration/test_ecg.jl | 36 +++++++++++++++--------- 3 files changed, 38 insertions(+), 29 deletions(-) diff --git a/src/ferrite-addons/transfer_operators.jl b/src/ferrite-addons/transfer_operators.jl index 2e0dacfcc..3d2ce9487 100644 --- a/src/ferrite-addons/transfer_operators.jl +++ b/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." diff --git a/src/modeling/electrophysiology/ecg.jl b/src/modeling/electrophysiology/ecg.jl index eabd7590b..80628b728 100644 --- a/src/modeling/electrophysiology/ecg.jl +++ b/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 diff --git a/test/integration/test_ecg.jl b/test/integration/test_ecg.jl index 0184b2441..0b35ca1a1 100644 --- a/test/integration/test_ecg.jl +++ b/test/integration/test_ecg.jl @@ -5,7 +5,7 @@ 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) @@ -277,14 +277,16 @@ using FerriteGmsh end @testset "FIMH2021" begin - heart_elsize = 5.0#0.1 - torso_elsize = 20.0#5.0 - λ = 0.2 / (0.01*140.0) + 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, 50.0, 50.0, 1) - heart_handle = gmsh.model.occ.add_disk(0.0, 0.0, 0.0, 2.0, 2.0, 2) + 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() @@ -302,13 +304,13 @@ using FerriteGmsh gmsh.finalize() grid = Grid(elements, nodes, cellsets=Dict(["heart" => cellsets["heart"], "torso" => cellsets["torso"]])) torso_grid = heart_grid = to_mesh(grid) - # VTKGridFile("FIMH2021-Debug", grid) do vtk - # Ferrite.write_cellset(vtk, grid) - # end - κ = ConstantCoefficient(SymmetricTensor{2,2,Float64,3}((λ, 0, λ))) + κ = 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) ≤ 2.0 ? SymmetricTensor{2,2,Float64,3}((λ, 0, λ)) : SymmetricTensor{2,2,Float64,3}((0.0, 0.0, 0.0)), + (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}() ) @@ -327,7 +329,7 @@ using FerriteGmsh heart_grid ) - ground_vertex = Vec(50.0, 0.0) + ground_vertex = Vec(torso_radius, 0.0) electrodes = [ ground_vertex, Vec(0., 10.), @@ -350,7 +352,6 @@ using FerriteGmsh plonsey_ecg = Thunderbolt.Plonsey1964ECGGaussCache(op, u, ["heart"]) Thunderbolt.update_ecg!(plonsey_ecg, u) - # plonsey_vals = Thunderbolt.evaluate_ecg(plonsey_ecg) plonsey_vals = Thunderbolt.evaluate_ecg(plonsey_ecg, electrodes[2], λ) poisson_ecg = Thunderbolt.PoissonECGReconstructionCache( @@ -367,7 +368,14 @@ using FerriteGmsh ) 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 - @test plonsey_vals ≈ poisson_vals + # TODO investigate where the 2π comes from + @test plonsey_vals*2π ≈ poisson_vals[1] - poisson_vals[2] rtol = 1e-1 end end From 1e70d51a412280710f014afbacdb752def17646b Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 1 Apr 2025 23:02:01 +0200 Subject: [PATCH 3/3] [skip ci] Make EP04 consistent --- docs/src/literate-tutorials/ep04_geselowitz-ecg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/literate-tutorials/ep04_geselowitz-ecg.jl b/docs/src/literate-tutorials/ep04_geselowitz-ecg.jl index 52f14d751..67488de1f 100644 --- a/docs/src/literate-tutorials/ep04_geselowitz-ecg.jl +++ b/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(