Electrophysiology Tutorial 4: Geselowitz ECG with Monodomain Model
This tutorial shows how to setup ECG problems with monodomain models as the source and compute the QRS complex in a simple toy problem.
Commented Program
using Thunderbolt, LinearAlgebra, StaticArrays, OrdinaryDiffEqOperatorSplittingThe initializer API is not yet finished and hence we deconstruct stuff here manually. Please note that this method is quite fragile w.r.t. to many changes you can make in the code below.
function steady_state_initializer!(u₀, f::GenericSplitFunction)
# TODO cleaner implementation. We need to extract this from the types or via dispatch.
heatfun = f.functions[1]
heat_dofrange = f.solution_indices[1]
odefun = f.functions[2]
ionic_model = odefun.ode
φ₀ = @view u₀[heat_dofrange];
# TODO extraction these via utility functions
dh = heatfun.dh
s₀flat = @view u₀[(ndofs(dh)+1):end];
# Should not be reshape but some array of arrays fun
s₀ = reshape(s₀flat, (ndofs(dh), Thunderbolt.num_states(ionic_model)-1));
default_values = Thunderbolt.default_initial_state(ionic_model)
φ₀ .= default_values[1]
for i ∈ 1:(Thunderbolt.num_states(ionic_model)-1)
s₀[:, i] .= default_values[i+1]
end
endsteady_state_initializer! (generic function with 1 method)We start by defining a custom activation function
Base.@kwdef struct UniformEndocardialActivation <: Function
transmural_depth::Float64 = 0.15
end
function (p::UniformEndocardialActivation)(x::Vec{3}, t)
τᶠ = 0.25
# TODO source for this
if t ≤ 2.0 && x[1] < p.transmural_depth
return 0.5/τᶠ * exp(t/τᶠ)
else
return 0.0
end
end
protocol = Thunderbolt.AnalyticalTransmembraneStimulationProtocol(
AnalyticalCoefficient(
UniformEndocardialActivation(),
CartesianCoordinateSystem{3}()
),
[SVector((-Inf, Inf))],
)AnalyticalTransmembraneStimulationProtocol{AnalyticalCoefficient{Main.var"Main".UniformEndocardialActivation, CartesianCoordinateSystem{3, Float32}}, Float64, Vector{StaticArraysCore.SVector{2, Float64}}}(AnalyticalCoefficient{Main.var"Main".UniformEndocardialActivation, CartesianCoordinateSystem{3, Float32}}(Main.var"Main".UniformEndocardialActivation(0.15), CartesianCoordinateSystem{3, Float32}()), StaticArraysCore.SVector{2, Float64}[[-Inf, Inf]])We also generate both meshes
num_elements_heart = (32,16,16)
heart_mesh = generate_mesh(Tetrahedron, num_elements_heart, Vec((1.5, 1.5, 0.0)), Vec((5.5, 3.5, 2.0)))
num_elements_torso = (56,40,28)
torso_mesh = generate_mesh(Hexahedron, num_elements_torso, Vec((0.0, 0.0, 0.0)), Vec((7.0, 5.0, 3.5)))Thunderbolt.SimpleMesh{3, Hexahedron, Float64} with 980 Hexahedron cells and 1320 nodes
Surface subdomains:
left 70 Hexahedron cells
bottom 140 Hexahedron cells
right 70 Hexahedron cells
back 98 Hexahedron cells
top 140 Hexahedron cells
front 98 Hexahedron cells
Then we place some electrodes and leads.
ground_vertex = Thunderbolt.get_closest_vertex(Vec(0.0, 0.0, 0.0), torso_mesh)
leads = [
[Vec( 0., 0., 1.5), Vec( 7., 0., 1.5)],
[Vec( 3.5, 0., 1.5), Vec( 3.5, 5., 1.5)],
]2-element Vector{Vector{Vec{3, Float64}}}:
[[0.0, 0.0, 1.5], [7.0, 0.0, 1.5]]
[[3.5, 0.0, 1.5], [3.5, 5.0, 1.5]]For our toy problem we use a very simple microstructure.
microstructure = OrthotropicMicrostructureModel(
ConstantCoefficient((Vec(0.0,0.0,1.0))),
ConstantCoefficient((Vec(0.0,1.0,0.0))),
ConstantCoefficient((Vec(1.0,0.0,0.0))),
)OrthotropicMicrostructureModel{ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}}(ConstantCoefficient{Vec{3, Float64}}([0.0, 0.0, 1.0]), ConstantCoefficient{Vec{3, Float64}}([0.0, 1.0, 0.0]), ConstantCoefficient{Vec{3, Float64}}([1.0, 0.0, 0.0]))With the microstructure we setup the diffusion tensor field in spectral form.
κ₁ = 0.17 * 0.62 / (0.17 + 0.62)
κᵣ = 0.019 * 0.24 / (0.019 + 0.24)
diffusion_tensor_field = SpectralTensorCoefficient(
microstructure,
ConstantCoefficient(SVector(κ₁, κᵣ, κᵣ))
)SpectralTensorCoefficient{OrthotropicMicrostructureModel{ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}}, ConstantCoefficient{StaticArraysCore.SVector{3, Float64}}}(OrthotropicMicrostructureModel{ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}, ConstantCoefficient{Vec{3, Float64}}}(ConstantCoefficient{Vec{3, Float64}}([0.0, 0.0, 1.0]), ConstantCoefficient{Vec{3, Float64}}([0.0, 1.0, 0.0]), ConstantCoefficient{Vec{3, Float64}}([1.0, 0.0, 0.0])), ConstantCoefficient{StaticArraysCore.SVector{3, Float64}}([0.13341772151898734, 0.017606177606177605, 0.017606177606177605]))Now we setup our monodomain solver as usual.
cellmodel = Thunderbolt.PCG2019()
heart_model = MonodomainModel(
ConstantCoefficient(1.0),
ConstantCoefficient(1.0),
diffusion_tensor_field,
protocol,
cellmodel,
:φₘ, :s
)
heart_odeform = semidiscretize(
ReactionDiffusionSplit(heart_model),
FiniteElementDiscretization(Dict(:φₘ => LagrangeCollection{1}())),
heart_mesh,
)
u₀ = zeros(Float64, solution_size(heart_odeform))
steady_state_initializer!(u₀, heart_odeform)
dt₀ = 0.01
dtvis = 0.5
Tₘₐₓ = 50.0
tspan = (0.0, Tₘₐₓ)
problem = OperatorSplittingProblem(heart_odeform, u₀, tspan)
timestepper = LieTrotterGodunov((
BackwardEulerSolver(),
ForwardEulerCellSolver(),
))
integrator = init(problem, timestepper, dt=dt₀, verbose=true)t: 0.0
u: 1575-element Vector{Float64}:
-85.0
-85.0
-85.0
-85.0
-85.0
-85.0
-85.0
-85.0
-85.0
-85.0
⋮
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137
0.00012530739244227137Now that the time integrator is ready we setup the ECG problem.
torso_mesh_κᵢ = ConstantCoefficient(1.0)
torso_mesh_κ = ConstantCoefficient(1.0)ConstantCoefficient{Float64}(1.0)geselowitz_ecg = Thunderbolt.Geselowitz1989ECGLeadCache(
heart_odeform,
torso_mesh,
torso_mesh_κᵢ,
torso_mesh_κ,
leads;
ground = Thunderbolt.OrderedSet([ground_vertex])
)Thunderbolt.Geselowitz1989ECGLeadCache{Matrix{Float64}, Thunderbolt.AssembledBilinearOperator{Thunderbolt.ThreadedSparseMatrixCSR{Float64, Int64}, Thunderbolt.ThreadedSparseMatrixCSR{Float64, Int64}, Thunderbolt.BilinearDiffusionIntegrator{ConstantCoefficient{Float64}, QuadratureRuleCollection{2}}, DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}, Thunderbolt.SequentialAssemblyStrategyCache{SequentialCPUDevice{Float64, Int64}}}, Thunderbolt.NodalIntergridInterpolation{PointEvalHandler{Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}, Float64}, DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}, DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}, Vector{Float64}, Vector{Vector{VertexIndex}}}(Thunderbolt.AssembledBilinearOperator{Thunderbolt.ThreadedSparseMatrixCSR{Float64, Int64}, Thunderbolt.ThreadedSparseMatrixCSR{Float64, Int64}, Thunderbolt.BilinearDiffusionIntegrator{ConstantCoefficient{Float64}, QuadratureRuleCollection{2}}, DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}, Thunderbolt.SequentialAssemblyStrategyCache{SequentialCPUDevice{Float64, Int64}}}([-0.16666666666666663 1.0842021724855044e-17 … 0.0 0.0; 1.0842021724855044e-17 -0.33333333333333326 … 0.0 0.0; … ; 0.0 0.0 … -0.33333333333333337 2.949029909160572e-17; 0.0 0.0 … 2.949029909160572e-17 -0.1666666666666669], [-0.16666666666666663 1.0842021724855044e-17 … 0.0 0.0; 1.0842021724855044e-17 -0.33333333333333326 … 0.0 0.0; … ; 0.0 0.0 … -0.33333333333333337 2.949029909160572e-17; 0.0 0.0 … 2.949029909160572e-17 -0.1666666666666669], Thunderbolt.BilinearDiffusionIntegrator{ConstantCoefficient{Float64}, QuadratureRuleCollection{2}}(ConstantCoefficient{Float64}(1.0), QuadratureRuleCollection{2}(), :φₘ), DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}(SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}[SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}(DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}(#= circular reference @-3 =#), OrderedCollections.OrderedSet{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 971, 972, 973, 974, 975, 976, 977, 978, 979, 980]), [:Z], Interpolation[Lagrange{RefHexahedron, 1}()], Int64[], 8)], [:Z], [1, 2, 3, 4, 5, 6, 7, 8, 2, 9 … 1319, 1318, 1139, 1140, 1155, 1154, 1304, 1305, 1320, 1319], [1, 9, 17, 25, 33, 41, 49, 57, 65, 73 … 7761, 7769, 7777, 7785, 7793, 7801, 7809, 7817, 7825, 7833], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}(Grid{3, Hexahedron, Float64}(Hexahedron[Hexahedron((1, 2, 17, 16, 166, 167, 182, 181)), Hexahedron((2, 3, 18, 17, 167, 168, 183, 182)), Hexahedron((3, 4, 19, 18, 168, 169, 184, 183)), Hexahedron((4, 5, 20, 19, 169, 170, 185, 184)), Hexahedron((5, 6, 21, 20, 170, 171, 186, 185)), Hexahedron((6, 7, 22, 21, 171, 172, 187, 186)), Hexahedron((7, 8, 23, 22, 172, 173, 188, 187)), Hexahedron((8, 9, 24, 23, 173, 174, 189, 188)), Hexahedron((9, 10, 25, 24, 174, 175, 190, 189)), Hexahedron((10, 11, 26, 25, 175, 176, 191, 190)) … Hexahedron((1130, 1131, 1146, 1145, 1295, 1296, 1311, 1310)), Hexahedron((1131, 1132, 1147, 1146, 1296, 1297, 1312, 1311)), Hexahedron((1132, 1133, 1148, 1147, 1297, 1298, 1313, 1312)), Hexahedron((1133, 1134, 1149, 1148, 1298, 1299, 1314, 1313)), Hexahedron((1134, 1135, 1150, 1149, 1299, 1300, 1315, 1314)), Hexahedron((1135, 1136, 1151, 1150, 1300, 1301, 1316, 1315)), Hexahedron((1136, 1137, 1152, 1151, 1301, 1302, 1317, 1316)), Hexahedron((1137, 1138, 1153, 1152, 1302, 1303, 1318, 1317)), Hexahedron((1138, 1139, 1154, 1153, 1303, 1304, 1319, 1318)), Hexahedron((1139, 1140, 1155, 1154, 1304, 1305, 1320, 1319))], Node{3, Float64}[Node{3, Float64}([0.0, 0.0, 0.0]), Node{3, Float64}([0.5, 0.0, 0.0]), Node{3, Float64}([1.0, 0.0, 0.0]), Node{3, Float64}([1.5, 0.0, 0.0]), Node{3, Float64}([2.0, 0.0, 0.0]), Node{3, Float64}([2.5, 0.0, 0.0]), Node{3, Float64}([3.0, 0.0, 0.0]), Node{3, Float64}([3.5, 0.0, 0.0]), Node{3, Float64}([4.0, 0.0, 0.0]), Node{3, Float64}([4.5, 0.0, 0.0]) … Node{3, Float64}([2.5, 5.0, 3.5]), Node{3, Float64}([3.0, 5.0, 3.5]), Node{3, Float64}([3.5, 5.0, 3.5]), Node{3, Float64}([4.0, 5.0, 3.5]), Node{3, Float64}([4.5, 5.0, 3.5]), Node{3, Float64}([5.0, 5.0, 3.5]), Node{3, Float64}([5.5, 5.0, 3.5]), Node{3, Float64}([6.0, 5.0, 3.5]), Node{3, Float64}([6.5, 5.0, 3.5]), Node{3, Float64}([7.0, 5.0, 3.5])], Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{FacetIndex}}("left" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 5)), FacetIndex((15, 5)), FacetIndex((29, 5)), FacetIndex((43, 5)), FacetIndex((57, 5)), FacetIndex((71, 5)), FacetIndex((85, 5)), FacetIndex((99, 5)), FacetIndex((113, 5)), FacetIndex((127, 5)) … FacetIndex((841, 5)), FacetIndex((855, 5)), FacetIndex((869, 5)), FacetIndex((883, 5)), FacetIndex((897, 5)), FacetIndex((911, 5)), FacetIndex((925, 5)), FacetIndex((939, 5)), FacetIndex((953, 5)), FacetIndex((967, 5))]), "bottom" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((2, 1)), FacetIndex((3, 1)), FacetIndex((4, 1)), FacetIndex((5, 1)), FacetIndex((6, 1)), FacetIndex((7, 1)), FacetIndex((8, 1)), FacetIndex((9, 1)), FacetIndex((10, 1)) … FacetIndex((131, 1)), FacetIndex((132, 1)), FacetIndex((133, 1)), FacetIndex((134, 1)), FacetIndex((135, 1)), FacetIndex((136, 1)), FacetIndex((137, 1)), FacetIndex((138, 1)), FacetIndex((139, 1)), FacetIndex((140, 1))]), "right" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((14, 3)), FacetIndex((28, 3)), FacetIndex((42, 3)), FacetIndex((56, 3)), FacetIndex((70, 3)), FacetIndex((84, 3)), FacetIndex((98, 3)), FacetIndex((112, 3)), FacetIndex((126, 3)), FacetIndex((140, 3)) … FacetIndex((854, 3)), FacetIndex((868, 3)), FacetIndex((882, 3)), FacetIndex((896, 3)), FacetIndex((910, 3)), FacetIndex((924, 3)), FacetIndex((938, 3)), FacetIndex((952, 3)), FacetIndex((966, 3)), FacetIndex((980, 3))]), "back" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((127, 4)), FacetIndex((128, 4)), FacetIndex((129, 4)), FacetIndex((130, 4)), FacetIndex((131, 4)), FacetIndex((132, 4)), FacetIndex((133, 4)), FacetIndex((134, 4)), FacetIndex((135, 4)), FacetIndex((136, 4)) … FacetIndex((971, 4)), FacetIndex((972, 4)), FacetIndex((973, 4)), FacetIndex((974, 4)), FacetIndex((975, 4)), FacetIndex((976, 4)), FacetIndex((977, 4)), FacetIndex((978, 4)), FacetIndex((979, 4)), FacetIndex((980, 4))]), "top" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((841, 6)), FacetIndex((842, 6)), FacetIndex((843, 6)), FacetIndex((844, 6)), FacetIndex((845, 6)), FacetIndex((846, 6)), FacetIndex((847, 6)), FacetIndex((848, 6)), FacetIndex((849, 6)), FacetIndex((850, 6)) … FacetIndex((971, 6)), FacetIndex((972, 6)), FacetIndex((973, 6)), FacetIndex((974, 6)), FacetIndex((975, 6)), FacetIndex((976, 6)), FacetIndex((977, 6)), FacetIndex((978, 6)), FacetIndex((979, 6)), FacetIndex((980, 6))]), "front" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 2)), FacetIndex((2, 2)), FacetIndex((3, 2)), FacetIndex((4, 2)), FacetIndex((5, 2)), FacetIndex((6, 2)), FacetIndex((7, 2)), FacetIndex((8, 2)), FacetIndex((9, 2)), FacetIndex((10, 2)) … FacetIndex((845, 2)), FacetIndex((846, 2)), FacetIndex((847, 2)), FacetIndex((848, 2)), FacetIndex((849, 2)), FacetIndex((850, 2)), FacetIndex((851, 2)), FacetIndex((852, 2)), FacetIndex((853, 2)), FacetIndex((854, 2))])), Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()), OrderedCollections.OrderedDict{Tuple{Int64, Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Tuple{Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Int64, Int64}(), OrderedCollections.OrderedDict{DataType, Int64}(Hexahedron => 980), OrderedCollections.OrderedDict{String, Thunderbolt.VolumetricSubdomainDesriptor}("" => Thunderbolt.VolumetricSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{CellIndex}}(Hexahedron => OrderedCollections.OrderedSet{CellIndex}(CellIndex[CellIndex(1), CellIndex(2), CellIndex(3), CellIndex(4), CellIndex(5), CellIndex(6), CellIndex(7), CellIndex(8), CellIndex(9), CellIndex(10) … CellIndex(971), CellIndex(972), CellIndex(973), CellIndex(974), CellIndex(975), CellIndex(976), CellIndex(977), CellIndex(978), CellIndex(979), CellIndex(980)])))), OrderedCollections.OrderedDict{String, Thunderbolt.SurfaceSubdomainDesriptor}("left" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 5)), FacetIndex((15, 5)), FacetIndex((29, 5)), FacetIndex((43, 5)), FacetIndex((57, 5)), FacetIndex((71, 5)), FacetIndex((85, 5)), FacetIndex((99, 5)), FacetIndex((113, 5)), FacetIndex((127, 5)) … FacetIndex((841, 5)), FacetIndex((855, 5)), FacetIndex((869, 5)), FacetIndex((883, 5)), FacetIndex((897, 5)), FacetIndex((911, 5)), FacetIndex((925, 5)), FacetIndex((939, 5)), FacetIndex((953, 5)), FacetIndex((967, 5))]))), "bottom" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 1)), FacetIndex((2, 1)), FacetIndex((3, 1)), FacetIndex((4, 1)), FacetIndex((5, 1)), FacetIndex((6, 1)), FacetIndex((7, 1)), FacetIndex((8, 1)), FacetIndex((9, 1)), FacetIndex((10, 1)) … FacetIndex((131, 1)), FacetIndex((132, 1)), FacetIndex((133, 1)), FacetIndex((134, 1)), FacetIndex((135, 1)), FacetIndex((136, 1)), FacetIndex((137, 1)), FacetIndex((138, 1)), FacetIndex((139, 1)), FacetIndex((140, 1))]))), "right" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((14, 3)), FacetIndex((28, 3)), FacetIndex((42, 3)), FacetIndex((56, 3)), FacetIndex((70, 3)), FacetIndex((84, 3)), FacetIndex((98, 3)), FacetIndex((112, 3)), FacetIndex((126, 3)), FacetIndex((140, 3)) … FacetIndex((854, 3)), FacetIndex((868, 3)), FacetIndex((882, 3)), FacetIndex((896, 3)), FacetIndex((910, 3)), FacetIndex((924, 3)), FacetIndex((938, 3)), FacetIndex((952, 3)), FacetIndex((966, 3)), FacetIndex((980, 3))]))), "back" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((127, 4)), FacetIndex((128, 4)), FacetIndex((129, 4)), FacetIndex((130, 4)), FacetIndex((131, 4)), FacetIndex((132, 4)), FacetIndex((133, 4)), FacetIndex((134, 4)), FacetIndex((135, 4)), FacetIndex((136, 4)) … FacetIndex((971, 4)), FacetIndex((972, 4)), FacetIndex((973, 4)), FacetIndex((974, 4)), FacetIndex((975, 4)), FacetIndex((976, 4)), FacetIndex((977, 4)), FacetIndex((978, 4)), FacetIndex((979, 4)), FacetIndex((980, 4))]))), "top" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((841, 6)), FacetIndex((842, 6)), FacetIndex((843, 6)), FacetIndex((844, 6)), FacetIndex((845, 6)), FacetIndex((846, 6)), FacetIndex((847, 6)), FacetIndex((848, 6)), FacetIndex((849, 6)), FacetIndex((850, 6)) … FacetIndex((971, 6)), FacetIndex((972, 6)), FacetIndex((973, 6)), FacetIndex((974, 6)), FacetIndex((975, 6)), FacetIndex((976, 6)), FacetIndex((977, 6)), FacetIndex((978, 6)), FacetIndex((979, 6)), FacetIndex((980, 6))]))), "front" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 2)), FacetIndex((2, 2)), FacetIndex((3, 2)), FacetIndex((4, 2)), FacetIndex((5, 2)), FacetIndex((6, 2)), FacetIndex((7, 2)), FacetIndex((8, 2)), FacetIndex((9, 2)), FacetIndex((10, 2)) … FacetIndex((845, 2)), FacetIndex((846, 2)), FacetIndex((847, 2)), FacetIndex((848, 2)), FacetIndex((849, 2)), FacetIndex((850, 2)), FacetIndex((851, 2)), FacetIndex((852, 2)), FacetIndex((853, 2)), FacetIndex((854, 2))])))), OrderedCollections.OrderedDict{String, Thunderbolt.InterfaceSubdomainDesriptor}()), 1320), Thunderbolt.SequentialAssemblyStrategyCache{SequentialCPUDevice{Float64, Int64}}(SequentialCPUDevice{Float64, Int64}())), Thunderbolt.NodalIntergridInterpolation{PointEvalHandler{Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}, Float64}, DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}, DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}(PointEvalHandler{Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}, Float64}(Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}(Grid{3, Tetrahedron, Float64}(Tetrahedron[Tetrahedron((1, 2, 10, 55)), Tetrahedron((1, 46, 2, 55)), Tetrahedron((2, 11, 10, 55)), Tetrahedron((2, 56, 11, 55)), Tetrahedron((2, 46, 47, 55)), Tetrahedron((2, 47, 56, 55)), Tetrahedron((2, 3, 11, 56)), Tetrahedron((2, 47, 3, 56)), Tetrahedron((3, 12, 11, 56)), Tetrahedron((3, 57, 12, 56)) … Tetrahedron((170, 179, 178, 223)), Tetrahedron((170, 224, 179, 223)), Tetrahedron((170, 214, 215, 223)), Tetrahedron((170, 215, 224, 223)), Tetrahedron((170, 171, 179, 224)), Tetrahedron((170, 215, 171, 224)), Tetrahedron((171, 180, 179, 224)), Tetrahedron((171, 225, 180, 224)), Tetrahedron((171, 215, 216, 224)), Tetrahedron((171, 216, 225, 224))], Node{3, Float64}[Node{3, Float64}([1.5, 1.5, 0.0]), Node{3, Float64}([2.0, 1.5, 0.0]), Node{3, Float64}([2.5, 1.5, 0.0]), Node{3, Float64}([3.0, 1.5, 0.0]), Node{3, Float64}([3.5, 1.5, 0.0]), Node{3, Float64}([4.0, 1.5, 0.0]), Node{3, Float64}([4.5, 1.5, 0.0]), Node{3, Float64}([5.0, 1.5, 0.0]), Node{3, Float64}([5.5, 1.5, 0.0]), Node{3, Float64}([1.5, 2.0, 0.0]) … Node{3, Float64}([5.5, 3.0, 2.0]), Node{3, Float64}([1.5, 3.5, 2.0]), Node{3, Float64}([2.0, 3.5, 2.0]), Node{3, Float64}([2.5, 3.5, 2.0]), Node{3, Float64}([3.0, 3.5, 2.0]), Node{3, Float64}([3.5, 3.5, 2.0]), Node{3, Float64}([4.0, 3.5, 2.0]), Node{3, Float64}([4.5, 3.5, 2.0]), Node{3, Float64}([5.0, 3.5, 2.0]), Node{3, Float64}([5.5, 3.5, 2.0])], Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{FacetIndex}}("left" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 4)), FacetIndex((2, 2)), FacetIndex((49, 4)), FacetIndex((50, 2)), FacetIndex((97, 4)), FacetIndex((98, 2)), FacetIndex((145, 4)), FacetIndex((146, 2)), FacetIndex((193, 4)), FacetIndex((194, 2)) … FacetIndex((529, 4)), FacetIndex((530, 2)), FacetIndex((577, 4)), FacetIndex((578, 2)), FacetIndex((625, 4)), FacetIndex((626, 2)), FacetIndex((673, 4)), FacetIndex((674, 2)), FacetIndex((721, 4)), FacetIndex((722, 2))]), "bottom" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((7, 1)), FacetIndex((9, 1)), FacetIndex((13, 1)), FacetIndex((15, 1)), FacetIndex((19, 1)), FacetIndex((21, 1)), FacetIndex((25, 1)), FacetIndex((27, 1)) … FacetIndex((163, 1)), FacetIndex((165, 1)), FacetIndex((169, 1)), FacetIndex((171, 1)), FacetIndex((175, 1)), FacetIndex((177, 1)), FacetIndex((181, 1)), FacetIndex((183, 1)), FacetIndex((187, 1)), FacetIndex((189, 1))]), "right" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((46, 1)), FacetIndex((48, 1)), FacetIndex((94, 1)), FacetIndex((96, 1)), FacetIndex((142, 1)), FacetIndex((144, 1)), FacetIndex((190, 1)), FacetIndex((192, 1)), FacetIndex((238, 1)), FacetIndex((240, 1)) … FacetIndex((574, 1)), FacetIndex((576, 1)), FacetIndex((622, 1)), FacetIndex((624, 1)), FacetIndex((670, 1)), FacetIndex((672, 1)), FacetIndex((718, 1)), FacetIndex((720, 1)), FacetIndex((766, 1)), FacetIndex((768, 1))]), "back" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((147, 3)), FacetIndex((148, 3)), FacetIndex((153, 3)), FacetIndex((154, 3)), FacetIndex((159, 3)), FacetIndex((160, 3)), FacetIndex((165, 3)), FacetIndex((166, 3)), FacetIndex((171, 3)), FacetIndex((172, 3)) … FacetIndex((741, 3)), FacetIndex((742, 3)), FacetIndex((747, 3)), FacetIndex((748, 3)), FacetIndex((753, 3)), FacetIndex((754, 3)), FacetIndex((759, 3)), FacetIndex((760, 3)), FacetIndex((765, 3)), FacetIndex((766, 3))]), "top" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((581, 3)), FacetIndex((582, 3)), FacetIndex((587, 3)), FacetIndex((588, 3)), FacetIndex((593, 3)), FacetIndex((594, 3)), FacetIndex((599, 3)), FacetIndex((600, 3)), FacetIndex((605, 3)), FacetIndex((606, 3)) … FacetIndex((743, 3)), FacetIndex((744, 3)), FacetIndex((749, 3)), FacetIndex((750, 3)), FacetIndex((755, 3)), FacetIndex((756, 3)), FacetIndex((761, 3)), FacetIndex((762, 3)), FacetIndex((767, 3)), FacetIndex((768, 3))]), "front" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((2, 1)), FacetIndex((5, 1)), FacetIndex((8, 1)), FacetIndex((11, 1)), FacetIndex((14, 1)), FacetIndex((17, 1)), FacetIndex((20, 1)), FacetIndex((23, 1)), FacetIndex((26, 1)), FacetIndex((29, 1)) … FacetIndex((596, 1)), FacetIndex((599, 1)), FacetIndex((602, 1)), FacetIndex((605, 1)), FacetIndex((608, 1)), FacetIndex((611, 1)), FacetIndex((614, 1)), FacetIndex((617, 1)), FacetIndex((620, 1)), FacetIndex((623, 1))])), Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()), OrderedCollections.OrderedDict{Tuple{Int64, Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Tuple{Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Int64, Int64}(), OrderedCollections.OrderedDict{DataType, Int64}(Tetrahedron => 768), OrderedCollections.OrderedDict{String, Thunderbolt.VolumetricSubdomainDesriptor}("" => Thunderbolt.VolumetricSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{CellIndex}}(Tetrahedron => OrderedCollections.OrderedSet{CellIndex}(CellIndex[CellIndex(1), CellIndex(2), CellIndex(3), CellIndex(4), CellIndex(5), CellIndex(6), CellIndex(7), CellIndex(8), CellIndex(9), CellIndex(10) … CellIndex(759), CellIndex(760), CellIndex(761), CellIndex(762), CellIndex(763), CellIndex(764), CellIndex(765), CellIndex(766), CellIndex(767), CellIndex(768)])))), OrderedCollections.OrderedDict{String, Thunderbolt.SurfaceSubdomainDesriptor}("left" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 4)), FacetIndex((2, 2)), FacetIndex((49, 4)), FacetIndex((50, 2)), FacetIndex((97, 4)), FacetIndex((98, 2)), FacetIndex((145, 4)), FacetIndex((146, 2)), FacetIndex((193, 4)), FacetIndex((194, 2)) … FacetIndex((529, 4)), FacetIndex((530, 2)), FacetIndex((577, 4)), FacetIndex((578, 2)), FacetIndex((625, 4)), FacetIndex((626, 2)), FacetIndex((673, 4)), FacetIndex((674, 2)), FacetIndex((721, 4)), FacetIndex((722, 2))]))), "bottom" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((7, 1)), FacetIndex((9, 1)), FacetIndex((13, 1)), FacetIndex((15, 1)), FacetIndex((19, 1)), FacetIndex((21, 1)), FacetIndex((25, 1)), FacetIndex((27, 1)) … FacetIndex((163, 1)), FacetIndex((165, 1)), FacetIndex((169, 1)), FacetIndex((171, 1)), FacetIndex((175, 1)), FacetIndex((177, 1)), FacetIndex((181, 1)), FacetIndex((183, 1)), FacetIndex((187, 1)), FacetIndex((189, 1))]))), "right" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((46, 1)), FacetIndex((48, 1)), FacetIndex((94, 1)), FacetIndex((96, 1)), FacetIndex((142, 1)), FacetIndex((144, 1)), FacetIndex((190, 1)), FacetIndex((192, 1)), FacetIndex((238, 1)), FacetIndex((240, 1)) … FacetIndex((574, 1)), FacetIndex((576, 1)), FacetIndex((622, 1)), FacetIndex((624, 1)), FacetIndex((670, 1)), FacetIndex((672, 1)), FacetIndex((718, 1)), FacetIndex((720, 1)), FacetIndex((766, 1)), FacetIndex((768, 1))]))), "back" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((147, 3)), FacetIndex((148, 3)), FacetIndex((153, 3)), FacetIndex((154, 3)), FacetIndex((159, 3)), FacetIndex((160, 3)), FacetIndex((165, 3)), FacetIndex((166, 3)), FacetIndex((171, 3)), FacetIndex((172, 3)) … FacetIndex((741, 3)), FacetIndex((742, 3)), FacetIndex((747, 3)), FacetIndex((748, 3)), FacetIndex((753, 3)), FacetIndex((754, 3)), FacetIndex((759, 3)), FacetIndex((760, 3)), FacetIndex((765, 3)), FacetIndex((766, 3))]))), "top" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((581, 3)), FacetIndex((582, 3)), FacetIndex((587, 3)), FacetIndex((588, 3)), FacetIndex((593, 3)), FacetIndex((594, 3)), FacetIndex((599, 3)), FacetIndex((600, 3)), FacetIndex((605, 3)), FacetIndex((606, 3)) … FacetIndex((743, 3)), FacetIndex((744, 3)), FacetIndex((749, 3)), FacetIndex((750, 3)), FacetIndex((755, 3)), FacetIndex((756, 3)), FacetIndex((761, 3)), FacetIndex((762, 3)), FacetIndex((767, 3)), FacetIndex((768, 3))]))), "front" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((2, 1)), FacetIndex((5, 1)), FacetIndex((8, 1)), FacetIndex((11, 1)), FacetIndex((14, 1)), FacetIndex((17, 1)), FacetIndex((20, 1)), FacetIndex((23, 1)), FacetIndex((26, 1)), FacetIndex((29, 1)) … FacetIndex((596, 1)), FacetIndex((599, 1)), FacetIndex((602, 1)), FacetIndex((605, 1)), FacetIndex((608, 1)), FacetIndex((611, 1)), FacetIndex((614, 1)), FacetIndex((617, 1)), FacetIndex((620, 1)), FacetIndex((623, 1))])))), OrderedCollections.OrderedDict{String, Thunderbolt.InterfaceSubdomainDesriptor}()), Union{Nothing, Int64}[nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing … nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing], Union{Nothing, Vec{1, Float64}, Vec{2, Float64}, Vec{3, Float64}}[nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing … nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]), DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}(SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}}[SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}}(DofHandler{3, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}}(#= circular reference @-3 =#), OrderedCollections.OrderedSet{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 759, 760, 761, 762, 763, 764, 765, 766, 767, 768]), [:φₘ], Interpolation[Lagrange{RefTetrahedron, 1}()], Int64[], 4)], [:φₘ], [1, 2, 3, 4, 1, 5, 2, 4, 2, 6 … 180, 224, 171, 215, 216, 224, 171, 216, 225, 224], [1, 5, 9, 13, 17, 21, 25, 29, 33, 37 … 3033, 3037, 3041, 3045, 3049, 3053, 3057, 3061, 3065, 3069], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, Thunderbolt.SimpleMesh{3, Tetrahedron, Float64}(Grid{3, Tetrahedron, Float64}(Tetrahedron[Tetrahedron((1, 2, 10, 55)), Tetrahedron((1, 46, 2, 55)), Tetrahedron((2, 11, 10, 55)), Tetrahedron((2, 56, 11, 55)), Tetrahedron((2, 46, 47, 55)), Tetrahedron((2, 47, 56, 55)), Tetrahedron((2, 3, 11, 56)), Tetrahedron((2, 47, 3, 56)), Tetrahedron((3, 12, 11, 56)), Tetrahedron((3, 57, 12, 56)) … Tetrahedron((170, 179, 178, 223)), Tetrahedron((170, 224, 179, 223)), Tetrahedron((170, 214, 215, 223)), Tetrahedron((170, 215, 224, 223)), Tetrahedron((170, 171, 179, 224)), Tetrahedron((170, 215, 171, 224)), Tetrahedron((171, 180, 179, 224)), Tetrahedron((171, 225, 180, 224)), Tetrahedron((171, 215, 216, 224)), Tetrahedron((171, 216, 225, 224))], Node{3, Float64}[Node{3, Float64}([1.5, 1.5, 0.0]), Node{3, Float64}([2.0, 1.5, 0.0]), Node{3, Float64}([2.5, 1.5, 0.0]), Node{3, Float64}([3.0, 1.5, 0.0]), Node{3, Float64}([3.5, 1.5, 0.0]), Node{3, Float64}([4.0, 1.5, 0.0]), Node{3, Float64}([4.5, 1.5, 0.0]), Node{3, Float64}([5.0, 1.5, 0.0]), Node{3, Float64}([5.5, 1.5, 0.0]), Node{3, Float64}([1.5, 2.0, 0.0]) … Node{3, Float64}([5.5, 3.0, 2.0]), Node{3, Float64}([1.5, 3.5, 2.0]), Node{3, Float64}([2.0, 3.5, 2.0]), Node{3, Float64}([2.5, 3.5, 2.0]), Node{3, Float64}([3.0, 3.5, 2.0]), Node{3, Float64}([3.5, 3.5, 2.0]), Node{3, Float64}([4.0, 3.5, 2.0]), Node{3, Float64}([4.5, 3.5, 2.0]), Node{3, Float64}([5.0, 3.5, 2.0]), Node{3, Float64}([5.5, 3.5, 2.0])], Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{FacetIndex}}("left" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 4)), FacetIndex((2, 2)), FacetIndex((49, 4)), FacetIndex((50, 2)), FacetIndex((97, 4)), FacetIndex((98, 2)), FacetIndex((145, 4)), FacetIndex((146, 2)), FacetIndex((193, 4)), FacetIndex((194, 2)) … FacetIndex((529, 4)), FacetIndex((530, 2)), FacetIndex((577, 4)), FacetIndex((578, 2)), FacetIndex((625, 4)), FacetIndex((626, 2)), FacetIndex((673, 4)), FacetIndex((674, 2)), FacetIndex((721, 4)), FacetIndex((722, 2))]), "bottom" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((7, 1)), FacetIndex((9, 1)), FacetIndex((13, 1)), FacetIndex((15, 1)), FacetIndex((19, 1)), FacetIndex((21, 1)), FacetIndex((25, 1)), FacetIndex((27, 1)) … FacetIndex((163, 1)), FacetIndex((165, 1)), FacetIndex((169, 1)), FacetIndex((171, 1)), FacetIndex((175, 1)), FacetIndex((177, 1)), FacetIndex((181, 1)), FacetIndex((183, 1)), FacetIndex((187, 1)), FacetIndex((189, 1))]), "right" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((46, 1)), FacetIndex((48, 1)), FacetIndex((94, 1)), FacetIndex((96, 1)), FacetIndex((142, 1)), FacetIndex((144, 1)), FacetIndex((190, 1)), FacetIndex((192, 1)), FacetIndex((238, 1)), FacetIndex((240, 1)) … FacetIndex((574, 1)), FacetIndex((576, 1)), FacetIndex((622, 1)), FacetIndex((624, 1)), FacetIndex((670, 1)), FacetIndex((672, 1)), FacetIndex((718, 1)), FacetIndex((720, 1)), FacetIndex((766, 1)), FacetIndex((768, 1))]), "back" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((147, 3)), FacetIndex((148, 3)), FacetIndex((153, 3)), FacetIndex((154, 3)), FacetIndex((159, 3)), FacetIndex((160, 3)), FacetIndex((165, 3)), FacetIndex((166, 3)), FacetIndex((171, 3)), FacetIndex((172, 3)) … FacetIndex((741, 3)), FacetIndex((742, 3)), FacetIndex((747, 3)), FacetIndex((748, 3)), FacetIndex((753, 3)), FacetIndex((754, 3)), FacetIndex((759, 3)), FacetIndex((760, 3)), FacetIndex((765, 3)), FacetIndex((766, 3))]), "top" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((581, 3)), FacetIndex((582, 3)), FacetIndex((587, 3)), FacetIndex((588, 3)), FacetIndex((593, 3)), FacetIndex((594, 3)), FacetIndex((599, 3)), FacetIndex((600, 3)), FacetIndex((605, 3)), FacetIndex((606, 3)) … FacetIndex((743, 3)), FacetIndex((744, 3)), FacetIndex((749, 3)), FacetIndex((750, 3)), FacetIndex((755, 3)), FacetIndex((756, 3)), FacetIndex((761, 3)), FacetIndex((762, 3)), FacetIndex((767, 3)), FacetIndex((768, 3))]), "front" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((2, 1)), FacetIndex((5, 1)), FacetIndex((8, 1)), FacetIndex((11, 1)), FacetIndex((14, 1)), FacetIndex((17, 1)), FacetIndex((20, 1)), FacetIndex((23, 1)), FacetIndex((26, 1)), FacetIndex((29, 1)) … FacetIndex((596, 1)), FacetIndex((599, 1)), FacetIndex((602, 1)), FacetIndex((605, 1)), FacetIndex((608, 1)), FacetIndex((611, 1)), FacetIndex((614, 1)), FacetIndex((617, 1)), FacetIndex((620, 1)), FacetIndex((623, 1))])), Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()), OrderedCollections.OrderedDict{Tuple{Int64, Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Tuple{Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Int64, Int64}(), OrderedCollections.OrderedDict{DataType, Int64}(Tetrahedron => 768), OrderedCollections.OrderedDict{String, Thunderbolt.VolumetricSubdomainDesriptor}("" => Thunderbolt.VolumetricSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{CellIndex}}(Tetrahedron => OrderedCollections.OrderedSet{CellIndex}(CellIndex[CellIndex(1), CellIndex(2), CellIndex(3), CellIndex(4), CellIndex(5), CellIndex(6), CellIndex(7), CellIndex(8), CellIndex(9), CellIndex(10) … CellIndex(759), CellIndex(760), CellIndex(761), CellIndex(762), CellIndex(763), CellIndex(764), CellIndex(765), CellIndex(766), CellIndex(767), CellIndex(768)])))), OrderedCollections.OrderedDict{String, Thunderbolt.SurfaceSubdomainDesriptor}("left" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 4)), FacetIndex((2, 2)), FacetIndex((49, 4)), FacetIndex((50, 2)), FacetIndex((97, 4)), FacetIndex((98, 2)), FacetIndex((145, 4)), FacetIndex((146, 2)), FacetIndex((193, 4)), FacetIndex((194, 2)) … FacetIndex((529, 4)), FacetIndex((530, 2)), FacetIndex((577, 4)), FacetIndex((578, 2)), FacetIndex((625, 4)), FacetIndex((626, 2)), FacetIndex((673, 4)), FacetIndex((674, 2)), FacetIndex((721, 4)), FacetIndex((722, 2))]))), "bottom" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((7, 1)), FacetIndex((9, 1)), FacetIndex((13, 1)), FacetIndex((15, 1)), FacetIndex((19, 1)), FacetIndex((21, 1)), FacetIndex((25, 1)), FacetIndex((27, 1)) … FacetIndex((163, 1)), FacetIndex((165, 1)), FacetIndex((169, 1)), FacetIndex((171, 1)), FacetIndex((175, 1)), FacetIndex((177, 1)), FacetIndex((181, 1)), FacetIndex((183, 1)), FacetIndex((187, 1)), FacetIndex((189, 1))]))), "right" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((46, 1)), FacetIndex((48, 1)), FacetIndex((94, 1)), FacetIndex((96, 1)), FacetIndex((142, 1)), FacetIndex((144, 1)), FacetIndex((190, 1)), FacetIndex((192, 1)), FacetIndex((238, 1)), FacetIndex((240, 1)) … FacetIndex((574, 1)), FacetIndex((576, 1)), FacetIndex((622, 1)), FacetIndex((624, 1)), FacetIndex((670, 1)), FacetIndex((672, 1)), FacetIndex((718, 1)), FacetIndex((720, 1)), FacetIndex((766, 1)), FacetIndex((768, 1))]))), "back" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((147, 3)), FacetIndex((148, 3)), FacetIndex((153, 3)), FacetIndex((154, 3)), FacetIndex((159, 3)), FacetIndex((160, 3)), FacetIndex((165, 3)), FacetIndex((166, 3)), FacetIndex((171, 3)), FacetIndex((172, 3)) … FacetIndex((741, 3)), FacetIndex((742, 3)), FacetIndex((747, 3)), FacetIndex((748, 3)), FacetIndex((753, 3)), FacetIndex((754, 3)), FacetIndex((759, 3)), FacetIndex((760, 3)), FacetIndex((765, 3)), FacetIndex((766, 3))]))), "top" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((581, 3)), FacetIndex((582, 3)), FacetIndex((587, 3)), FacetIndex((588, 3)), FacetIndex((593, 3)), FacetIndex((594, 3)), FacetIndex((599, 3)), FacetIndex((600, 3)), FacetIndex((605, 3)), FacetIndex((606, 3)) … FacetIndex((743, 3)), FacetIndex((744, 3)), FacetIndex((749, 3)), FacetIndex((750, 3)), FacetIndex((755, 3)), FacetIndex((756, 3)), FacetIndex((761, 3)), FacetIndex((762, 3)), FacetIndex((767, 3)), FacetIndex((768, 3))]))), "front" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Tetrahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((2, 1)), FacetIndex((5, 1)), FacetIndex((8, 1)), FacetIndex((11, 1)), FacetIndex((14, 1)), FacetIndex((17, 1)), FacetIndex((20, 1)), FacetIndex((23, 1)), FacetIndex((26, 1)), FacetIndex((29, 1)) … FacetIndex((596, 1)), FacetIndex((599, 1)), FacetIndex((602, 1)), FacetIndex((605, 1)), FacetIndex((608, 1)), FacetIndex((611, 1)), FacetIndex((614, 1)), FacetIndex((617, 1)), FacetIndex((620, 1)), FacetIndex((623, 1))])))), OrderedCollections.OrderedDict{String, Thunderbolt.InterfaceSubdomainDesriptor}()), 225), DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}(SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}[SubDofHandler{DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}}(DofHandler{3, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}}(#= circular reference @-3 =#), OrderedCollections.OrderedSet{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 971, 972, 973, 974, 975, 976, 977, 978, 979, 980]), [:Z], Interpolation[Lagrange{RefHexahedron, 1}()], Int64[], 8)], [:Z], [1, 2, 3, 4, 5, 6, 7, 8, 2, 9 … 1319, 1318, 1139, 1140, 1155, 1154, 1304, 1305, 1320, 1319], [1, 9, 17, 25, 33, 41, 49, 57, 65, 73 … 7761, 7769, 7777, 7785, 7793, 7801, 7809, 7817, 7825, 7833], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, Thunderbolt.SimpleMesh{3, Hexahedron, Float64}(Grid{3, Hexahedron, Float64}(Hexahedron[Hexahedron((1, 2, 17, 16, 166, 167, 182, 181)), Hexahedron((2, 3, 18, 17, 167, 168, 183, 182)), Hexahedron((3, 4, 19, 18, 168, 169, 184, 183)), Hexahedron((4, 5, 20, 19, 169, 170, 185, 184)), Hexahedron((5, 6, 21, 20, 170, 171, 186, 185)), Hexahedron((6, 7, 22, 21, 171, 172, 187, 186)), Hexahedron((7, 8, 23, 22, 172, 173, 188, 187)), Hexahedron((8, 9, 24, 23, 173, 174, 189, 188)), Hexahedron((9, 10, 25, 24, 174, 175, 190, 189)), Hexahedron((10, 11, 26, 25, 175, 176, 191, 190)) … Hexahedron((1130, 1131, 1146, 1145, 1295, 1296, 1311, 1310)), Hexahedron((1131, 1132, 1147, 1146, 1296, 1297, 1312, 1311)), Hexahedron((1132, 1133, 1148, 1147, 1297, 1298, 1313, 1312)), Hexahedron((1133, 1134, 1149, 1148, 1298, 1299, 1314, 1313)), Hexahedron((1134, 1135, 1150, 1149, 1299, 1300, 1315, 1314)), Hexahedron((1135, 1136, 1151, 1150, 1300, 1301, 1316, 1315)), Hexahedron((1136, 1137, 1152, 1151, 1301, 1302, 1317, 1316)), Hexahedron((1137, 1138, 1153, 1152, 1302, 1303, 1318, 1317)), Hexahedron((1138, 1139, 1154, 1153, 1303, 1304, 1319, 1318)), Hexahedron((1139, 1140, 1155, 1154, 1304, 1305, 1320, 1319))], Node{3, Float64}[Node{3, Float64}([0.0, 0.0, 0.0]), Node{3, Float64}([0.5, 0.0, 0.0]), Node{3, Float64}([1.0, 0.0, 0.0]), Node{3, Float64}([1.5, 0.0, 0.0]), Node{3, Float64}([2.0, 0.0, 0.0]), Node{3, Float64}([2.5, 0.0, 0.0]), Node{3, Float64}([3.0, 0.0, 0.0]), Node{3, Float64}([3.5, 0.0, 0.0]), Node{3, Float64}([4.0, 0.0, 0.0]), Node{3, Float64}([4.5, 0.0, 0.0]) … Node{3, Float64}([2.5, 5.0, 3.5]), Node{3, Float64}([3.0, 5.0, 3.5]), Node{3, Float64}([3.5, 5.0, 3.5]), Node{3, Float64}([4.0, 5.0, 3.5]), Node{3, Float64}([4.5, 5.0, 3.5]), Node{3, Float64}([5.0, 5.0, 3.5]), Node{3, Float64}([5.5, 5.0, 3.5]), Node{3, Float64}([6.0, 5.0, 3.5]), Node{3, Float64}([6.5, 5.0, 3.5]), Node{3, Float64}([7.0, 5.0, 3.5])], Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{FacetIndex}}("left" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 5)), FacetIndex((15, 5)), FacetIndex((29, 5)), FacetIndex((43, 5)), FacetIndex((57, 5)), FacetIndex((71, 5)), FacetIndex((85, 5)), FacetIndex((99, 5)), FacetIndex((113, 5)), FacetIndex((127, 5)) … FacetIndex((841, 5)), FacetIndex((855, 5)), FacetIndex((869, 5)), FacetIndex((883, 5)), FacetIndex((897, 5)), FacetIndex((911, 5)), FacetIndex((925, 5)), FacetIndex((939, 5)), FacetIndex((953, 5)), FacetIndex((967, 5))]), "bottom" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((2, 1)), FacetIndex((3, 1)), FacetIndex((4, 1)), FacetIndex((5, 1)), FacetIndex((6, 1)), FacetIndex((7, 1)), FacetIndex((8, 1)), FacetIndex((9, 1)), FacetIndex((10, 1)) … FacetIndex((131, 1)), FacetIndex((132, 1)), FacetIndex((133, 1)), FacetIndex((134, 1)), FacetIndex((135, 1)), FacetIndex((136, 1)), FacetIndex((137, 1)), FacetIndex((138, 1)), FacetIndex((139, 1)), FacetIndex((140, 1))]), "right" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((14, 3)), FacetIndex((28, 3)), FacetIndex((42, 3)), FacetIndex((56, 3)), FacetIndex((70, 3)), FacetIndex((84, 3)), FacetIndex((98, 3)), FacetIndex((112, 3)), FacetIndex((126, 3)), FacetIndex((140, 3)) … FacetIndex((854, 3)), FacetIndex((868, 3)), FacetIndex((882, 3)), FacetIndex((896, 3)), FacetIndex((910, 3)), FacetIndex((924, 3)), FacetIndex((938, 3)), FacetIndex((952, 3)), FacetIndex((966, 3)), FacetIndex((980, 3))]), "back" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((127, 4)), FacetIndex((128, 4)), FacetIndex((129, 4)), FacetIndex((130, 4)), FacetIndex((131, 4)), FacetIndex((132, 4)), FacetIndex((133, 4)), FacetIndex((134, 4)), FacetIndex((135, 4)), FacetIndex((136, 4)) … FacetIndex((971, 4)), FacetIndex((972, 4)), FacetIndex((973, 4)), FacetIndex((974, 4)), FacetIndex((975, 4)), FacetIndex((976, 4)), FacetIndex((977, 4)), FacetIndex((978, 4)), FacetIndex((979, 4)), FacetIndex((980, 4))]), "top" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((841, 6)), FacetIndex((842, 6)), FacetIndex((843, 6)), FacetIndex((844, 6)), FacetIndex((845, 6)), FacetIndex((846, 6)), FacetIndex((847, 6)), FacetIndex((848, 6)), FacetIndex((849, 6)), FacetIndex((850, 6)) … FacetIndex((971, 6)), FacetIndex((972, 6)), FacetIndex((973, 6)), FacetIndex((974, 6)), FacetIndex((975, 6)), FacetIndex((976, 6)), FacetIndex((977, 6)), FacetIndex((978, 6)), FacetIndex((979, 6)), FacetIndex((980, 6))]), "front" => OrderedCollections.OrderedSet{FacetIndex}([FacetIndex((1, 2)), FacetIndex((2, 2)), FacetIndex((3, 2)), FacetIndex((4, 2)), FacetIndex((5, 2)), FacetIndex((6, 2)), FacetIndex((7, 2)), FacetIndex((8, 2)), FacetIndex((9, 2)), FacetIndex((10, 2)) … FacetIndex((845, 2)), FacetIndex((846, 2)), FacetIndex((847, 2)), FacetIndex((848, 2)), FacetIndex((849, 2)), FacetIndex((850, 2)), FacetIndex((851, 2)), FacetIndex((852, 2)), FacetIndex((853, 2)), FacetIndex((854, 2))])), Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()), OrderedCollections.OrderedDict{Tuple{Int64, Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Tuple{Int64, Int64}, Int64}(), OrderedCollections.OrderedDict{Int64, Int64}(), OrderedCollections.OrderedDict{DataType, Int64}(Hexahedron => 980), OrderedCollections.OrderedDict{String, Thunderbolt.VolumetricSubdomainDesriptor}("" => Thunderbolt.VolumetricSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{CellIndex}}(Hexahedron => OrderedCollections.OrderedSet{CellIndex}(CellIndex[CellIndex(1), CellIndex(2), CellIndex(3), CellIndex(4), CellIndex(5), CellIndex(6), CellIndex(7), CellIndex(8), CellIndex(9), CellIndex(10) … CellIndex(971), CellIndex(972), CellIndex(973), CellIndex(974), CellIndex(975), CellIndex(976), CellIndex(977), CellIndex(978), CellIndex(979), CellIndex(980)])))), OrderedCollections.OrderedDict{String, Thunderbolt.SurfaceSubdomainDesriptor}("left" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 5)), FacetIndex((15, 5)), FacetIndex((29, 5)), FacetIndex((43, 5)), FacetIndex((57, 5)), FacetIndex((71, 5)), FacetIndex((85, 5)), FacetIndex((99, 5)), FacetIndex((113, 5)), FacetIndex((127, 5)) … FacetIndex((841, 5)), FacetIndex((855, 5)), FacetIndex((869, 5)), FacetIndex((883, 5)), FacetIndex((897, 5)), FacetIndex((911, 5)), FacetIndex((925, 5)), FacetIndex((939, 5)), FacetIndex((953, 5)), FacetIndex((967, 5))]))), "bottom" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 1)), FacetIndex((2, 1)), FacetIndex((3, 1)), FacetIndex((4, 1)), FacetIndex((5, 1)), FacetIndex((6, 1)), FacetIndex((7, 1)), FacetIndex((8, 1)), FacetIndex((9, 1)), FacetIndex((10, 1)) … FacetIndex((131, 1)), FacetIndex((132, 1)), FacetIndex((133, 1)), FacetIndex((134, 1)), FacetIndex((135, 1)), FacetIndex((136, 1)), FacetIndex((137, 1)), FacetIndex((138, 1)), FacetIndex((139, 1)), FacetIndex((140, 1))]))), "right" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((14, 3)), FacetIndex((28, 3)), FacetIndex((42, 3)), FacetIndex((56, 3)), FacetIndex((70, 3)), FacetIndex((84, 3)), FacetIndex((98, 3)), FacetIndex((112, 3)), FacetIndex((126, 3)), FacetIndex((140, 3)) … FacetIndex((854, 3)), FacetIndex((868, 3)), FacetIndex((882, 3)), FacetIndex((896, 3)), FacetIndex((910, 3)), FacetIndex((924, 3)), FacetIndex((938, 3)), FacetIndex((952, 3)), FacetIndex((966, 3)), FacetIndex((980, 3))]))), "back" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((127, 4)), FacetIndex((128, 4)), FacetIndex((129, 4)), FacetIndex((130, 4)), FacetIndex((131, 4)), FacetIndex((132, 4)), FacetIndex((133, 4)), FacetIndex((134, 4)), FacetIndex((135, 4)), FacetIndex((136, 4)) … FacetIndex((971, 4)), FacetIndex((972, 4)), FacetIndex((973, 4)), FacetIndex((974, 4)), FacetIndex((975, 4)), FacetIndex((976, 4)), FacetIndex((977, 4)), FacetIndex((978, 4)), FacetIndex((979, 4)), FacetIndex((980, 4))]))), "top" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((841, 6)), FacetIndex((842, 6)), FacetIndex((843, 6)), FacetIndex((844, 6)), FacetIndex((845, 6)), FacetIndex((846, 6)), FacetIndex((847, 6)), FacetIndex((848, 6)), FacetIndex((849, 6)), FacetIndex((850, 6)) … FacetIndex((971, 6)), FacetIndex((972, 6)), FacetIndex((973, 6)), FacetIndex((974, 6)), FacetIndex((975, 6)), FacetIndex((976, 6)), FacetIndex((977, 6)), FacetIndex((978, 6)), FacetIndex((979, 6)), FacetIndex((980, 6))]))), "front" => Thunderbolt.SurfaceSubdomainDesriptor(OrderedCollections.OrderedDict{Type, OrderedCollections.OrderedSet{FacetIndex}}(Hexahedron => OrderedCollections.OrderedSet{FacetIndex}(FacetIndex[FacetIndex((1, 2)), FacetIndex((2, 2)), FacetIndex((3, 2)), FacetIndex((4, 2)), FacetIndex((5, 2)), FacetIndex((6, 2)), FacetIndex((7, 2)), FacetIndex((8, 2)), FacetIndex((9, 2)), FacetIndex((10, 2)) … FacetIndex((845, 2)), FacetIndex((846, 2)), FacetIndex((847, 2)), FacetIndex((848, 2)), FacetIndex((849, 2)), FacetIndex((850, 2)), FacetIndex((851, 2)), FacetIndex((852, 2)), FacetIndex((853, 2)), FacetIndex((854, 2))])))), OrderedCollections.OrderedDict{String, Thunderbolt.InterfaceSubdomainDesriptor}()), 1320), [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 1311, 1312, 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320], Dict(1144 => 1144, 1175 => 1175, 719 => 719, 1028 => 1028, 699 => 699, 831 => 831, 1299 => 1299, 1074 => 1074, 319 => 319, 687 => 687…), :φₘ, :Z), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.3228124288530138 0.31328170147136813 … -0.07492340442627253 -0.07689159318238695; 0.042265844904757495 0.044573905360520555 … -0.040074719627719725 -0.038400489638246896], Vector{VertexIndex}[[VertexIndex((281, 5)), VertexIndex((294, 6))], [VertexIndex((287, 6)), VertexIndex((413, 7))]])We compute the ECG online as follows.
io = ParaViewWriter("ep04_ecg")
for (u, t) in TimeChoiceIterator(integrator, tspan[1]:dtvis:tspan[2])
dh = heart_odeform.functions[1].dh
φ = u[heart_odeform.solution_indices[1]]
store_timestep!(io, t, dh.grid) do file
Thunderbolt.store_timestep_field!(file, t, dh, φ, :φₘ)
end
# To compute the ECG we just need to update the ecg cache
Thunderbolt.update_ecg!(geselowitz_ecg, φ)
# which then allows us to evaluate the leads like this
electrode_values = Thunderbolt.evaluate_ecg(geselowitz_ecg)
@info "$t: Lead 1=$(electrode_values[1]) | Lead 2= $(electrode_values[2])"
end[ Info: 0.0: Lead 1=8.662647823456897e-16 | Lead 2= -2.154897274716401e-15
[ Info: 0.5: Lead 1=-1.711156208573285e-5 | Lead 2= 3.193950729256937e-6References
Plain program
Here follows a version of the program without any comments. The file is also available here: ep04_geselowitz-ecg.jl.
using Thunderbolt, LinearAlgebra, StaticArrays, OrdinaryDiffEqOperatorSplitting
function steady_state_initializer!(u₀, f::GenericSplitFunction)
# TODO cleaner implementation. We need to extract this from the types or via dispatch.
heatfun = f.functions[1]
heat_dofrange = f.solution_indices[1]
odefun = f.functions[2]
ionic_model = odefun.ode
φ₀ = @view u₀[heat_dofrange];
# TODO extraction these via utility functions
dh = heatfun.dh
s₀flat = @view u₀[(ndofs(dh)+1):end];
# Should not be reshape but some array of arrays fun
s₀ = reshape(s₀flat, (ndofs(dh), Thunderbolt.num_states(ionic_model)-1));
default_values = Thunderbolt.default_initial_state(ionic_model)
φ₀ .= default_values[1]
for i ∈ 1:(Thunderbolt.num_states(ionic_model)-1)
s₀[:, i] .= default_values[i+1]
end
end
Base.@kwdef struct UniformEndocardialActivation <: Function
transmural_depth::Float64 = 0.15
end
function (p::UniformEndocardialActivation)(x::Vec{3}, t)
τᶠ = 0.25
# TODO source for this
if t ≤ 2.0 && x[1] < p.transmural_depth
return 0.5/τᶠ * exp(t/τᶠ)
else
return 0.0
end
end
protocol = Thunderbolt.AnalyticalTransmembraneStimulationProtocol(
AnalyticalCoefficient(
UniformEndocardialActivation(),
CartesianCoordinateSystem{3}()
),
[SVector((-Inf, Inf))],
)
num_elements_heart = (32,16,16)
num_elements_heart = (8,4,4) # hide
heart_mesh = generate_mesh(Tetrahedron, num_elements_heart, Vec((1.5, 1.5, 0.0)), Vec((5.5, 3.5, 2.0)))
num_elements_torso = (56,40,28)
num_elements_torso = (14,10,7) # hide
torso_mesh = generate_mesh(Hexahedron, num_elements_torso, Vec((0.0, 0.0, 0.0)), Vec((7.0, 5.0, 3.5)))
ground_vertex = Thunderbolt.get_closest_vertex(Vec(0.0, 0.0, 0.0), torso_mesh)
leads = [
[Vec( 0., 0., 1.5), Vec( 7., 0., 1.5)],
[Vec( 3.5, 0., 1.5), Vec( 3.5, 5., 1.5)],
]
microstructure = OrthotropicMicrostructureModel(
ConstantCoefficient((Vec(0.0,0.0,1.0))),
ConstantCoefficient((Vec(0.0,1.0,0.0))),
ConstantCoefficient((Vec(1.0,0.0,0.0))),
)
κ₁ = 0.17 * 0.62 / (0.17 + 0.62)
κᵣ = 0.019 * 0.24 / (0.019 + 0.24)
diffusion_tensor_field = SpectralTensorCoefficient(
microstructure,
ConstantCoefficient(SVector(κ₁, κᵣ, κᵣ))
)
cellmodel = Thunderbolt.PCG2019()
heart_model = MonodomainModel(
ConstantCoefficient(1.0),
ConstantCoefficient(1.0),
diffusion_tensor_field,
protocol,
cellmodel,
:φₘ, :s
)
heart_odeform = semidiscretize(
ReactionDiffusionSplit(heart_model),
FiniteElementDiscretization(Dict(:φₘ => LagrangeCollection{1}())),
heart_mesh,
)
u₀ = zeros(Float64, solution_size(heart_odeform))
steady_state_initializer!(u₀, heart_odeform)
dt₀ = 0.01
dtvis = 0.5
Tₘₐₓ = 50.0
Tₘₐₓ = dtvis # hide
tspan = (0.0, Tₘₐₓ)
problem = OperatorSplittingProblem(heart_odeform, u₀, tspan)
timestepper = LieTrotterGodunov((
BackwardEulerSolver(),
ForwardEulerCellSolver(),
))
integrator = init(problem, timestepper, dt=dt₀, verbose=true)
torso_mesh_κᵢ = ConstantCoefficient(1.0)
torso_mesh_κ = ConstantCoefficient(1.0)
geselowitz_ecg = Thunderbolt.Geselowitz1989ECGLeadCache(
heart_odeform,
torso_mesh,
torso_mesh_κᵢ,
torso_mesh_κ,
leads;
ground = Thunderbolt.OrderedSet([ground_vertex])
)
io = ParaViewWriter("ep04_ecg")
for (u, t) in TimeChoiceIterator(integrator, tspan[1]:dtvis:tspan[2])
dh = heart_odeform.functions[1].dh
φ = u[heart_odeform.solution_indices[1]]
store_timestep!(io, t, dh.grid) do file
Thunderbolt.store_timestep_field!(file, t, dh, φ, :φₘ)
end
# To compute the ECG we just need to update the ecg cache
Thunderbolt.update_ecg!(geselowitz_ecg, φ)
# which then allows us to evaluate the leads like this
electrode_values = Thunderbolt.evaluate_ecg(geselowitz_ecg)
@info "$t: Lead 1=$(electrode_values[1]) | Lead 2= $(electrode_values[2])"
endThis page was generated using Literate.jl.