From cfe3863c5b2553a042ad25a6f379c847e1354b24 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 11 Dec 2023 01:23:55 +0100 Subject: [PATCH] More flexibility for internal numeric types (#802) --- Project.toml | 3 +- src/Dofs/ConstraintHandler.jl | 17 +++-- src/FEValues/FaceValues.jl | 2 +- src/FEValues/PointValues.jl | 2 +- src/assembler.jl | 6 +- src/interpolations.jl | 129 +++++++++++++++++----------------- src/iterators.jl | 4 +- test/runtests.jl | 11 +++ test/test_cellvalues.jl | 7 ++ test/test_dofs.jl | 1 + test/test_interpolations.jl | 43 ++++++++---- 11 files changed, 134 insertions(+), 91 deletions(-) diff --git a/Project.toml b/Project.toml index 0d44eb9d16..c18ee1634a 100644 --- a/Project.toml +++ b/Project.toml @@ -39,6 +39,7 @@ FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" @@ -49,4 +50,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [targets] -test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"] +test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "JET", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"] diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 8d9b9f1cec..9d013aa2fc 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -69,9 +69,10 @@ struct AffineConstraint{T} end """ - ConstraintHandler + ConstraintHandler([T=Float64], dh::AbstractDofHandler) -Collection of constraints. +A collection of constraints associated with the dof handler `dh`. +`T` is the numeric type for stored values. """ struct ConstraintHandler{DH<:AbstractDofHandler,T} dbcs::Vector{Dirichlet} @@ -90,11 +91,13 @@ struct ConstraintHandler{DH<:AbstractDofHandler,T} closed::ScalarWrapper{Bool} end -function ConstraintHandler(dh::AbstractDofHandler) +ConstraintHandler(dh::AbstractDofHandler) = ConstraintHandler(Float64, dh) + +function ConstraintHandler(::Type{T}, dh::AbstractDofHandler) where T <: Number @assert isclosed(dh) ConstraintHandler( - Dirichlet[], Int[], Int[], Float64[], Union{Nothing,Float64}[], Union{Nothing,DofCoefficients{Float64}}[], - Dict{Int,Int}(), BCValues{Float64}[], dh, ScalarWrapper(false), + Dirichlet[], Int[], Int[], T[], Union{Nothing,T}[], Union{Nothing,DofCoefficients{T}}[], + Dict{Int,Int}(), BCValues{T}[], dh, ScalarWrapper(false), ) end @@ -407,7 +410,7 @@ function update!(ch::ConstraintHandler, time::Real=0.0) end # for vertices, faces and edges -function _update!(inhomogeneities::Vector{Float64}, f::Function, boundary_entities::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int}, +function _update!(inhomogeneities::Vector{T}, f::Function, boundary_entities::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int}, components::Vector{Int}, dh::AbstractDofHandler, boundaryvalues::BCValues, dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where {T} @@ -444,7 +447,7 @@ function _update!(inhomogeneities::Vector{Float64}, f::Function, boundary_entiti end # for nodes -function _update!(inhomogeneities::Vector{Float64}, f::Function, ::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int}, +function _update!(inhomogeneities::Vector{T}, f::Function, ::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int}, components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues, dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T counter = 1 diff --git a/src/FEValues/FaceValues.jl b/src/FEValues/FaceValues.jl index 69bf96e29f..64d944e2df 100644 --- a/src/FEValues/FaceValues.jl +++ b/src/FEValues/FaceValues.jl @@ -6,7 +6,7 @@ values of nodal functions, gradients and divergences of nodal functions etc. on **Arguments:** -* `T`: an optional argument to determine the type the internal data is stored as. +* `T`: an optional argument (default to `Float64`) to determine the type the internal data is stored as. * `quad_rule`: an instance of a [`FaceQuadratureRule`](@ref) * `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function * `geom_interpol`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry. diff --git a/src/FEValues/PointValues.jl b/src/FEValues/PointValues.jl index 62107b4343..ab0f98e353 100644 --- a/src/FEValues/PointValues.jl +++ b/src/FEValues/PointValues.jl @@ -1,6 +1,6 @@ """ PointValues(cv::CellValues) - PointValues(ip_f::Interpolation, ip_g::Interpolation=ip_f) + PointValues([::Type{T}], func_interpol::Interpolation, [geom_interpol::Interpolation]) Similar to `CellValues` but with a single updateable "quadrature point". `PointValues` are used for evaluation of functions/gradients in diff --git a/src/assembler.jl b/src/assembler.jl index 69307a0782..573f4437dc 100644 --- a/src/assembler.jl +++ b/src/assembler.jl @@ -4,10 +4,12 @@ struct Assembler{T} V::Vector{T} end -function Assembler(N) +Assembler(N) = Assembler(Float64, N) + +function Assembler(::Type{T}, N) where T I = Int[] J = Int[] - V = Float64[] + V = T[] sizehint!(I, N) sizehint!(J, N) sizehint!(V, N) diff --git a/src/interpolations.jl b/src/interpolations.jl index 0d74256e44..537759b15c 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -545,8 +545,8 @@ end function shape_value(ip::Lagrange{RefLine, 1}, ξ::Vec{1}, i::Int) ξ_x = ξ[1] - i == 1 && return (1 - ξ_x) * 0.5 - i == 2 && return (1 + ξ_x) * 0.5 + i == 1 && return (1 - ξ_x) / 2 + i == 2 && return (1 + ξ_x) / 2 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -566,8 +566,8 @@ end function shape_value(ip::Lagrange{RefLine, 2}, ξ::Vec{1}, i::Int) ξ_x = ξ[1] - i == 1 && return ξ_x * (ξ_x - 1) * 0.5 - i == 2 && return ξ_x * (ξ_x + 1) * 0.5 + i == 1 && return ξ_x * (ξ_x - 1) / 2 + i == 2 && return ξ_x * (ξ_x + 1) / 2 i == 3 && return 1 - ξ_x^2 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -589,10 +589,10 @@ end function shape_value(ip::Lagrange{RefQuadrilateral, 1}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] - i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25 - i == 2 && return (1 + ξ_x) * (1 - ξ_y) * 0.25 - i == 3 && return (1 + ξ_x) * (1 + ξ_y) * 0.25 - i == 4 && return (1 - ξ_x) * (1 + ξ_y) * 0.25 + i == 1 && return (1 - ξ_x) * (1 - ξ_y) / 4 + i == 2 && return (1 + ξ_x) * (1 - ξ_y) / 4 + i == 3 && return (1 + ξ_x) * (1 + ξ_y) / 4 + i == 4 && return (1 - ξ_x) * (1 + ξ_y) / 4 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -620,14 +620,14 @@ end function shape_value(ip::Lagrange{RefQuadrilateral, 2}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] - i == 1 && return (ξ_x^2 - ξ_x) * (ξ_y^2 - ξ_y) * 0.25 - i == 2 && return (ξ_x^2 + ξ_x) * (ξ_y^2 - ξ_y) * 0.25 - i == 3 && return (ξ_x^2 + ξ_x) * (ξ_y^2 + ξ_y) * 0.25 - i == 4 && return (ξ_x^2 - ξ_x) * (ξ_y^2 + ξ_y) * 0.25 - i == 5 && return (1 - ξ_x^2) * (ξ_y^2 - ξ_y) * 0.5 - i == 6 && return (ξ_x^2 + ξ_x) * (1 - ξ_y^2) * 0.5 - i == 7 && return (1 - ξ_x^2) * (ξ_y^2 + ξ_y) * 0.5 - i == 8 && return (ξ_x^2 - ξ_x) * (1 - ξ_y^2) * 0.5 + i == 1 && return (ξ_x^2 - ξ_x) * (ξ_y^2 - ξ_y) / 4 + i == 2 && return (ξ_x^2 + ξ_x) * (ξ_y^2 - ξ_y) / 4 + i == 3 && return (ξ_x^2 + ξ_x) * (ξ_y^2 + ξ_y) / 4 + i == 4 && return (ξ_x^2 - ξ_x) * (ξ_y^2 + ξ_y) / 4 + i == 5 && return (1 - ξ_x^2) * (ξ_y^2 - ξ_y) / 2 + i == 6 && return (ξ_x^2 + ξ_x) * (1 - ξ_y^2) / 2 + i == 7 && return (1 - ξ_x^2) * (ξ_y^2 + ξ_y) / 2 + i == 8 && return (ξ_x^2 - ξ_x) * (1 - ξ_y^2) / 2 i == 9 && return (1 - ξ_x^2) * (1 - ξ_y^2) throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -663,8 +663,8 @@ end function shape_value(ip::Lagrange{RefQuadrilateral, 3}, ξ::Vec{2}, i::Int) # See https://defelement.com/elements/examples/quadrilateral-Q-3.html # Transform domain from [-1, 1] × [-1, 1] to [0, 1] × [0, 1] - ξ_x = ξ[1]*0.5 + 0.5 - ξ_y = ξ[2]*0.5 + 0.5 + ξ_x = (ξ[1]+1)/2 + ξ_y = (ξ[2]+1)/2 i == 1 && return (81*ξ_x^3*ξ_y^3)/4 - (81*ξ_x^3*ξ_y^2)/2 + (99*ξ_x^3*ξ_y)/4 - (9*ξ_x^3)/2 - (81*ξ_x^2*ξ_y^3)/2 + (81*ξ_x^2*ξ_y^2) - (99*ξ_x^2*ξ_y)/2 + (9*ξ_x^2) + (99*ξ_x*ξ_y^3)/4 - (99*ξ_x*ξ_y^2)/2 + (121*ξ_x*ξ_y)/4 - (11*ξ_x)/2 - (9*ξ_y^3)/2 + 9*ξ_y^2 - (11*ξ_y)/2 + 1 i == 2 && return (ξ_x*( - 81*ξ_x^2*ξ_y^3 + 162*ξ_x^2*ξ_y^2 - 99*ξ_x^2*ξ_y + 18*ξ_x^2 + 81*ξ_x*ξ_y^3 - 162*ξ_x*ξ_y^2 + 99*ξ_x*ξ_y - 18*ξ_x - 18*ξ_y^3 + 36*ξ_y^2 - 22*ξ_y + 4))/4 i == 4 && return (ξ_y*( - 81*ξ_x^3*ξ_y^2 + 81*ξ_x^3*ξ_y - 18*ξ_x^3 + 162*ξ_x^2*ξ_y^2 - 162*ξ_x^2*ξ_y + 36*ξ_x^2 - 99*ξ_x*ξ_y^2 + 99*ξ_x*ξ_y - 22*ξ_x + 18*ξ_y^2 - 18*ξ_y + 4))/4 @@ -702,7 +702,7 @@ function shape_value(ip::Lagrange{RefTriangle, 1}, ξ::Vec{2}, i::Int) ξ_y = ξ[2] i == 1 && return ξ_x i == 2 && return ξ_y - i == 3 && return 1. - ξ_x - ξ_y + i == 3 && return 1 - ξ_x - ξ_y throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -726,7 +726,7 @@ end function shape_value(ip::Lagrange{RefTriangle, 2}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] - γ = 1. - ξ_x - ξ_y + γ = 1 - ξ_x - ξ_y i == 1 && return ξ_x * (2ξ_x - 1) i == 2 && return ξ_y * (2ξ_y - 1) i == 3 && return γ * (2γ - 1) @@ -850,7 +850,7 @@ function shape_value(ip::Lagrange{RefTetrahedron, 1}, ξ::Vec{3}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] ξ_z = ξ[3] - i == 1 && return 1.0 - ξ_x - ξ_y - ξ_z + i == 1 && return 1 - ξ_x - ξ_y - ξ_z i == 2 && return ξ_x i == 3 && return ξ_y i == 4 && return ξ_z @@ -921,14 +921,14 @@ function shape_value(ip::Lagrange{RefHexahedron, 1}, ξ::Vec{3}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] ξ_z = ξ[3] - i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) - i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) - i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) - i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) - i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) - i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) - i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) - i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) + i == 1 && return (1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) / 8 + i == 2 && return (1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) / 8 + i == 3 && return (1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) / 8 + i == 4 && return (1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) / 8 + i == 5 && return (1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) / 8 + i == 6 && return (1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) / 8 + i == 7 && return (1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) / 8 + i == 8 && return (1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) / 8 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -1007,9 +1007,9 @@ end function shape_value(ip::Lagrange{RefHexahedron, 2}, ξ::Vec{3, T}, i::Int) where {T} # Some local helpers. - @inline φ₁(x::T) = -0.5*x*(1-x) + @inline φ₁(x::T) = -x*(1-x)/2 @inline φ₂(x::T) = (1+x)*(1-x) - @inline φ₃(x::T) = 0.5*x*(1+x) + @inline φ₃(x::T) = x*(1+x)/2 (ξ_x, ξ_y, ξ_z) = ξ # vertices i == 1 && return φ₁(ξ_x) * φ₁(ξ_y) * φ₁(ξ_z) @@ -1361,14 +1361,14 @@ end function shape_value(ip::Serendipity{RefQuadrilateral,2}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] - i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25(-ξ_x - ξ_y - 1) - i == 2 && return (1 + ξ_x) * (1 - ξ_y) * 0.25( ξ_x - ξ_y - 1) - i == 3 && return (1 + ξ_x) * (1 + ξ_y) * 0.25( ξ_x + ξ_y - 1) - i == 4 && return (1 - ξ_x) * (1 + ξ_y) * 0.25(-ξ_x + ξ_y - 1) - i == 5 && return 0.5(1 - ξ_x * ξ_x) * (1 - ξ_y) - i == 6 && return 0.5(1 + ξ_x) * (1 - ξ_y * ξ_y) - i == 7 && return 0.5(1 - ξ_x * ξ_x) * (1 + ξ_y) - i == 8 && return 0.5(1 - ξ_x) * (1 - ξ_y * ξ_y) + i == 1 && return (1 - ξ_x) * (1 - ξ_y) * (-ξ_x - ξ_y - 1) / 4 + i == 2 && return (1 + ξ_x) * (1 - ξ_y) * ( ξ_x - ξ_y - 1) / 4 + i == 3 && return (1 + ξ_x) * (1 + ξ_y) * ( ξ_x + ξ_y - 1) / 4 + i == 4 && return (1 - ξ_x) * (1 + ξ_y) * (-ξ_x + ξ_y - 1) / 4 + i == 5 && return (1 - ξ_x * ξ_x) * (1 - ξ_y) / 2 + i == 6 && return (1 + ξ_x) * (1 - ξ_y * ξ_y) / 2 + i == 7 && return (1 - ξ_x * ξ_x) * (1 + ξ_y) / 2 + i == 8 && return (1 - ξ_x) * (1 - ξ_y * ξ_y) / 2 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -1429,30 +1429,31 @@ function reference_coordinates(::Serendipity{RefHexahedron,2}) Vec{3, Float64}((-1.0, 1.0, 0.0)),] end -function shape_value(ip::Serendipity{RefHexahedron, 2}, ξ::Vec{3}, i::Int) +# Inlined to resolve the recursion properly +@inline function shape_value(ip::Serendipity{RefHexahedron, 2}, ξ::Vec{3}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] ξ_z = ξ[3] - i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 12) + shape_value(ip, ξ, 9) + shape_value(ip, ξ, 17)) - i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 9) + shape_value(ip, ξ, 10) + shape_value(ip, ξ, 18)) - i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 10) + shape_value(ip, ξ, 11) + shape_value(ip, ξ, 19)) - i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 11) + shape_value(ip, ξ, 12) + shape_value(ip, ξ, 20)) - i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 16) + shape_value(ip, ξ, 13) + shape_value(ip, ξ, 17)) - i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 13) + shape_value(ip, ξ, 14) + shape_value(ip, ξ, 18)) - i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 14) + shape_value(ip, ξ, 15) + shape_value(ip, ξ, 19)) - i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 15) + shape_value(ip, ξ, 16) + shape_value(ip, ξ, 20)) - i == 9 && return 0.25(1 - ξ_x^2) * (1 - ξ_y) * (1 - ξ_z) - i == 10 && return 0.25(1 + ξ_x) * (1 - ξ_y^2) * (1 - ξ_z) - i == 11 && return 0.25(1 - ξ_x^2) * (1 + ξ_y) * (1 - ξ_z) - i == 12 && return 0.25(1 - ξ_x) * (1 - ξ_y^2) * (1 - ξ_z) - i == 13 && return 0.25(1 - ξ_x^2) * (1 - ξ_y) * (1 + ξ_z) - i == 14 && return 0.25(1 + ξ_x) * (1 - ξ_y^2) * (1 + ξ_z) - i == 15 && return 0.25(1 - ξ_x^2) * (1 + ξ_y) * (1 + ξ_z) - i == 16 && return 0.25(1 - ξ_x) * (1 - ξ_y^2) * (1 + ξ_z) - i == 17 && return 0.25(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z^2) - i == 18 && return 0.25(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z^2) - i == 19 && return 0.25(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z^2) - i == 20 && return 0.25(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z^2) + i == 1 && return (1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) / 8 - (shape_value(ip, ξ, 12) + shape_value(ip, ξ, 9) + shape_value(ip, ξ, 17)) / 2 + i == 2 && return (1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) / 8 - (shape_value(ip, ξ, 9) + shape_value(ip, ξ, 10) + shape_value(ip, ξ, 18)) / 2 + i == 3 && return (1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) / 8 - (shape_value(ip, ξ, 10) + shape_value(ip, ξ, 11) + shape_value(ip, ξ, 19)) / 2 + i == 4 && return (1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) / 8 - (shape_value(ip, ξ, 11) + shape_value(ip, ξ, 12) + shape_value(ip, ξ, 20)) / 2 + i == 5 && return (1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) / 8 - (shape_value(ip, ξ, 16) + shape_value(ip, ξ, 13) + shape_value(ip, ξ, 17)) / 2 + i == 6 && return (1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) / 8 - (shape_value(ip, ξ, 13) + shape_value(ip, ξ, 14) + shape_value(ip, ξ, 18)) / 2 + i == 7 && return (1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) / 8 - (shape_value(ip, ξ, 14) + shape_value(ip, ξ, 15) + shape_value(ip, ξ, 19)) / 2 + i == 8 && return (1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) / 8 - (shape_value(ip, ξ, 15) + shape_value(ip, ξ, 16) + shape_value(ip, ξ, 20)) / 2 + i == 9 && return (1 - ξ_x^2) * (1 - ξ_y) * (1 - ξ_z) / 4 + i == 10 && return (1 + ξ_x) * (1 - ξ_y^2) * (1 - ξ_z) / 4 + i == 11 && return (1 - ξ_x^2) * (1 + ξ_y) * (1 - ξ_z) / 4 + i == 12 && return (1 - ξ_x) * (1 - ξ_y^2) * (1 - ξ_z) / 4 + i == 13 && return (1 - ξ_x^2) * (1 - ξ_y) * (1 + ξ_z) / 4 + i == 14 && return (1 + ξ_x) * (1 - ξ_y^2) * (1 + ξ_z) / 4 + i == 15 && return (1 - ξ_x^2) * (1 + ξ_y) * (1 + ξ_z) / 4 + i == 16 && return (1 - ξ_x) * (1 - ξ_y^2) * (1 + ξ_z) / 4 + i == 17 && return (1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z^2) / 4 + i == 18 && return (1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z^2) / 4 + i == 19 && return (1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z^2) / 4 + i == 20 && return (1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z^2) / 4 throw(ArgumentError("no shape function $i for interpolation $ip")) end @@ -1486,12 +1487,12 @@ function reference_coordinates(::CrouzeixRaviart) Vec{2, Float64}((0.5, 0.0))] end -function shape_value(ip::CrouzeixRaviart, ξ::Vec{2}, i::Int) +function shape_value(ip::CrouzeixRaviart{RefTriangle, 1}, ξ::Vec{2}, i::Int) ξ_x = ξ[1] ξ_y = ξ[2] - i == 1 && return 2*ξ_x + 2*ξ_y - 1.0 - i == 2 && return 1.0 - 2*ξ_x - i == 3 && return 1.0 - 2*ξ_y + i == 1 && return 2*ξ_x + 2*ξ_y - 1 + i == 2 && return 1 - 2*ξ_x + i == 3 && return 1 - 2*ξ_y throw(ArgumentError("no shape function $i for interpolation $ip")) end diff --git a/src/iterators.jl b/src/iterators.jl index ddab91900b..b73a2bbd2e 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -67,8 +67,8 @@ function CellCache(dh::DofHandler{dim}, flags::UpdateFlags=UpdateFlags()) where end function CellCache(sdh::SubDofHandler, flags::UpdateFlags=UpdateFlags()) - dim = getdim(sdh.dh.grid) - CellCache(flags, sdh.dh.grid, ScalarWrapper(-1), Int[], Vec{dim,Float64}[], sdh, Int[]) + Tv = get_coordinate_type(sdh.dh.grid) + CellCache(flags, sdh.dh.grid, ScalarWrapper(-1), Int[], Tv[], sdh, Int[]) end function reinit!(cc::CellCache, i::Int) diff --git a/test/runtests.jl b/test/runtests.jl index e31028893a..72f7d71e53 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,17 @@ if HAS_EXTENSIONS && MODULE_CAN_BE_TYPE_PARAMETER import Metis end +const RUN_JET_TESTS = VERSION >= v"1.9" + +if RUN_JET_TESTS + using JET: @test_call +else + # Just eat the macro on incompatible versions + macro test_call(args...) + nothing + end +end + include("test_utils.jl") # Unit tests diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 052251c1f7..f414b4c6cb 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -28,6 +28,7 @@ x, n = valid_coordinates_and_normals(func_interpol) reinit!(cv, x) + @test_call reinit!(cv, x) # We test this by applying a given deformation gradient on all the nodes. # Since this is a linear deformation we should get back the exact values @@ -199,6 +200,7 @@ end ## sdim = 2, Consistency with 1D csv2 = CellValues(qr, ip, ip_base^2) reinit!(csv2, [Vec((0.0, 0.0)), Vec((1.0, 0.0))]) + @test_call skip=true reinit!(csv2, [Vec((0.0, 0.0)), Vec((1.0, 0.0))]) # External error in pinv # Test spatial interpolation @test spatial_coordinate(csv2, 1, [Vec((0.0, 0.0)), Vec((1.0, 0.0))]) == Vec{2}((0.5, 0.0)) # Test volume @@ -219,6 +221,7 @@ end ## sdim = 3, Consistency with 1D csv3 = CellValues(qr, ip, ip_base^3) reinit!(csv3, [Vec((0.0, 0.0, 0.0)), Vec((1.0, 0.0, 0.0))]) + @test_call skip=true reinit!(csv3, [Vec((0.0, 0.0, 0.0)), Vec((1.0, 0.0, 0.0))]) # External error in pinv # Test spatial interpolation @test spatial_coordinate(csv3, 1, [Vec((0.0, 0.0, 0.0)), Vec((1.0, 0.0, 0.0))]) == Vec{3}((0.5, 0.0, 0.0)) # Test volume @@ -281,7 +284,9 @@ end csv2 = CellValues(qr, ip) csv3 = CellValues(qr, ip, ip_base^3) reinit!(csv2, [Vec((-1.0,-1.0)), Vec((1.0,-1.0)), Vec((1.0,1.0)), Vec((-1.0,1.0))]) + @test_call skip=true reinit!(csv2, [Vec((-1.0,-1.0)), Vec((1.0,-1.0)), Vec((1.0,1.0)), Vec((-1.0,1.0))]) # External error in pinv reinit!(csv3, [Vec((-1.0,-1.0,0.0)), Vec((1.0,-1.0,0.0)), Vec((1.0,1.0,0.0)), Vec((-1.0,1.0,0.0))]) + @test_call skip=true reinit!(csv3, [Vec((-1.0,-1.0,0.0)), Vec((1.0,-1.0,0.0)), Vec((1.0,1.0,0.0)), Vec((-1.0,1.0,0.0))]) # External error in pinv # Test spatial interpolation @test spatial_coordinate(csv2, 1, [Vec((-1.0,-1.0)), Vec((1.0,-1.0)), Vec((1.0,1.0)), Vec((-1.0,1.0))]) == Vec{2}((0.0, 0.0)) @test spatial_coordinate(csv3, 1, [Vec((-1.0,-1.0,0.0)), Vec((1.0,-1.0,0.0)), Vec((1.0,1.0,0.0)), Vec((-1.0,1.0,0.0))]) == Vec{3}((0.0, 0.0, 0.0)) @@ -328,6 +333,8 @@ end @test Ferrite.shape_gradient_type(cv) == grad_type(Float32) @test Ferrite.geometric_interpolation(cv) == scalar_ip(geo_ip) end + x = Ferrite.reference_coordinates(fun_ip) + @test_call reinit!(cv, x) end end diff --git a/test/test_dofs.jl b/test/test_dofs.jl index 676fa1c6d6..8264cb1d1b 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -432,6 +432,7 @@ end # Full coupling (default) K = create_sparsity_pattern(dh) + @test eltype(K) == Float64 for j in 1:ndofs(dh), i in 1:ndofs(dh) @test is_stored(K, i, j) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index b67e0599a1..25c6905154 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -1,6 +1,7 @@ @testset "interpolations" begin -@testset "$interpolation" for interpolation in ( +@testset "Value Type $value_type" for value_type in (Float32, Float64) + @testset "Correctness of $interpolation" for interpolation in ( Lagrange{RefLine, 1}(), Lagrange{RefLine, 2}(), Lagrange{RefQuadrilateral, 1}(), @@ -75,11 +76,11 @@ #TODO prefer this test style after 1.6 is removed from CI # @testset let x = sample_random_point(ref_shape) # not compatible with Julia 1.6 - x = sample_random_point(ref_shape) + x = Vec{ref_dim,value_type}(sample_random_point(ref_shape)) random_point_testset = @testset "Random point test" begin # Check gradient evaluation @test vec(ForwardDiff.jacobian(f, Array(x))') ≈ - reinterpret(Float64, [shape_gradient(interpolation, x, i) for i in 1:n_basefuncs]) + reinterpret(value_type, [shape_gradient(interpolation, x, i) for i in 1:n_basefuncs]) # Check partition of unity at random point. @test sum([shape_value(interpolation, x, i) for i in 1:n_basefuncs]) ≈ 1.0 # Check if the important functions are consistent @@ -112,7 +113,7 @@ totaldofs += sum(length.(Ferrite.facedof_interior_indices(interpolation));init=0) end if ref_dim > 2 - totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) + totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) end totaldofs += length(Ferrite.celldof_interior_indices(interpolation)) @test totaldofs == n_basefuncs @@ -128,6 +129,12 @@ end @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.celldof_interior_indices(interpolation)]) + # Check for evaluation type correctness of interpolation + @testset "return type correctness dof $dof" for dof in 1:n_basefuncs + @test (@inferred shape_value(interpolation, x, dof)) isa value_type + @test (@inferred shape_gradient(interpolation, x, dof)) isa Vec{ref_dim, value_type} + end + # Check for dirac delta property of interpolation @testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs for k in 1:n_basefuncs @@ -136,7 +143,7 @@ @test N_dof ≈ 1.0 else factor = interpolation isa Lagrange{RefQuadrilateral, 3} ? 200 : 4 - @test N_dof ≈ 0.0 atol = factor * eps(Float64) + @test N_dof ≈ 0.0 atol = factor * eps(value_type) end end end @@ -183,14 +190,22 @@ end end - # VectorizedInterpolation - v_interpolation_1 = interpolation^2 - v_interpolation_2 = (d = 2; interpolation^d) - @test getnbasefunctions(v_interpolation_1) == getnbasefunctions(v_interpolation_2) == - getnbasefunctions(interpolation) * 2 - # pretty printing - @test repr("text/plain", v_interpolation_1) == repr(v_interpolation_1.ip) * "^2" -end + # VectorizedInterpolation + v_interpolation_1 = interpolation^2 + v_interpolation_2 = (d = 2; interpolation^d) + @test getnbasefunctions(v_interpolation_1) == + getnbasefunctions(v_interpolation_2) == + getnbasefunctions(interpolation) * 2 + # pretty printing + @test repr("text/plain", v_interpolation_1) == repr(v_interpolation_1.ip) * "^2" + + # Check for evaluation type correctness of vectorized interpolation + v_interpolation_3 = interpolation^ref_dim + @testset "vectorized case of return type correctness of dof $dof" for dof in 1:n_basefuncs + @test @inferred(shape_value(v_interpolation_1, x, dof)) isa Vec{2, value_type} + @test @inferred(shape_gradient(v_interpolation_3, x, dof)) isa Tensor{2, ref_dim, value_type} + end + end # correctness testset @test Ferrite.reference_coordinates(DiscontinuousLagrange{RefTriangle,0}()) ≈ [Vec{2,Float64}((1/3,1/3))] @test Ferrite.reference_coordinates(DiscontinuousLagrange{RefQuadrilateral,0}()) ≈ [Vec{2,Float64}((0,0))] @@ -209,3 +224,5 @@ end @test Ferrite.is_discontinuous(d_ip) == true @test Ferrite.is_discontinuous(d_ip_t) == true end + +end # testset