Quadrature Point Iterator (#53)

termi-official · Dec 12, 2023 · 38c0de7 · 38c0de7
1 parent e04bd2f
commit 38c0de7
Show file tree

Hide file tree

Showing 14 changed files with 98 additions and 74 deletions.
diff --git a/docs/src/api-reference/utility.md b/docs/src/api-reference/utility.md
@@ -21,13 +21,13 @@ FaceValueCollection
 
 ## Iteration
 
-TODO QuadraturePointIterator https://github.com/termi-official/Thunderbolt.jl/issues/27
-TODO TimeChoiceIterator https://github.com/termi-official/Thunderbolt.jl/issues/32
-
-```@devdocs
+```@docs
 QuadraturePoint
+QuadratureIterator
 ```
 
+TODO TimeChoiceIterator https://github.com/termi-official/Thunderbolt.jl/issues/32
+
 ## IO
 
 ```@docs

diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl
@@ -145,6 +145,7 @@ export
     calculate_volume_deformed_mesh,
     elementtypes,
     QuadraturePoint,
+    QuadratureIterator,
     # IO
     ParaViewWriter,
     JLD2Writer,

diff --git a/src/io.jl b/src/io.jl
@@ -66,7 +66,7 @@ function store_coefficient!(io::ParaViewWriter, t, coefficient::ConstantCoeffici
     qrc = QuadratureRuleCollection(1)
     for cell_cache in CellIterator(dh)
         qr = getquadraturerule(qr_collection, getcells(get_grid(dh), cellid(cell_cache)))
-        data[cellid(cell_cache)] = evaluate_coefficient(coefficient, cell_cache, QuadraturePoint(1, getpoints(qr)[1]), t)
+        data[cellid(cell_cache)] = evaluate_coefficient(coefficient, cell_cache, first(QuadratureIterator(qr)), t)
     end
     vtk_cell_data(io.current_file, data, name)
 end
@@ -82,8 +82,7 @@ function _store_coefficient!(::Union{Type{<:Tuple{T}},Type{<:SVector{T}}}, tlen:
     data = zeros(T, getncells(grid), tlen)
     for cell_cache in CellIterator(dh) # TODO subdomain support
         qr = getquadraturerule(qr_collection, getcells(get_grid(dh), cellid(cell_cache)))
-        for qpᵢ ∈ 1:getnquadpoints(qr) # TODO iterator
-            qp = QuadraturePoint(qpᵢ, getpoints(qr)[qpᵢ])
+        for qp ∈ QuadratureIterator(qr)
             tval = evaluate_coefficient(coefficient, cell_cache, qp, t)
             for i ∈ 1:tlen
                 data[cellid(cell_cache), i] += tval[I]
@@ -100,8 +99,7 @@ function _store_coefficient!(T::Type, tlen::Int, io::ParaViewWriter, dh, coeffic
     data = zeros(T, getncells(grid))
     for cell_cache in CellIterator(dh) # TODO subdomain support
         qr = getquadraturerule(qr_collection, getcells(get_grid(dh), cellid(cell_cache)))
-        for qpᵢ ∈ 1:getnquadpoints(qr) # TODO iterator
-            qp = QuadraturePoint(qpᵢ, getpoints(qr)[qpᵢ])
+        for qp ∈ QuadratureIterator(qr)
             data[cellid(cell_cache)] += evaluate_coefficient(coefficient, cell_cache, qp, t)
         end
         data[cellid(cell_cache)] /= getnquadpoints(qr)
@@ -125,9 +123,7 @@ function store_green_lagrange!(io::ParaViewWriter, dh, u::AbstractVector, a_coef
         global_dofs = celldofs(cell_cache)
         field_dofs  = dof_range(sdh, field_idx)
         uₑ = @view u[global_dofs] # element dofs
-        nqp = getnquadpoints(cv)
-        for qpᵢ in 1:nqp
-            qp = QuadraturePoint(qpᵢ, cv.qr.points[qpᵢ])
+        for qp in QuadratureIterator(cv)
             ∇u = function_gradient(cv, qpᵢ, uₑ)
 
             F = one(∇u) + ∇u
@@ -139,7 +135,7 @@ function store_green_lagrange!(io::ParaViewWriter, dh, u::AbstractVector, a_coef
 
             E[cellid(cell)] += a ⋅ E ⋅ b
         end
-        E[cellid(cell)] /= nqp
+        E[cellid(cell)] /= getnquadpoints(cv)
     end
     vtk_cell_data(io.current_file, E, name)
 end

diff --git a/src/mesh/coordinate_systems.jl b/src/mesh/coordinate_systems.jl
@@ -49,13 +49,13 @@ function compute_LV_coordinate_system(grid::AbstractGrid, ip_geo::Interpolation{
 
         reinit!(cellvalues, cell)
 
-        for q_point in 1:getnquadpoints(cellvalues)
-            dΩ = getdetJdV(cellvalues, q_point)
+        for qp in QuadratureIterator(cellvalues)
+            dΩ = getdetJdV(cellvalues, qp)
 
             for i in 1:n_basefuncs
-                ∇v = shape_gradient(cellvalues, q_point, i)
+                ∇v = shape_gradient(cellvalues, qp, i)
                 for j in 1:n_basefuncs
-                    ∇u = shape_gradient(cellvalues, q_point, j)
+                    ∇u = shape_gradient(cellvalues, qp, j)
                     Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                 end
             end
@@ -132,13 +132,13 @@ function compute_midmyocardial_section_coordinate_system(grid::AbstractGrid,ip_g
 
         reinit!(cellvalues, cell)
 
-        for q_point in 1:getnquadpoints(cellvalues)
-            dΩ = getdetJdV(cellvalues, q_point)
+        for qp in QuadratureIterator(cellvalues)
+            dΩ = getdetJdV(cellvalues, qp)
 
             for i in 1:n_basefuncs
-                ∇v = shape_gradient(cellvalues, q_point, i)
+                ∇v = shape_gradient(cellvalues, qp, i)
                 for j in 1:n_basefuncs
-                    ∇u = shape_gradient(cellvalues, q_point, j)
+                    ∇u = shape_gradient(cellvalues, qp, j)
                     Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                 end
             end

diff --git a/src/modeling/core/boundary_conditions.jl b/src/modeling/core/boundary_conditions.jl
@@ -48,7 +48,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -76,7 +76,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Robin
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -103,7 +103,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -131,7 +131,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Norma
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -158,7 +158,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -190,7 +190,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Bendi
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
         N = getnormal(fv, qp)
 
@@ -220,9 +220,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
-        ξ = Ferrite.getpoints(fv.qr, face[2])[qp]
-
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
 
         n₀ = getnormal(fv, qp)
@@ -237,7 +235,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
         # Add contribution to the residual from this test function
         cofF = transpose(inv(F))
         # TODO fix the "nothing" here
-        neumann_term = evaluate_coefficient(p, nothing, QuadraturePoint(qp, ξ), time) * det(F) * cofF
+        neumann_term = evaluate_coefficient(p, nothing, qp, time) * det(F) * cofF
         for i in 1:ndofs_face
             δuᵢ = shape_value(fv, qp, i)
             residualₑ[i] += neumann_term ⋅ n₀ ⋅ δuᵢ * dΓ
@@ -261,9 +259,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Press
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
-        ξ = Ferrite.getpoints(fv.qr, face[2])[qp]
-
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
 
         n₀ = getnormal(fv, qp)
@@ -278,7 +274,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Press
         # Add contribution to the residual from this test function
         cofF = transpose(inv(F))
         # TODO fix the "nothing" here
-        neumann_term = evaluate_coefficient(p, nothing, QuadraturePoint(qp, ξ), time) * det(F) * cofF
+        neumann_term = evaluate_coefficient(p, nothing, qp, time) * det(F) * cofF
         for i in 1:ndofs_face
             δuᵢ = shape_value(fv, qp, i)
 
@@ -301,7 +297,7 @@ function assemble_face!(Kₑ::Matrix, residualₑ, uₑ, face, cache::SimpleFace
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
 
         n₀ = getnormal(fv, qp)
@@ -339,7 +335,7 @@ function assemble_face!(Kₑ::Matrix, uₑ, face, cache::SimpleFaceCache{<:Const
     reinit!(fv, face[1], face[2])
 
     ndofs_face = getnbasefunctions(fv)
-    for qp in 1:getnquadpoints(fv)
+    for qp in QuadratureIterator(fv)
         dΓ = getdetJdV(fv, qp)
 
         n₀ = getnormal(fv, qp)

diff --git a/src/modeling/core/coefficients.jl b/src/modeling/core/coefficients.jl
@@ -68,7 +68,7 @@ end
 function evaluate_coefficient(coeff::CartesianCoordinateSystemCoefficient{<:VectorizedInterpolation{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 += shape_value(coeff.ip.ip, qp.ξ, i) * cell_cache.coords[i]
+        x += Ferrite.shape_value(coeff.ip.ip, qp.ξ, i) * cell_cache.coords[i]
     end
     return x
 end

diff --git a/src/modeling/core/diffusion.jl b/src/modeling/core/diffusion.jl
@@ -22,15 +22,13 @@ function assemble_element!(Kₑ, cell, element_cache::CACHE, time) where {CACHE
 
     reinit!(cellvalues, cell)
 
-    for q_point in 1:getnquadpoints(cellvalues)
-        ξ = cellvalues.qr.points[q_point]
-        qp = QuadraturePoint(q_point, ξ)
+    for qp in QuadratureIterator(cellvalues)
         D_loc = evaluate_coefficient(element_cache.integrator.D, cell, qp, time)
-        dΩ = getdetJdV(cellvalues, q_point)
+        dΩ = getdetJdV(cellvalues, qp)
         for i in 1:n_basefuncs
-            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)
+            ∇Nᵢ = shape_gradient(cellvalues, qp, i)
             for j in 1:n_basefuncs
-                ∇Nⱼ = shape_gradient(cellvalues, q_point, j)
+                ∇Nⱼ = shape_gradient(cellvalues, qp, j)
                 Kₑ[i,j] -= ((D_loc ⋅ ∇Nᵢ) ⋅ ∇Nⱼ) * dΩ
             end
         end

diff --git a/src/modeling/core/mass.jl b/src/modeling/core/mass.jl
@@ -21,15 +21,13 @@ function assemble_element!(Mₑ, cell, element_cache::CACHE, time) where {CACHE
     @unpack cellvalues = element_cache
     reinit!(element_cache.cellvalues, cell)
     n_basefuncs = getnbasefunctions(cellvalues)
-    for q_point in 1:getnquadpoints(cellvalues)
-        ξ = cellvalues.qr.points[q_point]
-        qp = QuadraturePoint(q_point, ξ)
+    for qp in QuadratureIterator(cellvalues)
         ρ = evaluate_coefficient(element_cache.integrator.ρ, cell, qp, time)
-        dΩ = getdetJdV(cellvalues, q_point)
+        dΩ = getdetJdV(cellvalues, qp)
         for i in 1:n_basefuncs
-            Nᵢ = shape_value(cellvalues, q_point, i)
+            Nᵢ = shape_value(cellvalues, qp, i)
             for j in 1:n_basefuncs
-                Nⱼ = shape_value(cellvalues, q_point, j)
+                Nⱼ = shape_value(cellvalues, qp, j)
                 Mₑ[i,j] += ρ * Nᵢ * Nⱼ * dΩ 
             end
         end

diff --git a/src/modeling/coupler/fsi.jl b/src/modeling/coupler/fsi.jl
@@ -60,7 +60,7 @@ function assemble_interface_coupling_contribution!(C, r, dh, u, setname, method:
         ddofs = @view celldofs(face)[drange]
         uₑ = @view u[ddofs]
 
-        for qp in 1:getnquadpoints(fv)
+        for qp in QuadratureIterator(fv)
             dΓ = getdetJdV(fv, qp)
             N = getnormal(fv, qp)
 
@@ -103,7 +103,7 @@ function compute_chamber_volume(dh, u, setname, method)
         ddofs = @view celldofs(face)[drange]
         uₑ = @view u[ddofs]
 
-        for qp in 1:getnquadpoints(fv)
+        for qp in QuadratureIterator(fv)
             dΓ = getdetJdV(fv, qp)
             N = getnormal(fv, qp)
 

diff --git a/src/modeling/microstructure.jl b/src/modeling/microstructure.jl
@@ -78,7 +78,7 @@ function create_simple_microstructure_model(coordinate_system, ip::VectorInterpo
         reinit!(cv, cell)
         dof_indices = celldofs(cell)
 
-        for qp in 1:getnquadpoints(cv)
+        for qp in QuadratureIterator(cv)
             # TODO grab these via some interface!
             apicobasal_direction = function_gradient(cv, qp, coordinate_system.u_apicobasal[dof_indices])
             apicobasal_direction /= norm(apicobasal_direction)

diff --git a/src/modeling/solid/elements.jl b/src/modeling/solid/elements.jl
@@ -30,13 +30,11 @@ function assemble_element!(Kₑ::Matrix, residualₑ, uₑ, geometry_cache, elem
 
     reinit!(cv, geometry_cache)
 
-    @inbounds for qpᵢ in 1:getnquadpoints(cv)
-        ξ = cv.qr.points[qpᵢ]
-        qp = QuadraturePoint(qpᵢ, ξ)
-        dΩ = getdetJdV(cv, qpᵢ)
+    @inbounds for qp ∈ 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
 
         # Compute stress and tangent
@@ -45,14 +43,14 @@ function assemble_element!(Kₑ::Matrix, residualₑ, uₑ, geometry_cache, elem
 
         # Loop over test functions
         for i in 1:ndofs
-            ∇δui = shape_gradient(cv, qpᵢ, i)
+            ∇δui = shape_gradient(cv, qp, i)
 
             # Add contribution to the residual from this test function
             residualₑ[i] += ∇δui ⊡ P * dΩ
 
             ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
             for j in 1:ndofs
-                ∇δuj = shape_gradient(cv, qpᵢ, j)
+                ∇δuj = shape_gradient(cv, qp, j)
                 # Add contribution to the tangent
                 Kₑ[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
             end

diff --git a/src/quadrature_iterator.jl b/src/quadrature_iterator.jl
@@ -0,0 +1,44 @@
+"""
+    QuadraturePoint{dim, T}
+
+A simple helper to carry quadrature point information.
+"""
+struct QuadraturePoint{dim, T}
+    i::Int
+    ξ::Vec{dim, T}
+end
+
+"""
+    QuadratureIterator(::QuadratureRule)
+    QuadratureIterator(::FaceQuadratureRule, local_face_idx::Int)
+    QuadratureIterator(::CellValues)
+    QuadratureIterator(::FaceValues)
+
+A helper to loop over the quadrature points in some rule or cache with type [`QuadraturePoint`](@ref).
+"""
+struct QuadratureIterator{QR<:QuadratureRule}
+    qr::QR
+end
+QuadratureIterator(fqr::FaceQuadratureRule, local_face_idx::Int) = QuadratureIterator(fqr.face_rules[local_face_idx])
+QuadratureIterator(cv::CellValues) = QuadratureIterator(cv.qr)
+QuadratureIterator(fv::FaceValues) = QuadratureIterator(fv.fqr.face_rules[fv.current_face[]])
+
+function Base.iterate(iterator::QuadratureIterator, i = 1)
+    i > getnquadpoints(iterator.qr) && return nothing
+    return (QuadraturePoint(i,Ferrite.getpoints(iterator.qr)[i]), i + 1)
+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.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)
+Ferrite.function_value(cv::CellValues, qp::QuadraturePoint, ue) = Ferrite.function_value(cv, qp.i, ue)
+Ferrite.function_gradient(cv::CellValues, qp::QuadraturePoint, ue) = Ferrite.function_gradient(cv, qp.i, ue)
+
+Ferrite.getdetJdV(fv::FaceValues, qp::QuadraturePoint) = Ferrite.getdetJdV(fv, qp.i)
+Ferrite.shape_value(fv::FaceValues, qp::QuadraturePoint, base_fun_idx::Int) = Ferrite.shape_value(fv, qp.i, base_fun_idx)
+Ferrite.shape_gradient(fv::FaceValues, qp::QuadraturePoint, base_fun_idx::Int) = Ferrite.shape_gradient(fv, qp.i, base_fun_idx)
+Ferrite.function_value(fv::FaceValues, qp::QuadraturePoint, ue) = Ferrite.function_value(fv, qp.i, ue)
+Ferrite.function_gradient(fv::FaceValues, qp::QuadraturePoint, ue) = Ferrite.function_gradient(fv, qp.i, ue)
+Ferrite.getnormal(fv::FaceValues, qp::QuadraturePoint) = Ferrite.getnormal(fv, qp.i)