Modularize microstructure generator (#207)

Author: termi-official

docs/src/api-reference/models.md

@@ -20,8 +20,8 @@ evaluate_coefficient
 ```@docs
 AnisotropicPlanarMicrostructureModel
 OrthotropicMicrostructureModel
-create_simple_microstructure_model
-Thunderbolt.streeter_type_fsn
+create_microstructure_model
+ODB25LTMicrostructureParameters
 ```
 
 ## Boundary Conditions

docs/src/literate-tutorials/cm01_simple-active-stress.jl

@@ -43,11 +43,10 @@ coordinate_system = compute_lv_coordinate_system(mesh);
 
 # In this coordinate system we will now create a microstructure with linearly varying helix angle in transmural direction.
 # The compute microstructure field will be generated on the function space of piecewise continuous first order Lagrange polynomials.
-microstructure = create_simple_microstructure_model(
+microstructure = create_microstructure_model(
     coordinate_system,
-    LagrangeCollection{1}()^3;
-    endo_helix_angle = deg2rad(60.0),
-    epi_helix_angle = deg2rad(-60.0),
+    LagrangeCollection{1}()^3,
+    ODB25LTMicrostructureParameters(),
 );
 
 # Now we describe the model which we want to use.

docs/src/literate-tutorials/cm03_3d0d-coupling.jl

@@ -255,12 +255,11 @@ mesh = Thunderbolt.hexahedralize(mesh)
 #     The 3D0D coupling does not yet support multiple subdomains.
 
 coordinate_system = compute_lv_coordinate_system(mesh)
-microstructure    = create_simple_microstructure_model(
+microstructure    = create_microstructure_model(
     coordinate_system,
     LagrangeCollection{1}()^3,
-    endo_helix_angle = deg2rad(60.0),
-    epi_helix_angle = deg2rad(-60.0),
-)
+    ODB25LTMicrostructureParameters(),
+);
 passive_material_model = Guccione1991PassiveModel()
 active_material_model  = Guccione1993ActiveModel()
 function calcium_profile_function(x::LVCoordinate,t)

src/Thunderbolt.jl

@@ -192,7 +192,8 @@ export
     OrthotropicMicrostructure,
     TransverselyIsotropicMicrostructureModel,
     TransverselyIsotropicMicrostructure,
-    create_simple_microstructure_model,
+    ODB25LTMicrostructureParameters,
+    create_microstructure_model,
     # Coordinate system
     LVCoordinateSystem,
     LVCoordinate,

src/modeling/core/coordinate_systems.jl

@@ -46,7 +46,7 @@ LV only part of the universal ventricular coordinate, containing
     * apicobasal
     * rotational
 """
-struct LVCoordinate{T}
+Base.@kwdef struct LVCoordinate{T}
     transmural::T
     apicobasal::T
     rotational::T
@@ -327,7 +327,7 @@ Biventricular universal coordinate, containing
     * rotational
     * transventricular
 """
-struct BiVCoordinate{T}
+Base.@kwdef struct BiVCoordinate{T}
     transmural::T
     apicobasal::T
     rotational::T

src/modeling/microstructure.jl

@@ -127,72 +127,55 @@ function evaluate_coefficient(fsn::OrthotropicMicrostructureCache, cell_cache, q
     return OrthotropicMicrostructure(orthogonalize_system(f, s, n)...)
 end
 
-
 """
-    streeter_type_fsn(transmural_direction::Vec{3}, circumferential_direction::Vec{3}, apicobasal_direction::Vec{3}, helix_angle, transversal_angle, sheetlet_pseudo_angle, make_orthogonal=true)
-
-Compute fiber, sheetlet and normal direction from the transmural, circumferential, apicobasal directions
-in addition to given helix, transversal and sheetlet angles. The theory is based on the classical work by
-[StreSpoPatRosSon:1969:foc](@citet).
+Linear transmural distribution of the microstructure with three angles to describe all reachable physiological angles.
 """
-function streeter_type_fsn(transmural_direction, circumferential_direction, apicobasal_direction, helix_angle, transversal_angle, sheetlet_pseudo_angle, make_orthogonal=true)
-    # First we construct the helix rotation ...
-    f₀ = rotate_around(circumferential_direction, transmural_direction, helix_angle)
-    f₀ /= norm(f₀)
-    # ... followed by the transversal_angle ...
-    f₀ = rotate_around(f₀, apicobasal_direction, transversal_angle)
-    f₀ /= norm(f₀)
-
-    # Then we construct the the orthogonal sheetlet vector ...
-    s₀ = rotate_around(circumferential_direction, transmural_direction, helix_angle+π/2.0)
-    s₀ /= norm(f₀)
-    # FIXME this does not preserve the sheetlet angle
-    s₀ = unproject(s₀, -transmural_direction, sheetlet_pseudo_angle)
-    if make_orthogonal
-        s₀ = orthogonalize(s₀/norm(s₀), f₀)
-    end
-    # TODO replace above with an implementation of the following pseudocode
-    # 1. Compute plane via P = I - f₀ ⊗ f₀
-    # 2. Eigen decomposition E of P
-    # 3. Compute a generalized eigenvector s' from the non-zero eigenvalue (with ||s'||=1) from E such that <f₀,s'> minimal
-    # 4. Compute s₀ by rotating s' around f₀ such that cos(sheetlet angle) = <s',f₀>
-    s₀ /= norm(s₀)
-
-    # Compute normal :)
-    n₀ = f₀ × s₀
-    n₀ /= norm(n₀)
-
-    return OrthotropicMicrostructure(f₀, s₀, n₀)
+Base.@kwdef struct ODB25LTMicrostructureParameters{T}
+    αendo::T = deg2rad(60.0)
+    αepi::T  = deg2rad(-60.0)
+    βendo::T = 0.0
+    βepi::T  = 0.0
+    γendo::T = 0.0
+    γepi::T  = 0.0
 end
 
-
 """
-    abc_fsn(transmural_direction::Vec{3}, circumferential_direction::Vec{3}, apicobasal_direction::Vec{3}, α, β, γ, make_orthogonal=true)
+    compute_local_microstructure(params::..., x, axes)
 
 Compute fiber, sheetlet and normal direction from the transmural, circumferential, apicobasal directions
 in addition to given helix, transversal and sheetlet angles. The theory is based on the classical work by
 [StreSpoPatRosSon:1969:foc](@citet).
 """
-function abc_fsn(transmural_direction, circumferential_direction, apicobasal_direction, α, β, γ, make_orthogonal=true)
+function compute_local_microstructure(p::ODB25LTMicrostructureParameters, x::Union{LVCoordinate, BiVCoordinate}, axes::Union{LVCoordinate, BiVCoordinate})
+    (; αendo, αepi, βendo, βepi, γendo, γepi) = p
+
+    # linear interpolation of rotation angle
+    α = (1-x.transmural) * αendo + (x.transmural) * αepi
+    β = (1-x.transmural) * βendo + (x.transmural) * βepi
+    γ = (1-x.transmural) * γendo + (x.transmural) * γepi
+
+    circumferential_direction = axes.rotational
+    transmural_direction      = axes.transmural
+    apicobasal_direction      = axes.apicobasal
+
     # First we abc_fsn the helix rotation ...
     f₀ = rotate_around(circumferential_direction, transmural_direction, α)
-    f₀ /= norm(f₀)
+    f₀ = normalize(f₀)
     # ... followed by the transversal_angle ...
-    f₀ = rotate_around(f₀, apicobasal_direction, β)
-    f₀ /= norm(f₀)
+    f₀ = rotate_around(f₀, apicobasal_direction, -β)
+    f₀ = normalize(f₀)
 
     # Then we construct the the orthogonal sheetlet vector ...
     s₀ = rotate_around(circumferential_direction, transmural_direction, α+π/2.0)
-    s₀ /= norm(f₀)
-    if make_orthogonal
-        s₀ = orthogonalize(s₀/norm(s₀), f₀)
-    end
-    s₀ = rotate_around(s₀, f₀, γ)
-    s₀ /= norm(s₀)
+    s₀ = normalize(s₀)
+    s₀ = orthogonalize(s₀, f₀)
+    s₀ = normalize(s₀)
+    s₀ = rotate_around(s₀, f₀, -γ)
+    s₀ = normalize(s₀)
 
     # Compute normal :)
     n₀ = f₀ × s₀
-    n₀ /= norm(n₀)
+    n₀ = normalize(n₀)
 
     return OrthotropicMicrostructure(f₀, s₀, n₀)
 end
@@ -202,7 +185,16 @@ end
 
 Create a rotating fiber field by deducing the circumferential direction from apicobasal and transmural gradients.
 """
-function create_simple_microstructure_model(coordinate_system, ip_collection::VectorizedInterpolationCollection{3}; endo_helix_angle = deg2rad(80.0), epi_helix_angle = deg2rad(-65.0), endo_transversal_angle = 0.0, epi_transversal_angle = 0.0, endo_rot_angle = 0.0, epi_rot_angle = 0.0, make_orthogonal=true)
+function create_simple_microstructure_model(coordinate_system, ip_collection::VectorizedInterpolationCollection{3}; endo_helix_angle = deg2rad(80.0), epi_helix_angle = deg2rad(-65.0), endo_transversal_angle = 0.0, epi_transversal_angle = 0.0, endo_rot_angle = 0.0, epi_rot_angle = 0.0)
+    return create_microstructure_model(coordinate_system, ip_collection, ODB25LTMicrostructureParameters(endo_helix_angle, epi_helix_angle, endo_transversal_angle, epi_transversal_angle, endo_rot_angle, epi_rot_angle))
+end
+
+"""
+    create_microstructure_model(coordinate_system::CoordinateSystemCoefficient, ip::VectorInterpolationCollection, parameters)
+
+Create a rotating fiber field by deducing the circumferential direction from apicobasal and transmural gradients.
+"""
+function create_microstructure_model(coordinate_system::CoordinateSystemCoefficient, ip_collection::VectorizedInterpolationCollection{3}, parameters)
     @unpack dh = coordinate_system
 
     # TODO this storage is redundant, can we reduce the memory footprint?
@@ -234,15 +226,20 @@ function create_simple_microstructure_model(coordinate_system, ip_collection::Ve
                 apicobasal_direction = apicobasal_direction - (apicobasal_direction ⋅ transmural_direction) * transmural_direction # We do this fix to ensure local orthogonality
                 circumferential_direction = transmural_direction × apicobasal_direction
                 circumferential_direction /= norm(circumferential_direction)
+                axes = LVCoordinate(;
+                    transmural = transmural_direction,
+                    rotational = circumferential_direction,
+                    apicobasal = apicobasal_direction,
+                )
+                # TODO grab these via some interface!
+                x = LVCoordinate(;
+                    transmural = function_value(cv, qp, coordinate_system.u_transmural[dof_indices]),
+                    rotational = function_value(cv, qp, coordinate_system.u_rotational[dof_indices]),
+                    apicobasal = function_value(cv, qp, coordinate_system.u_apicobasal[dof_indices]),
+                )
 
-                transmural  = function_value(cv, qp, coordinate_system.u_transmural[dof_indices])
-
-                # linear interpolation of rotation angle
-                helix_angle       = (1-transmural) * endo_helix_angle + (transmural) * epi_helix_angle
-                transversal_angle = (1-transmural) * endo_transversal_angle + (transmural) * epi_transversal_angle
-                rot_angle         = (1-transmural) * endo_rot_angle + (transmural) * epi_rot_angle
+                coeff = compute_local_microstructure(parameters, x, axes)
 
-                coeff = abc_fsn(transmural_direction, circumferential_direction, apicobasal_direction, helix_angle, transversal_angle, rot_angle, make_orthogonal)
                 f_buf[qp.i, cellindex] = Tv(coeff.f)
                 s_buf[qp.i, cellindex] = Tv(coeff.s)
                 n_buf[qp.i, cellindex] = Tv(coeff.n)

test/integration/test_solid_mechanics.jl

@@ -221,7 +221,8 @@ end
     @test !any(isnan.(cs.u_apicobasal))
     @test !any(isnan.(cs.u_transmural))
     @test !any(isnan.(cs.u_rotational))
-    microstructure_model = create_simple_microstructure_model(cs, LagrangeCollection{1}()^3)
+    microstructure_parameters = ODB25LTMicrostructureParameters(αendo=deg2rad(80.0), αepi=deg2rad(-65.0))
+    microstructure_model      = create_microstructure_model(cs, LagrangeCollection{1}()^3, microstructure_parameters)
 
     test_solve_contractile_ideal_lv(grid, ExtendedHillModel(
         HolzapfelOgden2009Model(),

test/test_microstructures.jl

@@ -29,13 +29,15 @@
     end
 
     @testset "OrthotropicMicrostructureModel" begin
-        ms = create_simple_microstructure_model(ring_cs, ip_collection,
-            endo_helix_angle = deg2rad(0.0),
-            epi_helix_angle = deg2rad(0.0),
-            endo_transversal_angle = 0.0,
-            epi_transversal_angle = 0.0,
-            endo_rot_angle = deg2rad(0.0),
-            epi_rot_angle = deg2rad(0.0),
+        ms = create_microstructure_model(ring_cs, ip_collection, 
+            ODB25LTMicrostructureParameters(
+                αendo = 0.0,
+                αepi  = 0.0,
+                βendo = 0.0,
+                βepi  = 0.0,
+                γendo = 0.0,
+                γepi  = 0.0,
+            )
         )
         cache2 = Thunderbolt.setup_coefficient_cache(ms, qr, sdh)
         for cellcache in CellIterator(ring_grid)