Using collections in user-facing API

benchmarks/benchmarks-coefficients.jl

@@ -0,0 +1,64 @@
+using Thunderbolt, BenchmarkTools, StaticArrays
+
+cell_cache = Ferrite.CellCache(generate_grid(Line, (2,)))
+qp1 = QuadraturePoint(1, Vec((0.0,)))
+ip_collection = LagrangeCollection{1}()
+reinit!(cell_cache, 1)
+
+for val ∈ [1.0, one(Tensor{2,2})]
+    cc = ConstantCoefficient(val)
+    @btime evaluate_coefficient($cc, $cell_cache, $qp1, 0.0)
+end
+
+data_scalar = zeros(2,2,1)
+data_scalar[1,1] =  1.0
+data_scalar[1,2] = -1.0
+data_scalar[2,1] = -1.0
+fcs = FieldCoefficient(data_scalar, ip_collection)
+@btime evaluate_coefficient($fcs, $cell_cache, $qp1, 0.0)
+
+data_vector = zeros(Vec{2,Float64},2,2)
+data_vector[1,1] = Vec((1.0,0.0))
+data_vector[1,2] = Vec((-1.0,-0.0))
+data_vector[2,1] = Vec((0.0,-1.0))
+fcv = FieldCoefficient(data_vector, ip_collection^2)
+@btime evaluate_coefficient($fcv, $cell_cache, $qp1, 0.0)
+
+ccsc = CartesianCoordinateSystemCoefficient(ip_collection^1)
+@btime evaluate_coefficient($ccsc, $cell_cache, $qp1, 0.0)
+
+ac = AnalyticalCoefficient(
+    (x,t) -> norm(x)+t,
+    CartesianCoordinateSystemCoefficient(ip_collection^1)
+    )
+@btime evaluate_coefficient($ac, $cell_cache, $qp1, 0.0)
+
+eigvec = Vec((1.0,0.0))
+eigval = -1.0
+stc = SpectralTensorCoefficient(
+    ConstantCoefficient(SVector((eigvec,))),
+    ConstantCoefficient(SVector((eigval,))),
+    )
+st = Tensor{2,2}((-1.0,0.0,0.0,0.0))
+@btime evaluate_coefficient($stc, $cell_cache, $qp1, 0.0)
+
+shdc = SpatiallyHomogeneousDataField(
+    [1.0, 2.0],
+    [Vec((0.1,)), Vec((0.2,)), Vec((0.3,))]
+    )
+@btime evaluate_coefficient($shdc, $cell_cache, $qp1, 0.0)
+
+eigvec = Vec((1.0,0.0))
+eigval = -1.0
+stc = SpectralTensorCoefficient(
+    ConstantCoefficient(SVector((eigvec,))),
+    ConstantCoefficient(SVector((eigval,))),
+    )
+st = Tensor{2,2}((-1.0,0.0,0.0,0.0))
+ctdc = Thunderbolt.ConductivityToDiffusivityCoefficient(
+    stc,
+    ConstantCoefficient(2.0),
+    ConstantCoefficient(0.5),
+    )
+
+@btime evaluate_coefficient($stc, $cell_cache, $qp1, 0.0)

benchmarks/benchmarks-linear-form.jl

@@ -0,0 +1,20 @@
+using BenchmarkTools, Thunderbolt, StaticArrays
+celltype = Hexahedron
+cell_cache = Ferrite.CellCache(generate_grid(celltype, (1,1,1)))
+reinit!(cell_cache,1)
+ref_shape = RefHexahedron
+order = 1
+ip_collection = LagrangeCollection{order}()
+ip = getinterpolation(ip_collection, ref_shape)
+qr_collection = QuadratureRuleCollection(2)
+qr = getquadraturerule(celltype)
+cv = CellValues(qr, ip)
+ac = AnalyticalCoefficient(
+    (x,t) -> norm(x)+t,
+    CartesianCoordinateSystemCoefficient(
+        getinterpolation(LagrangeCollection{1}()^3)
+    )
+)
+b = zeros(8)
+element_cache = Thunderbolt.AnalyticalCoefficientElementCache(ac, [SVector((0.0,1.0))], cv)
+@btime Thunderbolt.assemble_element!($b, $cell_cache, $element_cache, 0.0)

examples/LV.jl

@@ -150,7 +150,7 @@ function (postproc::StandardMechanicalIOPostProcessor2)(t, problem, solver_cache
 end
 
 
-function solve_ideal_lv(name_base, constitutive_model, grid, coordinate_system, face_models, ip_mech::IPM, ip_geo::IPG, intorder::Int, Δt = 100.0, T = 1000.0) where {ref_shape, IPM <: Interpolation{ref_shape}, IPG <: Interpolation{ref_shape}}
+function solve_ideal_lv(name_base, constitutive_model, grid, coordinate_system, face_models, ip_collection::Thunderbolt.ScalarInterpolationCollection, ip_mech::IPM, ip_geo::IPG, qr_collection::QuadratureRuleCollection, Δt = 100.0, T = 1000.0) where {ref_shape, IPM <: Interpolation{ref_shape}, IPG <: Interpolation{ref_shape}}
     io = ParaViewWriter(name_base);
     # io = JLD2Writer(name_base);
 
@@ -167,7 +167,7 @@ function solve_ideal_lv(name_base, constitutive_model, grid, coordinate_system,
             ]
         )),
         FiniteElementDiscretization(
-            Dict(:displacement => LagrangeCollection{1}()^3),
+            Dict(:displacement => ip_collection^3),
             Dirichlet[],
             # [Dirichlet(:displacement, getfaceset(grid, "Myocardium"), (x,t) -> [0.0], [3])],
             # [Dirichlet(:displacement, getfaceset(grid, "left"), (x,t) -> [0.0], [1]),Dirichlet(:displacement, getfaceset(grid, "front"), (x,t) -> [0.0], [2]),Dirichlet(:displacement, getfaceset(grid, "bottom"), (x,t) -> [0.0], [3]), Dirichlet(:displacement, Set([1]), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3])],
@@ -176,7 +176,8 @@ function solve_ideal_lv(name_base, constitutive_model, grid, coordinate_system,
     )
 
     # Postprocessor
-    cv_post = CellValues(QuadratureRule{ref_shape}(intorder-1), ip_mech, ip_geo)
+    qr = getquadraturerule(qr_collection, elementtypes(grid)[1])
+    cv_post = CellValues(qr, ip_mech, ip_geo)
     standard_postproc = StandardMechanicalIOPostProcessor2(io, cv_post, [1], coordinate_system)
 
     # Create sparse matrix and residual vector
@@ -198,11 +199,14 @@ end
 LV_grid = Thunderbolt.hexahedralize(Thunderbolt.generate_ideal_lv_mesh(15,2,6))
 ref_shape = RefHexahedron
 order = 1
-ip = Lagrange{ref_shape, order}()^3
-ip_fiber = Lagrange{ref_shape, order}()
-ip_geo = Lagrange{ref_shape, order}()
-LV_cs = compute_LV_coordinate_system(LV_grid, ip)
-LV_fm = create_simple_fiber_model(LV_cs, ip_fiber, ip_geo, endo_helix_angle = -60.0, epi_helix_angle = 70.0, endo_transversal_angle = 10.0, epi_transversal_angle = -20.0)
+intorder = max(2*order-1,2)
+ip_collection = LagrangeCollection{order}()
+qr_collection = QuadratureRuleCollection(intorder-1)
+ip = getinterpolation(ip_collection^3, ref_shape)
+ip_fiber = getinterpolation(ip_collection, ref_shape)
+ip_geo = getinterpolation(ip_collection, ref_shape)
+LV_cs = compute_LV_coordinate_system(LV_grid, ip_collection)
+LV_fm = create_simple_microstructure_model(LV_cs, ip_collection^3, ip_collection^3, endo_helix_angle = -60.0, epi_helix_angle = 70.0, endo_transversal_angle = 10.0, epi_transversal_angle = -20.0)
 passive_ho_model = HolzapfelOgden2009Model(1.5806251396691438, 5.8010248271289395, 0.28504197825657906, 4.126552003938297, 0.0, 1.0, 0.0, 1.0, SimpleCompressionPenalty(4.0))
 solve_ideal_lv("lv_test",
     ActiveStressModel(
@@ -212,5 +216,7 @@ solve_ideal_lv("lv_test",
         LV_fm,
     ), LV_grid, LV_cs,
     [NormalSpringBC(0.001, "Epicardium")],
-    ip, ip_geo, max(2*order-1,2)
+    ip_collection,
+    ip, ip_geo,
+    qr_collection
 )

examples/conduction-velocity-benchmark.jl

@@ -39,7 +39,10 @@ end
 
 ######################################################
 
-cs = Thunderbolt.CartesianCoordinateSystemCoefficient(Lagrange{RefHexahedron,1}()^3) # TODO normalize
+order = 1
+ip_collection = LagrangeCollection{order}()
+
+cs = Thunderbolt.CartesianCoordinateSystemCoefficient(ip_collection^3) # TODO normalize
 
 κ₁ = 0.17 * 0.62 / (0.17 + 0.62)
 κᵣ = 0.019 * 0.24 / (0.019 + 0.24)
@@ -62,7 +65,7 @@ mesh = generate_mesh(Hexahedron, (80,28,12), Vec((0.0,0.0,0.0)), Vec((20.0,7.0,3
 
 problem = semidiscretize(
     ReactionDiffusionSplit(model),
-    FiniteElementDiscretization(Dict(:φₘ => LagrangeCollection{1}())),
+    FiniteElementDiscretization(Dict(:φₘ => ip_collection)),
     mesh
 )
 

examples/ring.jl

@@ -154,22 +154,23 @@ function (postproc::StandardMechanicalIOPostProcessor)(t, problem, solver_cache)
     # max_vol = max(max_vol, calculate_volume_deformed_mesh(uₙ,dh,cv));
 end
 
-function solve_test_ring(name_base, constitutive_model, grid, face_models::FM, ip_mech::IPM, ip_geo::IPG, intorder::Int, Δt = 100.0, T = 1000.0) where {ref_shape, IPM <: Interpolation{ref_shape}, IPG <: Interpolation{ref_shape}, FM}
+function solve_test_ring(name_base, constitutive_model, grid, face_models::FM, ip_collection::Thunderbolt.ScalarInterpolationCollection, ip_mech::IPM, ip_geo::IPG, qr_collection::QuadratureRuleCollection, Δt = 100.0, T = 1000.0) where {ref_shape, IPM <: Interpolation{ref_shape}, IPG <: Interpolation{ref_shape}, FM}
     io = ParaViewWriter(name_base);
     # io = JLD2Writer(name_base);
 
     problem = semidiscretize(
         StructuralModel(constitutive_model, face_models),
         FiniteElementDiscretization(
-            Dict(:displacement => LagrangeCollection{1}()^3),
+            Dict(:displacement => ip_collection^3),
             [Dirichlet(:displacement, getfaceset(grid, "Myocardium"), (x,t) -> [0.0], [3])],
             # [Dirichlet(:displacement, getfaceset(grid, "left"), (x,t) -> [0.0], [1]),Dirichlet(:displacement, getfaceset(grid, "front"), (x,t) -> [0.0], [2]),Dirichlet(:displacement, getfaceset(grid, "bottom"), (x,t) -> [0.0], [3]), Dirichlet(:displacement, Set([1]), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3])],
         ),
         grid
     )
 
     # Postprocessor
-    cv_post = CellValues(QuadratureRule{ref_shape}(intorder-1), ip_mech, ip_geo)
+    qr = getquadraturerule(qr_collection, elementtypes(grid)[1])
+    cv_post = CellValues(qr, ip_mech, ip_geo)
     standard_postproc = StandardMechanicalIOPostProcessor(io, cv_post, [1])
 
     # Create sparse matrix and residual vector
@@ -187,19 +188,21 @@ end
 
 ref_shape = RefHexahedron
 order = 1
-
-ip_fsn = Lagrange{ref_shape, order}()^3
-ip_u = Lagrange{ref_shape, order}()^3
-ip_geo = Lagrange{ref_shape, order}()^3
+intorder = 2*order
+ip_collection = LagrangeCollection{order}()
+qr_collection = QuadratureRuleCollection(intorder-1)
+ip_fsn = getinterpolation(ip_collection^3, ref_shape)
+ip_u = getinterpolation(ip_collection^3, ref_shape)
+ip_geo = getinterpolation(ip_collection^3, ref_shape)
 
 ring_grid = generate_ring_mesh(8,2,2)
-ring_cs = compute_midmyocardial_section_coordinate_system(ring_grid, ip_geo)
+ring_cs = compute_midmyocardial_section_coordinate_system(ring_grid, ip_collection^3)
 solve_test_ring("Debug",
     ActiveStressModel(
         Guccione1991PassiveModel(),
         Guccione1993ActiveModel(10.0),
         PelceSunLangeveld1995Model(;calcium_field=CalciumHatField()),
-        create_simple_microstructure_model(ring_cs, ip_fsn, ip_geo,
+        create_simple_microstructure_model(ring_cs, ip_collection^3, ip_collection^3,
             endo_helix_angle = deg2rad(0.0),
             epi_helix_angle = deg2rad(0.0),
             endo_transversal_angle = 0.0,
@@ -208,6 +211,7 @@ solve_test_ring("Debug",
         )
     ), ring_grid, 
     [NormalSpringBC(0.01, "Epicardium")],
-    ip_u, ip_geo, 2*order,
+    ip_collection,
+    ip_u, ip_geo, qr_collection,
     100.0
 )

src/collections.jl

@@ -37,7 +37,7 @@ getinterpolation(lc::LagrangeCollection{order}, ::Type{ref_shape}) where {order,
 
 A collection of fixed-order vectorized Lagrange interpolations across different cell types.
 """
-struct VectorizedInterpolationCollection{vdim, IPC <: ScalarInterpolationCollection} <: InterpolationCollection
+struct VectorizedInterpolationCollection{vdim, IPC <: ScalarInterpolationCollection} <: VectorInterpolationCollection
     base::IPC
     function VectorizedInterpolationCollection{vdim}(ip::SIPC) where {vdim, SIPC <: ScalarInterpolationCollection}
         return new{vdim, SIPC}(ip)
@@ -67,9 +67,13 @@ QuadratureRuleCollection(order::Int) = QuadratureRuleCollection(
     QuadratureRule{RefTetrahedron}(order)
 )
 
-getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Hexahedron) where {QRW, QRH, QRT} = qcr.qr_hex
-getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Wedge) where {QRW, QRH, QRT} = qcr.qr_wedge
-getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Tetrahedron) where {QRW, QRH, QRT} = qcr.qr_tet
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Hexahedron) where {QRW, QRH, QRT} = qrc.qr_hex
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Wedge) where {QRW, QRH, QRT} = qrc.qr_wedge
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Tetrahedron) where {QRW, QRH, QRT} = qrc.qr_tet
+
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Type{Hexahedron}) where {QRW, QRH, QRT} = qrc.qr_hex
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Type{Wedge}) where {QRW, QRH, QRT} = qrc.qr_wedge
+getquadraturerule(qrc::QuadratureRuleCollection{QRW, QRH, QRT}, cell::Type{Tetrahedron}) where {QRW, QRH, QRT} = qrc.qr_tet
 
 """
     CellValueCollection
@@ -82,10 +86,20 @@ struct CellValueCollection{CVW, CVH, CVT}
     cv_tet::CVT
 end
 
+CellValueCollection(qr:: QuadratureRuleCollection, ip::InterpolationCollection, ip_geo::InterpolationCollection = ip) = CellValueCollection(
+    CellValues(getquadraturerule(qr, Wedge), getinterpolation(ip, RefPrism), getinterpolation(ip_geo, RefPrism)), 
+    CellValues(getquadraturerule(qr, Hexahedron), getinterpolation(ip, RefHexahedron), getinterpolation(ip_geo, RefHexahedron)), 
+    CellValues(getquadraturerule(qr, Tetrahedron), getinterpolation(ip, RefTetrahedron), getinterpolation(ip_geo, RefTetrahedron)), 
+)
+
 getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Hexahedron) where {CVW, CVH, CVT} = cv.cv_hex
 getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Wedge) where {CVW, CVH, CVT} = cv.cv_wedge
 getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Tetrahedron) where {CVW, CVH, CVT} = cv.cv_tet
 
+# getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Type{Hexahedron}) where {CVW, CVH, CVT} = cv.cv_hex
+# getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Type{Wedge}) where {CVW, CVH, CVT} = cv.cv_wedge
+# getcellvalues(cv::CellValueCollection{CVW, CVH, CVT}, cell::Type{Tetrahedron}) where {CVW, CVH, CVT} = cv.cv_tet
+
 """
     FaceValueCollection
 

src/mesh/coordinate_systems.jl

@@ -32,10 +32,15 @@ Requires a grid with facesets
 and a nodeset
 * Apex
 """
-function compute_LV_coordinate_system(grid::AbstractGrid, ip_geo::Interpolation{ref_shape}) where {ref_shape <: AbstractRefShape{3}}
-    ip = Lagrange{ref_shape, 1}()
-    qr = QuadratureRule{ref_shape}(2)
-    cellvalues = CellValues(qr, ip, ip_geo);
+function compute_LV_coordinate_system(grid::AbstractGrid{dim}, ip_geo_collection::VectorInterpolationCollection) where {dim}
+    @assert dim == 3
+    @assert length(elementtypes(grid)) == 1
+    ref_shape = getrefshape(getcells(grid, 1))
+    ip_collection = LagrangeCollection{1}()
+    ip = getinterpolation(ip_collection, ref_shape)
+    qr_collection = QuadratureRuleCollection(2)
+    cv_collection = CellValueCollection(qr_collection, ip_collection, ip_geo_collection)
+    cellvalues = getcellvalues(cv_collection, getcells(grid, 1))
 
     dh = DofHandler(grid)
     Ferrite.add!(dh, :coordinates, ip)
@@ -50,7 +55,6 @@ function compute_LV_coordinate_system(grid::AbstractGrid, ip_geo::Interpolation{
     assembler = start_assemble(K)
     @inbounds for cell in CellIterator(dh)
         fill!(Ke, 0)
-
         reinit!(cellvalues, cell)
 
         for qp in QuadratureIterator(cellvalues)
@@ -115,10 +119,15 @@ end
 
 """
 """
-function compute_midmyocardial_section_coordinate_system(grid::AbstractGrid,ip_geo::Interpolation{ref_shape}) where {ref_shape <: AbstractRefShape{3}}
-    ip = Lagrange{ref_shape, 1}()
-    qr = QuadratureRule{ref_shape}(2)
-    cellvalues = CellValues(qr, ip, ip_geo);
+function compute_midmyocardial_section_coordinate_system(grid::AbstractGrid{dim}, ip_geo_collection::VectorInterpolationCollection) where {dim}
+    @assert dim == 3
+    @assert length(elementtypes(grid)) == 1
+    ref_shape = getrefshape(getcells(grid,1))
+    ip_collection = LagrangeCollection{1}()
+    ip = getinterpolation(ip_collection, ref_shape)
+    qr_collection = QuadratureRuleCollection(2)
+    cv_collection = CellValueCollection(qr_collection, ip_collection, ip_geo_collection)
+    cellvalues = getcellvalues(cv_collection, getcells(grid,1))
 
     dh = DofHandler(grid)
     Ferrite.add!(dh, :coordinates, ip)

src/modeling/core/coefficients.jl

@@ -3,16 +3,16 @@
 
 A constant in time data field, interpolated per element with a given interpolation.
 """
-struct FieldCoefficient{TA,IP<:Interpolation}
+struct FieldCoefficient{TA,IPC<:InterpolationCollection}
     # TODO data structure for this
     elementwise_data::TA #3d array (element_idx, base_fun_idx, dim)
-    ip::IP
+    ip_collection::IPC
     # TODO use CellValues
 end
 
-function evaluate_coefficient(coeff::FieldCoefficient{<:Any,<:ScalarInterpolation}, cell_cache, qp::QuadraturePoint{<:Any, T}, t) where T
-    @unpack elementwise_data, ip = coeff
-
+function evaluate_coefficient(coeff::FieldCoefficient{<:Any,<:ScalarInterpolationCollection}, cell_cache, qp::QuadraturePoint{<:Any, T}, t) where T
+    @unpack elementwise_data, ip_collection = coeff
+    ip = getinterpolation(ip_collection, getcells(cell_cache.grid, cellid(cell_cache)))
     n_base_funcs = Ferrite.getnbasefunctions(ip)
     val = zero(T)
 
@@ -22,9 +22,9 @@ function evaluate_coefficient(coeff::FieldCoefficient{<:Any,<:ScalarInterpolatio
     return val
 end
 
-function evaluate_coefficient(coeff::FieldCoefficient{<:Any,<:VectorizedInterpolation{vdim}}, cell_cache, qp::QuadraturePoint{<:Any, T}, t) where {vdim,T}
-    @unpack elementwise_data, ip = coeff
-
+function evaluate_coefficient(coeff::FieldCoefficient{<:Any,<:VectorizedInterpolationCollection{vdim}}, cell_cache, qp::QuadraturePoint{<:Any, T}, t) where {vdim,T}
+    @unpack elementwise_data, ip_collection = coeff
+    ip = getinterpolation(ip_collection, getcells(cell_cache.grid, cellid(cell_cache)))
     n_base_funcs = Ferrite.getnbasefunctions(ip.ip)
     val = zero(Vec{vdim, T})
 
@@ -61,14 +61,16 @@ function evaluate_coefficient(coeff::ConductivityToDiffusivityCoefficient, cell_
     return κ/(Cₘ*χ)
 end
 
-struct CartesianCoordinateSystemCoefficient{IP<:VectorizedInterpolation}
-    ip::IP
+struct CartesianCoordinateSystemCoefficient{IPC<:VectorizedInterpolationCollection}
+    ip_collection::IPC
 end
 
-function evaluate_coefficient(coeff::CartesianCoordinateSystemCoefficient{<:VectorizedInterpolation{sdim}}, cell_cache, qp::QuadraturePoint{<:Any,T}, t) where {sdim, T}
+function evaluate_coefficient(coeff::CartesianCoordinateSystemCoefficient{<:VectorizedInterpolationCollection{sdim}}, cell_cache, qp::QuadraturePoint{<:Any,T}, t) where {sdim, T}
     x = zero(Vec{sdim, T})
-    for i in 1:getnbasefunctions(coeff.ip.ip)
-        x += Ferrite.shape_value(coeff.ip.ip, qp.ξ, i) * cell_cache.coords[i]
+    ip_collection = coeff.ip_collection
+    ip = getinterpolation(ip_collection, getcells(cell_cache.grid, cellid(cell_cache)))
+    for i in 1:getnbasefunctions(ip.ip)
+        x += Ferrite.shape_value(ip.ip, qp.ξ, i) * cell_cache.coords[i]
     end
     return x
 end