Propagate coefficient API into spiral wave example.

examples/common-stuff.jl

@@ -40,24 +40,6 @@ import Thunderbolt: Ψ, U
 
 # ----------------------------------------------
 
-mutable struct LazyMicrostructureCache{MM, VT}
-    const microstructure_model::MM
-    const x_ref::Vector{VT}
-    cellid::Int
-end
-
-function Thunderbolt.directions(cache::LazyMicrostructureCache{MM}, qp::Int) where {MM}
-    return directions(cache.microstructure_model, cache.cellid, cache.x_ref[qp])
-end
-
-function setup_microstructure_cache(cv, model::OrthotropicMicrostructureModel{FiberCoefficientType, SheetletCoefficientType, NormalCoefficientType}) where {FiberCoefficientType, SheetletCoefficientType, NormalCoefficientType}
-    return LazyMicrostructureCache(model, cv.qr.points, -1)
-end
-
-function update_microstructure_cache!(cache::LazyMicrostructureCache{MM}, time::Float64, cell::CellCacheType, cv::CV) where {CellCacheType, CV, MM}
-    cache.cellid = cellid(cell)
-end
-
 
 # ----------------------------------------------
 
@@ -77,7 +59,7 @@ end
 function update_contraction_model_cache!(cache::PelceSunLangeveld1995Cache, time::Float64, cell::CellCacheType, cv::CV, calcium_field::CF) where {CellCacheType, CV, CF}
     for qp ∈ 1:getnquadpoints(cv)
         x_ref = cv.qr.points[qp]
-        cache.calcium_values[qp] = value(calcium_field, Ferrite.cellid(cell), x_ref, time)
+        cache.calcium_values[qp] = evaluate_coefficient(calcium_field, Ferrite.cellid(cell), x_ref, time)
     end
 end
 
@@ -454,7 +436,7 @@ struct CalciumHatField end
 
 """
 """
-value(coeff::CalciumHatField, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim} = t < 1.0 ? t : 2.0-t
+Thunderbolt.evaluate_coefficient(coeff::CalciumHatField, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim} = t < 1.0 ? t : 2.0-t
 
 """
 Parameterization from Vallespin paper.

examples/ring-refactored.jl

@@ -219,7 +219,8 @@ function (postproc::StandardMechanicalIOPostProcessor{IO, CV, MC})(problem, uₜ
     end
 
     # Save the solution
-    Thunderbolt.store_timestep!(io, t, dh, uₜ)
+    Thunderbolt.store_timestep!(io, t, dh)
+    Thunderbolt.store_timestep_field!(io, t, dh, uₜ, :displacement)
     Thunderbolt.store_timestep_celldata!(io, t, hcat(frefdata...),"Reference Fiber Data")
     Thunderbolt.store_timestep_celldata!(io, t, hcat(fdata...),"Current Fiber Data")
     Thunderbolt.store_timestep_celldata!(io, t, hcat(srefdata...),"Reference Sheet Data")

examples/spiral-wave.jl

@@ -110,20 +110,44 @@ function setup_mass_operator!(bifi::AssembledMassOperator{MT, CV}, dh::DH) where
     end
 end
 
+using StaticArrays
+struct PlanarDiffusionTensorCoefficient{MSC}
+    microstructure_cache::MSC
+    conductivities::SVector{2}
+end
+
+function Thunderbolt.evaluate_coefficient(coeff::PlanarDiffusionTensorCoefficient{MSC}, cell_id::Int, ξ::Vec{rdim}, t::Float64=0.0) where {MSC, rdim}
+    f₀, s₀ = directions(coeff.microstructure_cache, cell_id, ξ, t)
+    return coeff.conductivities[1] * f₀ ⊗ f₀ + coeff.conductivities[2] * s₀ ⊗ s₀
+end
+
+struct ConductivityToDiffusivityCoefficient{DTC, CC, STVC}
+    conductivity_tensor_coefficient::DTC
+    capacitance_coefficient::CC
+    χ_coefficient::STVC
+end
+
+function Thunderbolt.evaluate_coefficient(coeff::ConductivityToDiffusivityCoefficient{DTC, CC, STVC}, cell_id::Int, ξ::Vec{rdim}, t::Float64=0.0) where {DTC, CC, STVC, rdim}
+    κ = evaluate_coefficient(coeff.conductivity_tensor_coefficient, cell_id, ξ)
+    Cₘ = evaluate_coefficient(coeff.capacitance_coefficient, cell_id, ξ)
+    χ = evaluate_coefficient(coeff.χ_coefficient, cell_id, ξ)
+    return κ/(Cₘ*χ)
+end
+
 @doc raw"""
     AssembledDiffusionOperator{MT, DT, CV}
 
 Assembles the matrix associated to the bilinearform ``a(u,v) = -\int \nabla v(x) \cdot D(x) \nabla u(x) dx`` for a given diffusion tensor ``D(x)`` and ``u,v`` from the same function space.
 """
 struct AssembledDiffusionOperator{MT, DT, CV}
     K::MT
-    coefficient::DT
+    diffusion_coefficient::DT
     cv::CV
 end
 mul!(b, K::AssembledDiffusionOperator{MT, DT, CV}, uₙ₋₁) where {MT, DT, CV} = mul!(b, K.K, uₙ₋₁)
 
 function setup_diffusion_operator!(bifi::AssembledDiffusionOperator{MT, CV}, dh::DH) where {MT, CV, DH}
-    @unpack K, coefficient, cv = bifi
+    @unpack K, diffusion_coefficient, cv = bifi
 
     n_basefuncs = getnbasefunctions(cv)
     Kₑ = zeros(n_basefuncs, n_basefuncs)
@@ -137,16 +161,16 @@ function setup_diffusion_operator!(bifi::AssembledDiffusionOperator{MT, CV}, dh:
         reinit!(cv, cell)
 
         for q_point in 1:getnquadpoints(cv)
-            # TODO coefficient via AssembledDiffusionOperator::coefficient
-            D_loc = 0.001 # evaluate(coefficient, q_point)
+            ξ = cv.qr.points[q_point]
+            D_loc = evaluate_coefficient(diffusion_coefficient, cellid(cell), ξ)
             #based on the gauss point coordinates, we get the spatial dependent
             #material parameters
             dΩ = getdetJdV(cv, q_point)
             for i in 1:n_basefuncs
                 ∇Nᵢ = shape_gradient(cv, q_point, i)
                 for j in 1:n_basefuncs
                     ∇Nⱼ = shape_gradient(cv, q_point, j)
-                    Kₑ[i,j] -= ((D_loc * ∇Nᵢ) ⋅ ∇Nⱼ) * dΩ
+                    Kₑ[i,j] -= ((D_loc ⋅ ∇Nᵢ) ⋅ ∇Nⱼ) * dΩ
                 end
             end
         end
@@ -203,7 +227,7 @@ function setup_solver_caches(problem #= ::TransientHeatProblem =#, solver::Impli
         ),
         AssembledDiffusionOperator(
             create_sparsity_pattern(dh),
-            nothing,
+            problem.diffusion_tensor_field,
             cv
         ),
         create_sparsity_pattern(dh),
@@ -510,7 +534,7 @@ function semidiscretize(split::ReactionDiffusionSplit{MonodomainModel{A,B,C,D,E}
     #
     semidiscrete_problem = SplitProblem(
         TransientHeatProblem(
-            x -> epmodel.κ(x)/(epmodel.Cₘ(x)*epmodel.χ),
+            ConductivityToDiffusivityCoefficient(epmodel.κ, epmodel.Cₘ, epmodel.χ),
             epmodel.stim,
             dh
         ),
@@ -549,9 +573,9 @@ function spiral_wave_initializer(problem::SplitProblem)
 end
 
 model = MonodomainModel(
-    x->1.0,
-    x->1.0,
-    x->SymmetricTensor{2,2,Float64}((4.5e-5, 0, 2.0e-5)), # TODO coefficient API
+    ConstantCoefficient(1.0),
+    ConstantCoefficient(1.0),
+    ConstantCoefficient(SymmetricTensor{2,2,Float64}((4.5e-5, 0, 2.0e-5))),
     NoStimulationProtocol(),
     FHNModel()
 )

src/exports.jl

@@ -31,9 +31,14 @@ export
     NoStimulationProtocol,
     # Microstructure
     OrthotropicMicrostructureModel,
-    ConstantFieldCoefficient,
-    create_simple_fiber_model,
     directions,
+    FieldCoefficient,
+    ConstantCoefficient,
+    evaluate_coefficient,
+    create_simple_fiber_model,
+    update_microstructure_cache!,
+    setup_microstructure_cache,
+    LazyCoefficientCache,
     # Coordinate system
     LVCoordinateSystem,
     CartesianCoordinateSystem,

src/modeling/microstructure.jl

@@ -7,7 +7,7 @@ end
 
 """
 """
-function value(coeff::FieldCoefficient{TA,IP}, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim,TA,IP}
+function evaluate_coefficient(coeff::FieldCoefficient{TA,IP}, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim,TA,IP}
     @unpack elementwise_data, ip = coeff
 
     n_base_funcs = Ferrite.getnbasefunctions(ip)
@@ -21,15 +21,26 @@ end
 
 """
 """
-struct ConstantFieldCoefficient{T}
+struct ConstantCoefficient{T}
     val::T
 end
 
 """
 """
-value(coeff::ConstantFieldCoefficient{T}, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim,T} = coeff.val
+evaluate_coefficient(coeff::ConstantCoefficient{T}, cell_id::Int, ξ::Vec{dim}, t::Float64=0.0) where {dim,T} = coeff.val
 
 
+struct AnisotropicPlanarMicrostructureModel{FiberCoefficientType, SheetletCoefficientType}
+    fiber_coefficient::FiberCoefficientType
+    sheetlet_coefficient::SheetletCoefficientType
+end
+
+function directions(fsn::AnisotropicPlanarMicrostructureModel, cell_id::Int, ξ::Vec{dim}, t = 0.0) where {dim}
+    f₀ = evaluate_coefficient(fsn.fiber_coefficient, cell_id, ξ, t)
+    s₀ = evaluate_coefficient(fsn.sheetlet_coefficient, cell_id, ξ, t)
+
+    f₀, s₀
+end
 
 struct OrthotropicMicrostructureModel{FiberCoefficientType, SheetletCoefficientType, NormalCoefficientType}
     fiber_coefficient::FiberCoefficientType
@@ -38,14 +49,13 @@ struct OrthotropicMicrostructureModel{FiberCoefficientType, SheetletCoefficientT
 end
 
 function directions(fsn::OrthotropicMicrostructureModel, cell_id::Int, ξ::Vec{dim}, t = 0.0) where {dim}
-    f₀ = value(fsn.fiber_coefficient, cell_id, ξ, t)
-    s₀ = value(fsn.sheetlet_coefficient, cell_id, ξ, t)
-    n₀ = value(fsn.normal_coefficient, cell_id, ξ, t)
+    f₀ = evaluate_coefficient(fsn.fiber_coefficient, cell_id, ξ, t)
+    s₀ = evaluate_coefficient(fsn.sheetlet_coefficient, cell_id, ξ, t)
+    n₀ = evaluate_coefficient(fsn.normal_coefficient, cell_id, ξ, t)
 
     f₀, s₀, n₀
 end
 
-
 """
 """
 function generate_nodal_quadrature_rule(ip::Interpolation{ref_shape, order}) where {ref_shape, order}
@@ -131,3 +141,26 @@ function create_simple_fiber_model(coordinate_system, ip_fiber::ScalarInterpolat
         FieldCoefficient(elementwise_data_n, ip_fiber)
     )     
 end
+
+# TODO where to move this? Technically this is assembly infrastructure
+mutable struct LazyMicrostructureCache{MM, VT}
+    const microstructure_model::MM
+    const x_ref::Vector{VT}
+    cellid::Int
+end
+
+function directions(cache::LazyMicrostructureCache{MM}, qp::Int) where {MM}
+    return directions(cache.microstructure_model, cache.cellid, cache.x_ref[qp])
+end
+
+function setup_microstructure_cache(cv, model::AnisotropicPlanarMicrostructureModel{FiberCoefficientType, SheetletCoefficientType}) where {FiberCoefficientType, SheetletCoefficientType}
+    return LazyMicrostructureCache(model, cv.qr.points, -1)
+end
+
+function setup_microstructure_cache(cv, model::OrthotropicMicrostructureModel{FiberCoefficientType, SheetletCoefficientType, NormalCoefficientType}) where {FiberCoefficientType, SheetletCoefficientType, NormalCoefficientType}
+    return LazyMicrostructureCache(model, cv.qr.points, -1)
+end
+
+function update_microstructure_cache!(cache::LazyMicrostructureCache{MM}, time::Float64, cell::CellCacheType, cv::CV) where {CellCacheType, CV, MM}
+    cache.cellid = cellid(cell)
+end