Merge branch 'main' into do/pseudo-ecg

examples/LV.jl

@@ -59,12 +59,11 @@ function (postproc::StandardMechanicalIOPostProcessor2)(t, problem, solver_cache
             helixangleref_cell = 0.0
 
             nqp = getnquadpoints(cv)
-            for qpᵢ in 1:nqp
-                qp = QuadraturePoint(qpᵢ, cv.qr.points[qpᵢ])
-                dΩ = getdetJdV(cv, qpᵢ)
+            for qp in QuadratureIterator(cv)
+                dΩ = getdetJdV(cv, qp)
 
                 # Compute deformation gradient F
-                ∇u = function_gradient(cv, qpᵢ, uₑ)
+                ∇u = function_gradient(cv, qp, uₑ)
                 F = one(∇u) + ∇u
 
                 C = tdot(F)
@@ -80,7 +79,7 @@ function (postproc::StandardMechanicalIOPostProcessor2)(t, problem, solver_cache
                 s₀_current /= norm(s₀_current)
 
                 coords = getcoordinates(cell)
-                x_global = spatial_coordinate(cv, qpᵢ, coords)
+                x_global = spatial_coordinate(cv, qp, coords)
 
                 # v_longitudinal = function_gradient(cv_cs, qp, coordinate_system.u_apicobasal[celldofs(cell)])
                 # v_radial = function_gradient(cv_cs, qp, coordinate_system.u_transmural[celldofs(cell)])

examples/ring.jl

@@ -65,12 +65,11 @@ function (postproc::StandardMechanicalIOPostProcessor)(t, problem, solver_cache)
             helixangleref_cell = 0.0
 
             nqp = getnquadpoints(cv)
-            for qpᵢ in 1:nqp
-                qp = QuadraturePoint(qpᵢ, cv.qr.points[qpᵢ])
-                dΩ = getdetJdV(cv, qpᵢ)
+            for qp in QuadratureIterator(cv)
+                dΩ = getdetJdV(cv, qp)
 
                 # Compute deformation gradient F
-                ∇u = function_gradient(cv, qpᵢ, uₑ)
+                ∇u = function_gradient(cv, qp, uₑ)
                 F = one(∇u) + ∇u
 
                 C = tdot(F)
@@ -86,7 +85,7 @@ function (postproc::StandardMechanicalIOPostProcessor)(t, problem, solver_cache)
                 s₀_current /= norm(s₀_current)
 
                 coords = getcoordinates(cell)
-                x_global = spatial_coordinate(cv, qpᵢ, coords)
+                x_global = spatial_coordinate(cv, qp, coords)
 
                 # v_longitudinal = function_gradient(cv_cs, qp, coordinate_system.u_apicobasal[celldofs(cell)])
                 # v_radial = function_gradient(cv_cs, qp, coordinate_system.u_transmural[celldofs(cell)])

src/io.jl

@@ -124,7 +124,7 @@ function store_green_lagrange!(io::ParaViewWriter, dh, u::AbstractVector, a_coef
         field_dofs  = dof_range(sdh, field_idx)
         uₑ = @view u[global_dofs] # element dofs
         for qp in QuadratureIterator(cv)
-            ∇u = function_gradient(cv, qpᵢ, uₑ)
+            ∇u = function_gradient(cv, qp, uₑ)
 
             F = one(∇u) + ∇u
             C = tdot(F)

src/modeling/coupler/fsi.jl

@@ -34,7 +34,7 @@ end
 
 Compute the chamber volume as a surface integral via the integral
  - ∫ (x + d) det(F) adj(F) N ∂Ωendo
-where it is assumet that the chamber is convex, zero displacement in 
+where it is assumed that the chamber is convex, zero displacement in 
 apicobasal direction at the valvular plane occurs and the plane normal is aligned
 with the z axis, where the origin is at z=0.
 """
@@ -78,7 +78,7 @@ function assemble_interface_coupling_contribution!(C, r, dh, u, setname, method:
             ∂V∂F = Tensors.gradient(u -> volume_integral(x, d, u, N, method), F)
             for j ∈ 1:getnbasefunctions(fv)
                 δuⱼ = shape_value(fv, qp, j)
-                ∇δuj = shape_gradient(cv, qpᵢ, j)
+                ∇δuj = shape_gradient(cv, qp, j)
                 C[1, ddofs[j]] += (∂V∂u ⋅ δuⱼ + ∂V∂F ⊡ ∇δuj) * dΓ
             end
         end

src/modeling/microstructure.jl

@@ -25,6 +25,7 @@ function evaluate_coefficient(fsn::OrthotropicMicrostructureModel, cell_cache, q
     f, s, n
 end
 
+# TODO: citation
 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)
@@ -94,9 +95,9 @@ function create_simple_microstructure_model(coordinate_system, ip::VectorInterpo
             transversal_angle = (1-transmural) * endo_transversal_angle + (transmural) * epi_transversal_angle
 
             f₀, s₀, n₀ = streeter_type_fsn(transmural_direction, circumferential_direction, apicobasal_direction, helix_angle, transversal_angle, sheetlet_pseudo_angle, make_orthogonal)
-            elementwise_data_f[cellindex, qp] = f₀
-            elementwise_data_s[cellindex, qp] = s₀
-            elementwise_data_n[cellindex, qp] = n₀
+            elementwise_data_f[cellindex, qp.i] = f₀
+            elementwise_data_s[cellindex, qp.i] = s₀
+            elementwise_data_n[cellindex, qp.i] = n₀
         end
     end
 

src/quadrature_iterator.jl

@@ -30,6 +30,8 @@ end
 Base.eltype(::Type{<:QuadraturePoint{<:QuadratureRule{<:AbstractRefShape, dim, T}}}) where {dim, T} = QuadraturePoint{dim, T}
 Base.length(iterator::QuadratureIterator) = length(Ferrite.getnquadpoints(iterator.qr))
 
+Ferrite.spatial_coordinate(fe_v::Ferrite.AbstractValues, qp::QuadraturePoint, x::AbstractVector{<:Vec}) = Ferrite.spatial_coordinate(fe_v, qp.i, x)
+
 Ferrite.getdetJdV(cv::CellValues, qp::QuadraturePoint) = Ferrite.getdetJdV(cv, qp.i)
 Ferrite.shape_value(cv::CellValues, qp::QuadraturePoint, base_fun_idx::Int) = Ferrite.shape_value(cv, qp.i, base_fun_idx)
 Ferrite.shape_gradient(cv::CellValues, qp::QuadraturePoint, base_fun_idx::Int) = Ferrite.shape_gradient(cv, qp.i, base_fun_idx)

test/runtests.jl

@@ -5,6 +5,7 @@ include("test_operators.jl")
 include("test_type_stability.jl")
 include("test_mesh.jl")
 include("test_coefficients.jl")
+include("test_microstructures.jl")
 
 include("integration/test_contracting_cuboid.jl")
 include("integration/test_waveprop_cuboid.jl")

test/test_microstructures.jl

@@ -0,0 +1,53 @@
+@testset "microstructures" begin
+    ref_shape = RefHexahedron
+    order = 1
+
+    qr = QuadratureRule{ref_shape}(2)
+
+    ip_fsn = Lagrange{ref_shape, order}()^3
+    ip_geo = Lagrange{ref_shape, order}()^3
+    cv_fsn = CellValues(qr, ip_fsn, ip_geo)
+
+    ring_grid = generate_ring_mesh(80,1,1)
+    ring_cs = compute_midmyocardial_section_coordinate_system(ring_grid, ip_geo)
+
+    cv_cs = Thunderbolt.create_cellvalues(ring_cs, qr)
+
+    @testset "Midmyocardial coordinate system" begin
+        for cellcache in CellIterator(ring_cs.dh)
+            reinit!(cv_cs, cellcache)
+            dof_indices = celldofs(cellcache)
+            for qp in QuadratureIterator(cv_fsn)
+                coords = getcoordinates(cellcache)
+                sc = spatial_coordinate(cv_fsn, qp, coords)
+                transmural = 4*(norm(Vec(sc[1:2]..., 0.))-0.75)
+                @test transmural ≈ function_value(cv_cs, qp, ring_cs.u_transmural[dof_indices]) atol=0.01
+            end
+        end
+    end
+
+    @testset "OrthotropicMicrostructureModel" begin
+        ms = create_simple_microstructure_model(ring_cs, ip_fsn, ip_geo,
+            endo_helix_angle = deg2rad(0.0),
+            epi_helix_angle = deg2rad(0.0),
+            endo_transversal_angle = 0.0,
+            epi_transversal_angle = 0.0,
+            sheetlet_pseudo_angle = deg2rad(0)
+        )
+        for cellcache in CellIterator(ring_grid)
+            reinit!(cv_fsn, cellcache)
+            for qp in QuadratureIterator(cv_fsn)
+                coords = getcoordinates(cellcache)
+                sc = spatial_coordinate(cv_fsn, qp, coords)
+                # If we set all angles to 0, then the generator on the ring simply generates sheetlets which point in positive z direction, where the normal point radially outwards.
+                ndir = Vec(sc[1:2]..., 0.)/norm(sc[1:2])
+                sdir = Vec(0., 0., 1.)
+                fdir = sdir × ndir
+                fsn = evaluate_coefficient(ms, cellcache, qp, 0)
+                @test fsn[1] ≈ fdir atol=0.05
+                @test fsn[2] ≈ sdir
+                @test fsn[3] ≈ ndir atol=0.05
+            end
+        end
+    end
+end