[WIP] Reaction-Eikonal example by AbdAlazezAhmed · Pull Request #241 · JuliaHealth/Thunderbolt.jl

AbdAlazezAhmed · 2026-01-29T20:20:10Z

No description provided.

…al LV mesh

codecov · 2026-01-29T20:34:11Z

Codecov Report

❌ Patch coverage is 62.13992% with 92 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.59%. Comparing base (552693a) to head (ae6209d).

Files with missing lines	Patch %	Lines
src/mesh/tools.jl	79.19%	36 Missing ⚠️
src/discretization/eikonal.jl	0.00%	22 Missing ⚠️
src/modeling/electrophysiology.jl	5.55%	17 Missing ⚠️
src/discretization/fem.jl	31.81%	15 Missing ⚠️
src/modeling/functions.jl	75.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #241      +/-   ##
==========================================
- Coverage   73.03%   72.59%   -0.44%     
==========================================
  Files          78       79       +1     
  Lines        6865     7018     +153     
==========================================
+ Hits         5014     5095      +81     
- Misses       1851     1923      +72

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

termi-official · 2026-03-06T16:45:42Z

What is the state of this PR?
Also, where is the Eikonal package? And we might want to add some test beyond the integration test in the docs.

AbdAlazezAhmed · 2026-03-06T17:18:06Z

What is the state of this PR?

I'll write some tests and check documentation. Already working fine in another project using this branch.

Also, where is the Eikonal package? And we might want to add some test beyond the integration test in the docs.

It's in the noc git instance. Will make it public and link to it after adding docs

…olt.jl into ahm/eikonal

This reverts commit 51e4a08.

This reverts commit 81e0d2f.

termi-official

There are literally hundreds of random ass changes related to code style. Can you revert these please?

termi-official · 2026-03-12T21:47:31Z

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"


Why? I explicitly try to avoid pulling all solvers under the sun in.

src/mesh/tools.jl

src/discretization/fem.jl

termi-official · 2026-03-12T21:50:16Z

src/discretization/eikonal.jl

+            end
+        end
+    end
+    cells = getproperty.(cells_ferrite, :nodes)


What is this?

In the FastIterativeMethod.jl I didn't want to depend on Ferrite's definition of a cell and used just a tuple of ints (which is what Ferrite cell wraps)

termi-official · 2026-03-12T21:50:42Z

src/discretization/eikonal.jl

+end
+
+function semidiscretize(
+    ::EikonalModel,


Don't we need to pass the coefficient?

For what specifically? if you mean diffusivity I think we can take it from the inner MonodomainModel for consistency (or remove it). For now it's provided at solve

termi-official · 2026-03-12T21:51:19Z

src/discretization/eikonal.jl

+
+function semidiscretize(
+    ::EikonalModel,
+    discretization::SimplicialEikonalDiscretization,


Shouldn't this be FastIterativeMethod or so? Also, can't this be an extension?

I'm not sure, since here we define the discretization not the -solver?- method. It's just a current constraint that the FIM implementaiton we have now is limited to linear tets

termi-official · 2026-03-12T21:54:12Z

src/discretization/eikonal.jl

+"""
+Coordinates used for defining activation points when solving the eikonal equation.
+"""
+struct ActivationCoordinate{CoordT <: Union{LVCoordinate, BiVCoordinate}, T}
+    coord::CoordT
+    time_offset::T
+end
+
+abstract type AbstractEikonalActivationProtocol end
+
+"""
+Activation protocol for activating specific points with time offsets to mimic the behavior of a Purkinje network.
+"""
+struct NodalEikonalActivationProtocol{VectorT <: AbstractVector{<:ActivationCoordinate}} <:
+       AbstractEikonalActivationProtocol
+    nodes::VectorT
+end
+
+"""
+Activation protocol for uniformally activating the endocardium with zero wave arrival time.
+"""
+struct UniformEndocardialEikonalActivationProtocol <: AbstractEikonalActivationProtocol end
+
+"""
+Activation protocol for uniformally activating the endocardium with zero wave arrival time.
+"""
+struct AnalyticalEikonalActivationProtocol{fT} <: AbstractEikonalActivationProtocol
+    f::fT
+end
+


This isn't discretization info, but model, right?.

termi-official · 2026-03-12T21:54:59Z

src/discretization/eikonal.jl

+end
+
+function get_nodes(
+    ::UniformEndocardialEikonalActivationProtocol,


Why is this special to Eikonal?

termi-official · 2026-03-12T21:56:08Z

src/discretization/fem.jl

+A wrapper function containing the cell model dynamics function, as a pointwise ODE,
+and the eiknal function that needs to be solved once before solving the cell dynamics part.
+"""
+struct EikonalCoupledODEFunction{ODEFunctionT, EikonalFunctionT}


Isn't this more like ReactionEikonalFunction (to properly attribute the origin)?

termi-official · 2026-03-12T21:57:27Z

src/mesh/tools.jl

    enids = new_edge_offset .+ global_edge_indices
    fnids = new_face_offset .+ global_face_indices
-    cnid  = new_face_offset + num_faces(mgrid) + cell_idx
+    cnid = new_face_offset + num_faces(mgrid) + cell_idx


Can you please revert all the noise here? It makes really hard to track down what the actual relevant changes are...

Ah I thought last commits reverted them. Will do.

termi-official · 2026-03-13T15:05:37Z

src/mesh/tools.jl

+function tetrahedralize(grid::Grid{3, Hexahedron})
+    return _tetrahedralize(to_mesh(grid))
+end
+
+function tetrahedralize(mesh::SimpleMesh{3, Hexahedron})
+    grid = _tetrahedralize(mesh)
+    return to_mesh(grid)
+end


Why restrict to Hex. We allow Tet below, right?

termi-official · 2026-03-13T15:05:55Z

src/mesh/tools.jl

+    for (cell_idx, cell) in enumerate(cells)
+        push!(cell_offsets, length(new_cells))
+        start_idx = length(new_cells) + 1
+        if isa(cell, Hexahedron)


Use dispatch

termi-official · 2026-03-13T17:49:06Z

src/modeling/electrophysiology.jl


+"""
+    ReactionEikonalSplit(model)
+    ReactionEikonalSplit(model, coeff)


Suggested change

ReactionEikonalSplit(model, coeff)

ReactionEikonalSplit(model, cs)

termi-official · 2026-03-13T17:49:38Z

src/modeling/electrophysiology.jl

+Annotation for the reaction-Eikonal split of a given model. The
+second argument is a coefficient describing the input `x` for the reaction model rhs,


Suggested change

Annotation for the reaction-Eikonal split of a given model. The

second argument is a coefficient describing the input `x` for the reaction model rhs,

Annotation for the reaction-Eikonal split of a given model. The

argument `cs` is a coefficient describing the input `x` for the reaction model rhs,

termi-official · 2026-03-13T17:51:03Z

src/modeling/electrophysiology.jl

+        end
+    else
+        if m.stim_offset ≤ t ≤ m.stim_offset + m.stim_length
+            τᶠ = 0.25


This is an input parameter, right?

termi-official · 2026-03-13T18:10:03Z

src/mesh/tools.jl

+                t1 = orient_to_positive((v1, v3, v2, cnid), new_nodes)
+                t2 = orient_to_positive((v1, v4, v3, cnid), new_nodes)
+                push!(new_cells, Tetrahedron(t1))
+                push!(new_cells, Tetrahedron(t2))


What kind of algorithm is this and what is the insertion order (in terms of space-filling curves)?

termi-official · 2026-03-13T18:13:43Z

src/mesh/tools.jl

+            triangles =
+                length(fnodes) == 3 ? [(fnodes[1], fnodes[2], fnodes[3])] :
+                [(fnodes[1], fnodes[3], fnodes[2]), (fnodes[1], fnodes[4], fnodes[3])]
+            rng = cell_ranges[cellidx]
+            for tri in triangles
+                tri_set = Set(tri)
+                matched = false
+                for newcid in rng
+                    for localf = 1:length(Ferrite.faces(new_cells[newcid]))
+                        if Set(Ferrite.faces(new_cells[newcid])[localf]) == tri_set
+                            push!(s, FacetIndex(newcid, localf))
+                            matched = true
+                            break
+                        end
+                    end
+                    matched && break
+                end
+                if !matched
+                    throw(
+                        error(
+                            "Could not map facet (cell=$cellidx, face=$lfi, tri=$tri) into tetrahedral mesh",
+                        ),
+                    )
+                end
+            end


This is an overengineered solution. The mapping is predetermined if you remove the orient_to_positive function above and replace it with a deterministic insertion which follows a given convention. There really should be a single loop up table for this which should be always true. This kind of refactor is something Claude/Codex/whatever (or DIY) can do locally (if you give them also access to Ferrite an and LSP). Also, use dispatches.

As a test case you can use a single element (one Hex, one Tet) and mark all faces individually with different names. Then you can test if the face center of the new faces is where you expect them.

termi-official · 2026-03-13T18:16:38Z

src/mesh/tools.jl

+            center = zero(new_nodes[vnids[1]].x)
+            for vid in vnids
+                center += new_nodes[vid].x
+            end
+            center /= length(vnids)


Not necessarily true in general. There should be a function to compute the center which we should use to compute the center with the ansatz functions (which probably right now implements the same approximation scheme).

termi-official · 2026-03-13T18:17:25Z

src/mesh/tools.jl

+end
+
+function _tetrahedralize(mgrid::SimpleMesh{3, <:Any, T}) where {T}
+    materialize_all_entities!(mgrid)


Why all and where do you use them?

termi-official · 2026-03-13T18:18:14Z

src/mesh/tools.jl

+
+    new_nodes = copy(grid.nodes)
+    new_cells = Tetrahedron[]
+    cell_ranges = UnitRange{Int}[]


Is there some advantage of a UnitRange over just an Int?

AbdAlazezAhmed added 7 commits January 20, 2026 14:39

homeoffice, lots of copy paste

a116331

tetrahedralize (largely AI generated slop), needed for Eikonal tutori…

5742e04

…al LV mesh

initial eikonal split

a3670f7

oopsie

aaeeee7

delete commented code

7dbf020

docs

40c414d

add tests

35c8411

AbdAlazezAhmed added 12 commits March 10, 2026 20:11

Add tests for Tetrahedralize

22dd391

homeoffice, writing tests

75f09ad

apply runic partially

1a236b4

Merge branch 'main' into ahm/eikonal

ae6209d

tests

e54196c

formatting

81e0d2f

formatting 2

51e4a08

apply formatting

fc057ff

Merge branch 'ahm/eikonal' of https://github.com/JuliaHealth/Thunderb…

09994c8

…olt.jl into ahm/eikonal

Revert "formatting 2"

3018d75

This reverts commit 51e4a08.

Revert "formatting"

99903ab

This reverts commit 81e0d2f.

FORMAT

f596977

termi-official requested changes Mar 12, 2026

View reviewed changes

AbdAlazezAhmed added 2 commits March 13, 2026 15:21

undo the not undoed formatting

bd505a9

ignore Eikonal CI for 1.10

5a45f96

termi-official requested changes Mar 13, 2026

View reviewed changes

		@@ -1,7 +1,9 @@
		[deps]
		DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"

	ReactionEikonalSplit(model, coeff)
	ReactionEikonalSplit(model, cs)

		Annotation for the reaction-Eikonal split of a given model. The
		second argument is a coefficient describing the input `x` for the reaction model rhs,

Conversation

AbdAlazezAhmed commented Jan 29, 2026

Uh oh!

codecov bot commented Jan 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

termi-official commented Mar 6, 2026

Uh oh!

AbdAlazezAhmed commented Mar 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

termi-official left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

codecov bot commented Jan 29, 2026 •

edited

Loading

AbdAlazezAhmed commented Mar 6, 2026 •

edited

Loading