From 4f8f85c89568ca4b24d8122b449475b25394a36b Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 9 Nov 2025 17:30:07 +0200 Subject: [PATCH] chore(basis): Regenerate basis functions for Lagrange{T,P} architecture MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Generated by: julia --project=. src/basis/lagrange_generator.jl Changes: - All 15 element types now use Lagrange{T,P} parametric type - Functions: get_reference_element_coordinates(), eval_basis!(), eval_dbasis!() - Reference coordinates now return tuples (zero-allocation) - Removed old Seg2Basis, Tri3Basis, Quad4Basis, etc. struct definitions - All methods work with both Type{Lagrange{T,P}} and Lagrange{T,P} instances Validated: - Triangle: Kronecker delta property holds (N_i(x_j) = δ_ij) - Quadrilateral, Tetrahedron, Hexahedron: First node evaluates to (1,0,0,...) - Derivatives: Correct gradients at reference coordinates --- src/basis/lagrange_generated.jl | 489 ++++++++++++++++---------------- 1 file changed, 246 insertions(+), 243 deletions(-) diff --git a/src/basis/lagrange_generated.jl b/src/basis/lagrange_generated.jl index 635319e..85ddcbe 100644 --- a/src/basis/lagrange_generated.jl +++ b/src/basis/lagrange_generated.jl @@ -20,428 +20,431 @@ # Generator: # src/basis/lagrange_generator.jl (symbolic engine) # -# Generated: 2025-11-09 07:28:58 +# Generated: 2025-11-09 17:25:56 # ============================================================================ -# Export all basis types -export Seg2Basis, Seg3Basis, Tri3Basis, Tri6Basis, Quad4Basis, Quad8Basis, Quad9Basis, Tet4Basis, Tet10Basis, Hex8Basis, Hex20Basis, Hex27Basis, Pyr5Basis, Wedge6Basis, Wedge15Basis +# This file generates methods for Lagrange{T,P} where: +# T = topology type (Segment, Triangle, Quadrilateral, etc.) +# P = polynomial degree (1, 2, 3, ...) +# +# Example: Lagrange{Triangle, 2} is a quadratic triangular element # ────────────────────────────────────────────────────────────────────────────── -# Seg2: 2-node linear segment element +# Lagrange{Segment, 1}: 2-node linear segment element +# (Old name: Seg2) # ────────────────────────────────────────────────────────────────────────────── - struct Seg2Basis <: AbstractBasis{1} - end - Base.@pure function Base.size(::Type{Seg2Basis}) - return (1, 2) - end - function Base.size(::Type{Seg2Basis}, j::Int) - j == 1 && return 1 - j == 2 && return 2 + function get_reference_element_coordinates(::Type{Lagrange{Segment, 1}}) + return (Vec{1, Float64}(tuple(-1.0)), Vec{1, Float64}(tuple(1.0))) end - Base.@pure function Base.length(::Type{Seg2Basis}) - return 2 - end - function get_reference_element_coordinates(::Type{Seg2Basis}) - return Vec{1, Float64}[[-1.0], [1.0]] + function get_reference_element_coordinates(::Lagrange{Segment, 1}) + return (Vec{1, Float64}(tuple(-1.0)), Vec{1, Float64}(tuple(1.0))) end - @inline function eval_basis!(::Type{Seg2Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Segment, 1}}, ::Type{T}, xi::Vec) where T (u,) = xi @inbounds return (0.5 + -0.5u, 0.5 + 0.5u) end - @inline function eval_dbasis!(::Type{Seg2Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Segment, 1}, ::Type{T}, xi::Vec) where T + (u,) = xi + @inbounds return (0.5 + -0.5u, 0.5 + 0.5u) + end + @inline function eval_dbasis!(::Type{Lagrange{Segment, 1}}, xi::Vec) + (u,) = xi + @inbounds return (Vec(float.(tuple(-0.5))), Vec(float.(tuple(0.5)))) + end + @inline function eval_dbasis!(::Lagrange{Segment, 1}, xi::Vec) (u,) = xi @inbounds return (Vec(float.(tuple(-0.5))), Vec(float.(tuple(0.5)))) end # ────────────────────────────────────────────────────────────────────────────── -# Seg3: 3-node quadratic segment element +# Lagrange{Segment, 2}: 3-node quadratic segment element +# (Old name: Seg3) # ────────────────────────────────────────────────────────────────────────────── - struct Seg3Basis <: AbstractBasis{1} - end - Base.@pure function Base.size(::Type{Seg3Basis}) - return (1, 3) - end - function Base.size(::Type{Seg3Basis}, j::Int) - j == 1 && return 1 - j == 2 && return 3 + function get_reference_element_coordinates(::Type{Lagrange{Segment, 2}}) + return (Vec{1, Float64}(tuple(-1.0)), Vec{1, Float64}(tuple(1.0)), Vec{1, Float64}(tuple(0.0))) end - Base.@pure function Base.length(::Type{Seg3Basis}) - return 3 - end - function get_reference_element_coordinates(::Type{Seg3Basis}) - return Vec{1, Float64}[[-1.0], [1.0], [0.0]] + function get_reference_element_coordinates(::Lagrange{Segment, 2}) + return (Vec{1, Float64}(tuple(-1.0)), Vec{1, Float64}(tuple(1.0)), Vec{1, Float64}(tuple(0.0))) end - @inline function eval_basis!(::Type{Seg3Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Segment, 2}}, ::Type{T}, xi::Vec) where T (u,) = xi @inbounds return (-0.5u + 0.5 * u ^ 2, 0.5u + 0.5 * u ^ 2, 1 + -1.0 * u ^ 2) end - @inline function eval_dbasis!(::Type{Seg3Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Segment, 2}, ::Type{T}, xi::Vec) where T + (u,) = xi + @inbounds return (-0.5u + 0.5 * u ^ 2, 0.5u + 0.5 * u ^ 2, 1 + -1.0 * u ^ 2) + end + @inline function eval_dbasis!(::Type{Lagrange{Segment, 2}}, xi::Vec) + (u,) = xi + @inbounds return (Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)))))) + end + @inline function eval_dbasis!(::Lagrange{Segment, 2}, xi::Vec) (u,) = xi @inbounds return (Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)))))) end # ────────────────────────────────────────────────────────────────────────────── -# Tri3: 3-node linear triangular element +# Lagrange{Triangle, 1}: 3-node linear triangular element +# (Old name: Tri3) # ────────────────────────────────────────────────────────────────────────────── - struct Tri3Basis <: AbstractBasis{2} - end - Base.@pure function Base.size(::Type{Tri3Basis}) - return (2, 3) - end - function Base.size(::Type{Tri3Basis}, j::Int) - j == 1 && return 2 - j == 2 && return 3 + function get_reference_element_coordinates(::Type{Lagrange{Triangle, 1}}) + return (Vec{2, Float64}(tuple(0.0, 0.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0))) end - Base.@pure function Base.length(::Type{Tri3Basis}) - return 3 - end - function get_reference_element_coordinates(::Type{Tri3Basis}) - return Vec{2, Float64}[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Triangle, 1}) + return (Vec{2, Float64}(tuple(0.0, 0.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0))) end - @inline function eval_basis!(::Type{Tri3Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Triangle, 1}}, ::Type{T}, xi::Vec) where T (u, v) = xi @inbounds return (1 + -1.0u + -1.0v, +u, +v) end - @inline function eval_dbasis!(::Type{Tri3Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Triangle, 1}, ::Type{T}, xi::Vec) where T + (u, v) = xi + @inbounds return (1 + -1.0u + -1.0v, +u, +v) + end + @inline function eval_dbasis!(::Type{Lagrange{Triangle, 1}}, xi::Vec) + (u, v) = xi + @inbounds return (Vec(float.(tuple(-1.0, -1.0))), Vec(float.(tuple(1, 0))), Vec(float.(tuple(0, 1)))) + end + @inline function eval_dbasis!(::Lagrange{Triangle, 1}, xi::Vec) (u, v) = xi @inbounds return (Vec(float.(tuple(-1.0, -1.0))), Vec(float.(tuple(1, 0))), Vec(float.(tuple(0, 1)))) end # ────────────────────────────────────────────────────────────────────────────── -# Tri6: 6-node quadratic triangular element +# Lagrange{Triangle, 2}: 6-node quadratic triangular element +# (Old name: Tri6) # ────────────────────────────────────────────────────────────────────────────── - struct Tri6Basis <: AbstractBasis{2} - end - Base.@pure function Base.size(::Type{Tri6Basis}) - return (2, 6) - end - function Base.size(::Type{Tri6Basis}, j::Int) - j == 1 && return 2 - j == 2 && return 6 + function get_reference_element_coordinates(::Type{Lagrange{Triangle, 2}}) + return (Vec{2, Float64}(tuple(0.0, 0.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(0.5, 0.0)), Vec{2, Float64}(tuple(0.5, 0.5)), Vec{2, Float64}(tuple(0.0, 0.5))) end - Base.@pure function Base.length(::Type{Tri6Basis}) - return 6 - end - function get_reference_element_coordinates(::Type{Tri6Basis}) - return Vec{2, Float64}[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.5], [0.0, 0.5]] + function get_reference_element_coordinates(::Lagrange{Triangle, 2}) + return (Vec{2, Float64}(tuple(0.0, 0.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(0.5, 0.0)), Vec{2, Float64}(tuple(0.5, 0.5)), Vec{2, Float64}(tuple(0.0, 0.5))) end - @inline function eval_basis!(::Type{Tri6Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Triangle, 2}}, ::Type{T}, xi::Vec) where T (u, v) = xi @inbounds return (1 + -3.0u + -3.0v + 2.0 * u ^ 2 + 4.0 * (u * v) + 2.0 * v ^ 2, -1.0u + 2.0 * u ^ 2, -1.0v + 2.0 * v ^ 2, 4.0u + -4.0 * u ^ 2 + -4.0 * (u * v), +(4.0 * (u * v)), 4.0v + -4.0 * (u * v) + -4.0 * v ^ 2) end - @inline function eval_dbasis!(::Type{Tri6Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Triangle, 2}, ::Type{T}, xi::Vec) where T + (u, v) = xi + @inbounds return (1 + -3.0u + -3.0v + 2.0 * u ^ 2 + 4.0 * (u * v) + 2.0 * v ^ 2, -1.0u + 2.0 * u ^ 2, -1.0v + 2.0 * v ^ 2, 4.0u + -4.0 * u ^ 2 + -4.0 * (u * v), +(4.0 * (u * v)), 4.0v + -4.0 * (u * v) + -4.0 * v ^ 2) + end + @inline function eval_dbasis!(::Type{Lagrange{Triangle, 2}}, xi::Vec) + (u, v) = xi + @inbounds return (Vec(float.(tuple(-3.0 + 2.0 * (2 * u ^ (2 - 1)) + 4.0v, -3.0 + 4.0u + 2.0 * (2 * v ^ (2 - 1))))), Vec(float.(tuple(-1.0 + 2.0 * (2 * u ^ (2 - 1)), 0))), Vec(float.(tuple(0, -1.0 + 2.0 * (2 * v ^ (2 - 1))))), Vec(float.(tuple(4.0 + -4.0 * (2 * u ^ (2 - 1)) + -4.0v, -4.0u))), Vec(float.(tuple(4.0v, 4.0u))), Vec(float.(tuple(-4.0v, 4.0 + -4.0u + -4.0 * (2 * v ^ (2 - 1)))))) + end + @inline function eval_dbasis!(::Lagrange{Triangle, 2}, xi::Vec) (u, v) = xi @inbounds return (Vec(float.(tuple(-3.0 + 2.0 * (2 * u ^ (2 - 1)) + 4.0v, -3.0 + 4.0u + 2.0 * (2 * v ^ (2 - 1))))), Vec(float.(tuple(-1.0 + 2.0 * (2 * u ^ (2 - 1)), 0))), Vec(float.(tuple(0, -1.0 + 2.0 * (2 * v ^ (2 - 1))))), Vec(float.(tuple(4.0 + -4.0 * (2 * u ^ (2 - 1)) + -4.0v, -4.0u))), Vec(float.(tuple(4.0v, 4.0u))), Vec(float.(tuple(-4.0v, 4.0 + -4.0u + -4.0 * (2 * v ^ (2 - 1)))))) end # ────────────────────────────────────────────────────────────────────────────── -# Quad4: 4-node bilinear quadrilateral element +# Lagrange{Quadrilateral, 1}: 4-node bilinear quadrilateral element +# (Old name: Quad4) # ────────────────────────────────────────────────────────────────────────────── - struct Quad4Basis <: AbstractBasis{2} - end - Base.@pure function Base.size(::Type{Quad4Basis}) - return (2, 4) - end - function Base.size(::Type{Quad4Basis}, j::Int) - j == 1 && return 2 - j == 2 && return 4 + function get_reference_element_coordinates(::Type{Lagrange{Quadrilateral, 1}}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0))) end - Base.@pure function Base.length(::Type{Quad4Basis}) - return 4 - end - function get_reference_element_coordinates(::Type{Quad4Basis}) - return Vec{2, Float64}[[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Quadrilateral, 1}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0))) end - @inline function eval_basis!(::Type{Quad4Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Quadrilateral, 1}}, ::Type{T}, xi::Vec) where T (u, v) = xi @inbounds return (0.25 + -0.25u + -0.25v + 0.25 * (u * v), 0.25 + 0.25u + -0.25v + -0.25 * (u * v), 0.25 + 0.25u + 0.25v + 0.25 * (u * v), 0.25 + -0.25u + 0.25v + -0.25 * (u * v)) end - @inline function eval_dbasis!(::Type{Quad4Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Quadrilateral, 1}, ::Type{T}, xi::Vec) where T + (u, v) = xi + @inbounds return (0.25 + -0.25u + -0.25v + 0.25 * (u * v), 0.25 + 0.25u + -0.25v + -0.25 * (u * v), 0.25 + 0.25u + 0.25v + 0.25 * (u * v), 0.25 + -0.25u + 0.25v + -0.25 * (u * v)) + end + @inline function eval_dbasis!(::Type{Lagrange{Quadrilateral, 1}}, xi::Vec) + (u, v) = xi + @inbounds return (Vec(float.(tuple(-0.25 + 0.25v, -0.25 + 0.25u))), Vec(float.(tuple(0.25 + -0.25v, -0.25 + -0.25u))), Vec(float.(tuple(0.25 + 0.25v, 0.25 + 0.25u))), Vec(float.(tuple(-0.25 + -0.25v, 0.25 + -0.25u)))) + end + @inline function eval_dbasis!(::Lagrange{Quadrilateral, 1}, xi::Vec) (u, v) = xi @inbounds return (Vec(float.(tuple(-0.25 + 0.25v, -0.25 + 0.25u))), Vec(float.(tuple(0.25 + -0.25v, -0.25 + -0.25u))), Vec(float.(tuple(0.25 + 0.25v, 0.25 + 0.25u))), Vec(float.(tuple(-0.25 + -0.25v, 0.25 + -0.25u)))) end # ────────────────────────────────────────────────────────────────────────────── -# Quad8: 8-node serendipity quadrilateral element +# Lagrange{Quadrilateral, 2}: 8-node serendipity quadrilateral element +# (Old name: Quad8) # ────────────────────────────────────────────────────────────────────────────── - struct Quad8Basis <: AbstractBasis{2} - end - Base.@pure function Base.size(::Type{Quad8Basis}) - return (2, 8) - end - function Base.size(::Type{Quad8Basis}, j::Int) - j == 1 && return 2 - j == 2 && return 8 + function get_reference_element_coordinates(::Type{Lagrange{Quadrilateral, 2}}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0)), Vec{2, Float64}(tuple(0.0, -1.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 0.0))) end - Base.@pure function Base.length(::Type{Quad8Basis}) - return 8 - end - function get_reference_element_coordinates(::Type{Quad8Basis}) - return Vec{2, Float64}[[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [0.0, -1.0], [1.0, 0.0], [0.0, 1.0], [-1.0, 0.0]] + function get_reference_element_coordinates(::Lagrange{Quadrilateral, 2}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0)), Vec{2, Float64}(tuple(0.0, -1.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 0.0))) end - @inline function eval_basis!(::Type{Quad8Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Quadrilateral, 2}}, ::Type{T}, xi::Vec) where T (u, v) = xi @inbounds return (-0.25 + 0.25 * u ^ 2 + 0.25 * (u * v) + 0.25 * v ^ 2 + -0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + -0.25 * (u * v) + 0.25 * v ^ 2 + -0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + 0.25 * (u * v) + 0.25 * v ^ 2 + 0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + -0.25 * (u * v) + 0.25 * v ^ 2 + 0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2), 0.5 + -0.5v + -0.5 * u ^ 2 + 0.5 * (u ^ 2 * v), 0.5 + 0.5u + -0.5 * v ^ 2 + -0.5 * (u * v ^ 2), 0.5 + 0.5v + -0.5 * u ^ 2 + -0.5 * (u ^ 2 * v), 0.5 + -0.5u + -0.5 * v ^ 2 + 0.5 * (u * v ^ 2)) end - @inline function eval_dbasis!(::Type{Quad8Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Quadrilateral, 2}, ::Type{T}, xi::Vec) where T + (u, v) = xi + @inbounds return (-0.25 + 0.25 * u ^ 2 + 0.25 * (u * v) + 0.25 * v ^ 2 + -0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + -0.25 * (u * v) + 0.25 * v ^ 2 + -0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + 0.25 * (u * v) + 0.25 * v ^ 2 + 0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2), -0.25 + 0.25 * u ^ 2 + -0.25 * (u * v) + 0.25 * v ^ 2 + 0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2), 0.5 + -0.5v + -0.5 * u ^ 2 + 0.5 * (u ^ 2 * v), 0.5 + 0.5u + -0.5 * v ^ 2 + -0.5 * (u * v ^ 2), 0.5 + 0.5v + -0.5 * u ^ 2 + -0.5 * (u ^ 2 * v), 0.5 + -0.5u + -0.5 * v ^ 2 + 0.5 * (u * v ^ 2)) + end + @inline function eval_dbasis!(::Type{Lagrange{Quadrilateral, 2}}, xi::Vec) + (u, v) = xi + @inbounds return (Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + 0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2, 0.25u + 0.25 * (2 * v ^ (2 - 1)) + -0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + -0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2, -0.25u + 0.25 * (2 * v ^ (2 - 1)) + -0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + 0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2, 0.25u + 0.25 * (2 * v ^ (2 - 1)) + 0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + -0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2, -0.25u + 0.25 * (2 * v ^ (2 - 1)) + 0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * (2 * u ^ (2 - 1)) + 0.5 * ((2 * u ^ (2 - 1)) * v), -0.5 + 0.5 * u ^ 2))), Vec(float.(tuple(0.5 + -0.5 * v ^ 2, -0.5 * (2 * v ^ (2 - 1)) + -0.5 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * (2 * u ^ (2 - 1)) + -0.5 * ((2 * u ^ (2 - 1)) * v), 0.5 + -0.5 * u ^ 2))), Vec(float.(tuple(-0.5 + 0.5 * v ^ 2, -0.5 * (2 * v ^ (2 - 1)) + 0.5 * (u * (2 * v ^ (2 - 1))))))) + end + @inline function eval_dbasis!(::Lagrange{Quadrilateral, 2}, xi::Vec) (u, v) = xi @inbounds return (Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + 0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2, 0.25u + 0.25 * (2 * v ^ (2 - 1)) + -0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + -0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2, -0.25u + 0.25 * (2 * v ^ (2 - 1)) + -0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + 0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2, 0.25u + 0.25 * (2 * v ^ (2 - 1)) + 0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25 * (2 * u ^ (2 - 1)) + -0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2, -0.25u + 0.25 * (2 * v ^ (2 - 1)) + 0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * (2 * u ^ (2 - 1)) + 0.5 * ((2 * u ^ (2 - 1)) * v), -0.5 + 0.5 * u ^ 2))), Vec(float.(tuple(0.5 + -0.5 * v ^ 2, -0.5 * (2 * v ^ (2 - 1)) + -0.5 * (u * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * (2 * u ^ (2 - 1)) + -0.5 * ((2 * u ^ (2 - 1)) * v), 0.5 + -0.5 * u ^ 2))), Vec(float.(tuple(-0.5 + 0.5 * v ^ 2, -0.5 * (2 * v ^ (2 - 1)) + 0.5 * (u * (2 * v ^ (2 - 1))))))) end # ────────────────────────────────────────────────────────────────────────────── -# Quad9: 9-node biquadratic quadrilateral element +# Lagrange{Quadrilateral, 2}: 9-node biquadratic quadrilateral element +# (Old name: Quad9) # ────────────────────────────────────────────────────────────────────────────── - struct Quad9Basis <: AbstractBasis{2} - end - Base.@pure function Base.size(::Type{Quad9Basis}) - return (2, 9) - end - function Base.size(::Type{Quad9Basis}, j::Int) - j == 1 && return 2 - j == 2 && return 9 + function get_reference_element_coordinates(::Type{Lagrange{Quadrilateral, 2}}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0)), Vec{2, Float64}(tuple(0.0, -1.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 0.0))) end - Base.@pure function Base.length(::Type{Quad9Basis}) - return 9 - end - function get_reference_element_coordinates(::Type{Quad9Basis}) - return Vec{2, Float64}[[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [0.0, -1.0], [1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, 0.0]] + function get_reference_element_coordinates(::Lagrange{Quadrilateral, 2}) + return (Vec{2, Float64}(tuple(-1.0, -1.0)), Vec{2, Float64}(tuple(1.0, -1.0)), Vec{2, Float64}(tuple(1.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 1.0)), Vec{2, Float64}(tuple(0.0, -1.0)), Vec{2, Float64}(tuple(1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 1.0)), Vec{2, Float64}(tuple(-1.0, 0.0)), Vec{2, Float64}(tuple(0.0, 0.0))) end - @inline function eval_basis!(::Type{Quad9Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Quadrilateral, 2}}, ::Type{T}, xi::Vec) where T (u, v) = xi @inbounds return (0.25 * (u * v) + -0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.25 * (u * v) + -0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), 0.25 * (u * v) + 0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.25 * (u * v) + 0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.5v + 0.5 * v ^ 2 + 0.5 * (u ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2), 0.5u + 0.5 * u ^ 2 + -0.5 * (u * v ^ 2) + -0.5 * (u ^ 2 * v ^ 2), 0.5v + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2), -0.5u + 0.5 * u ^ 2 + 0.5 * (u * v ^ 2) + -0.5 * (u ^ 2 * v ^ 2), 1 + -1.0 * u ^ 2 + -1.0 * v ^ 2 + u ^ 2 * v ^ 2) end - @inline function eval_dbasis!(::Type{Quad9Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Quadrilateral, 2}, ::Type{T}, xi::Vec) where T + (u, v) = xi + @inbounds return (0.25 * (u * v) + -0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.25 * (u * v) + -0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), 0.25 * (u * v) + 0.25 * (u ^ 2 * v) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.25 * (u * v) + 0.25 * (u ^ 2 * v) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2), -0.5v + 0.5 * v ^ 2 + 0.5 * (u ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2), 0.5u + 0.5 * u ^ 2 + -0.5 * (u * v ^ 2) + -0.5 * (u ^ 2 * v ^ 2), 0.5v + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2), -0.5u + 0.5 * u ^ 2 + 0.5 * (u * v ^ 2) + -0.5 * (u ^ 2 * v ^ 2), 1 + -1.0 * u ^ 2 + -1.0 * v ^ 2 + u ^ 2 * v ^ 2) + end + @inline function eval_dbasis!(::Type{Lagrange{Quadrilateral, 2}}, xi::Vec) + (u, v) = xi + @inbounds return (Vec(float.(tuple(0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.5 + 0.5 * (2 * v ^ (2 - 1)) + 0.5 * u ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1)) + -0.5 * v ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.5 * (u * (2 * v ^ (2 - 1))) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.5 + 0.5 * (2 * v ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1)) + 0.5 * v ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.5 * (u * (2 * v ^ (2 - 1))) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)) + (2 * u ^ (2 - 1)) * v ^ 2, -1.0 * (2 * v ^ (2 - 1)) + u ^ 2 * (2 * v ^ (2 - 1)))))) + end + @inline function eval_dbasis!(::Lagrange{Quadrilateral, 2}, xi::Vec) (u, v) = xi @inbounds return (Vec(float.(tuple(0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * (u * (2 * v ^ (2 - 1))) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.5 + 0.5 * (2 * v ^ (2 - 1)) + 0.5 * u ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1)) + -0.5 * v ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), -0.5 * (u * (2 * v ^ (2 - 1))) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.5 + 0.5 * (2 * v ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1)) + 0.5 * v ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2), 0.5 * (u * (2 * v ^ (2 - 1))) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)) + (2 * u ^ (2 - 1)) * v ^ 2, -1.0 * (2 * v ^ (2 - 1)) + u ^ 2 * (2 * v ^ (2 - 1)))))) end # ────────────────────────────────────────────────────────────────────────────── -# Tet4: 4-node linear tetrahedral element +# Lagrange{Tetrahedron, 1}: 4-node linear tetrahedral element +# (Old name: Tet4) # ────────────────────────────────────────────────────────────────────────────── - struct Tet4Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Tet4Basis}) - return (3, 4) - end - function Base.size(::Type{Tet4Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 4 + function get_reference_element_coordinates(::Type{Lagrange{Tetrahedron, 1}}) + return (Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0))) end - Base.@pure function Base.length(::Type{Tet4Basis}) - return 4 - end - function get_reference_element_coordinates(::Type{Tet4Basis}) - return Vec{3, Float64}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Tetrahedron, 1}) + return (Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0))) end - @inline function eval_basis!(::Type{Tet4Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Tetrahedron, 1}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (1 + -1.0u + -1.0v + -1.0w, +u, +v, +w) end - @inline function eval_dbasis!(::Type{Tet4Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Tetrahedron, 1}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (1 + -1.0u + -1.0v + -1.0w, +u, +v, +w) + end + @inline function eval_dbasis!(::Type{Lagrange{Tetrahedron, 1}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-1.0, -1.0, -1.0))), Vec(float.(tuple(1, 0, 0))), Vec(float.(tuple(0, 1, 0))), Vec(float.(tuple(0, 0, 1)))) + end + @inline function eval_dbasis!(::Lagrange{Tetrahedron, 1}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-1.0, -1.0, -1.0))), Vec(float.(tuple(1, 0, 0))), Vec(float.(tuple(0, 1, 0))), Vec(float.(tuple(0, 0, 1)))) end # ────────────────────────────────────────────────────────────────────────────── -# Tet10: 10-node quadratic tetrahedral element +# Lagrange{Tetrahedron, 2}: 10-node quadratic tetrahedral element +# (Old name: Tet10) # ────────────────────────────────────────────────────────────────────────────── - struct Tet10Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Tet10Basis}) - return (3, 10) - end - function Base.size(::Type{Tet10Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 10 + function get_reference_element_coordinates(::Type{Lagrange{Tetrahedron, 2}}) + return (Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.0, 0.0)), Vec{3, Float64}(tuple(0.5, 0.5, 0.0)), Vec{3, Float64}(tuple(0.0, 0.5, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.5)), Vec{3, Float64}(tuple(0.5, 0.0, 0.5)), Vec{3, Float64}(tuple(0.0, 0.5, 0.5))) end - Base.@pure function Base.length(::Type{Tet10Basis}) - return 10 - end - function get_reference_element_coordinates(::Type{Tet10Basis}) - return Vec{3, Float64}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]] + function get_reference_element_coordinates(::Lagrange{Tetrahedron, 2}) + return (Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.0, 0.0)), Vec{3, Float64}(tuple(0.5, 0.5, 0.0)), Vec{3, Float64}(tuple(0.0, 0.5, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.5)), Vec{3, Float64}(tuple(0.5, 0.0, 0.5)), Vec{3, Float64}(tuple(0.0, 0.5, 0.5))) end - @inline function eval_basis!(::Type{Tet10Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Tetrahedron, 2}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (1 + -3.0u + -3.0v + -3.0w + 2.0 * u ^ 2 + 2.0 * v ^ 2 + 2.0 * w ^ 2 + 4.0 * (u * v) + 4.0 * (u * w) + 4.0 * (v * w), -1.0u + 2.0 * u ^ 2, -1.0v + 2.0 * v ^ 2, -1.0w + 2.0 * w ^ 2, 4.0u + -4.0 * u ^ 2 + -4.0 * (u * v) + -4.0 * (u * w), +(4.0 * (u * v)), 4.0v + -4.0 * v ^ 2 + -4.0 * (u * v) + -4.0 * (v * w), 4.0w + -4.0 * w ^ 2 + -4.0 * (u * w) + -4.0 * (v * w), +(4.0 * (u * w)), +(4.0 * (v * w))) end - @inline function eval_dbasis!(::Type{Tet10Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Tetrahedron, 2}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (1 + -3.0u + -3.0v + -3.0w + 2.0 * u ^ 2 + 2.0 * v ^ 2 + 2.0 * w ^ 2 + 4.0 * (u * v) + 4.0 * (u * w) + 4.0 * (v * w), -1.0u + 2.0 * u ^ 2, -1.0v + 2.0 * v ^ 2, -1.0w + 2.0 * w ^ 2, 4.0u + -4.0 * u ^ 2 + -4.0 * (u * v) + -4.0 * (u * w), +(4.0 * (u * v)), 4.0v + -4.0 * v ^ 2 + -4.0 * (u * v) + -4.0 * (v * w), 4.0w + -4.0 * w ^ 2 + -4.0 * (u * w) + -4.0 * (v * w), +(4.0 * (u * w)), +(4.0 * (v * w))) + end + @inline function eval_dbasis!(::Type{Lagrange{Tetrahedron, 2}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-3.0 + 2.0 * (2 * u ^ (2 - 1)) + 4.0v + 4.0w, -3.0 + 2.0 * (2 * v ^ (2 - 1)) + 4.0u + 4.0w, -3.0 + 2.0 * (2 * w ^ (2 - 1)) + 4.0u + 4.0v))), Vec(float.(tuple(-1.0 + 2.0 * (2 * u ^ (2 - 1)), 0, 0))), Vec(float.(tuple(0, -1.0 + 2.0 * (2 * v ^ (2 - 1)), 0))), Vec(float.(tuple(0, 0, -1.0 + 2.0 * (2 * w ^ (2 - 1))))), Vec(float.(tuple(4.0 + -4.0 * (2 * u ^ (2 - 1)) + -4.0v + -4.0w, -4.0u, -4.0u))), Vec(float.(tuple(4.0v, 4.0u, 0))), Vec(float.(tuple(-4.0v, 4.0 + -4.0 * (2 * v ^ (2 - 1)) + -4.0u + -4.0w, -4.0v))), Vec(float.(tuple(-4.0w, -4.0w, 4.0 + -4.0 * (2 * w ^ (2 - 1)) + -4.0u + -4.0v))), Vec(float.(tuple(4.0w, 0, 4.0u))), Vec(float.(tuple(0, 4.0w, 4.0v)))) + end + @inline function eval_dbasis!(::Lagrange{Tetrahedron, 2}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-3.0 + 2.0 * (2 * u ^ (2 - 1)) + 4.0v + 4.0w, -3.0 + 2.0 * (2 * v ^ (2 - 1)) + 4.0u + 4.0w, -3.0 + 2.0 * (2 * w ^ (2 - 1)) + 4.0u + 4.0v))), Vec(float.(tuple(-1.0 + 2.0 * (2 * u ^ (2 - 1)), 0, 0))), Vec(float.(tuple(0, -1.0 + 2.0 * (2 * v ^ (2 - 1)), 0))), Vec(float.(tuple(0, 0, -1.0 + 2.0 * (2 * w ^ (2 - 1))))), Vec(float.(tuple(4.0 + -4.0 * (2 * u ^ (2 - 1)) + -4.0v + -4.0w, -4.0u, -4.0u))), Vec(float.(tuple(4.0v, 4.0u, 0))), Vec(float.(tuple(-4.0v, 4.0 + -4.0 * (2 * v ^ (2 - 1)) + -4.0u + -4.0w, -4.0v))), Vec(float.(tuple(-4.0w, -4.0w, 4.0 + -4.0 * (2 * w ^ (2 - 1)) + -4.0u + -4.0v))), Vec(float.(tuple(4.0w, 0, 4.0u))), Vec(float.(tuple(0, 4.0w, 4.0v)))) end # ────────────────────────────────────────────────────────────────────────────── -# Hex8: 8-node trilinear hexahedral element +# Lagrange{Hexahedron, 1}: 8-node trilinear hexahedral element +# (Old name: Hex8) # ────────────────────────────────────────────────────────────────────────────── - struct Hex8Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Hex8Basis}) - return (3, 8) - end - function Base.size(::Type{Hex8Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 8 + function get_reference_element_coordinates(::Type{Lagrange{Hexahedron, 1}}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0))) end - Base.@pure function Base.length(::Type{Hex8Basis}) - return 8 - end - function get_reference_element_coordinates(::Type{Hex8Basis}) - return Vec{3, Float64}[[-1.0, -1.0, -1.0], [1.0, -1.0, -1.0], [1.0, 1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Hexahedron, 1}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0))) end - @inline function eval_basis!(::Type{Hex8Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Hexahedron, 1}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (0.125 + -0.125u + -0.125v + -0.125w + 0.125 * (u * v) + 0.125 * (u * w) + 0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + 0.125u + -0.125v + -0.125w + -0.125 * (u * v) + -0.125 * (u * w) + 0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + 0.125u + 0.125v + -0.125w + 0.125 * (u * v) + -0.125 * (u * w) + -0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + -0.125u + 0.125v + -0.125w + -0.125 * (u * v) + 0.125 * (u * w) + -0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + -0.125u + -0.125v + 0.125w + 0.125 * (u * v) + -0.125 * (u * w) + -0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + 0.125u + -0.125v + 0.125w + -0.125 * (u * v) + 0.125 * (u * w) + -0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + 0.125u + 0.125v + 0.125w + 0.125 * (u * v) + 0.125 * (u * w) + 0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + -0.125u + 0.125v + 0.125w + -0.125 * (u * v) + -0.125 * (u * w) + 0.125 * (v * w) + -0.125 * (u * v * w)) end - @inline function eval_dbasis!(::Type{Hex8Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Hexahedron, 1}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (0.125 + -0.125u + -0.125v + -0.125w + 0.125 * (u * v) + 0.125 * (u * w) + 0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + 0.125u + -0.125v + -0.125w + -0.125 * (u * v) + -0.125 * (u * w) + 0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + 0.125u + 0.125v + -0.125w + 0.125 * (u * v) + -0.125 * (u * w) + -0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + -0.125u + 0.125v + -0.125w + -0.125 * (u * v) + 0.125 * (u * w) + -0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + -0.125u + -0.125v + 0.125w + 0.125 * (u * v) + -0.125 * (u * w) + -0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + 0.125u + -0.125v + 0.125w + -0.125 * (u * v) + 0.125 * (u * w) + -0.125 * (v * w) + -0.125 * (u * v * w), 0.125 + 0.125u + 0.125v + 0.125w + 0.125 * (u * v) + 0.125 * (u * w) + 0.125 * (v * w) + 0.125 * (u * v * w), 0.125 + -0.125u + 0.125v + 0.125w + -0.125 * (u * v) + -0.125 * (u * w) + 0.125 * (v * w) + -0.125 * (u * v * w)) + end + @inline function eval_dbasis!(::Type{Lagrange{Hexahedron, 1}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-0.125 + 0.125v + 0.125w + -0.125 * (v * w), -0.125 + 0.125u + 0.125w + -0.125 * (u * w), -0.125 + 0.125u + 0.125v + -0.125 * (u * v)))), Vec(float.(tuple(0.125 + -0.125v + -0.125w + 0.125 * (v * w), -0.125 + -0.125u + 0.125w + 0.125 * (u * w), -0.125 + -0.125u + 0.125v + 0.125 * (u * v)))), Vec(float.(tuple(0.125 + 0.125v + -0.125w + -0.125 * (v * w), 0.125 + 0.125u + -0.125w + -0.125 * (u * w), -0.125 + -0.125u + -0.125v + -0.125 * (u * v)))), Vec(float.(tuple(-0.125 + -0.125v + 0.125w + 0.125 * (v * w), 0.125 + -0.125u + -0.125w + 0.125 * (u * w), -0.125 + 0.125u + -0.125v + 0.125 * (u * v)))), Vec(float.(tuple(-0.125 + 0.125v + -0.125w + 0.125 * (v * w), -0.125 + 0.125u + -0.125w + 0.125 * (u * w), 0.125 + -0.125u + -0.125v + 0.125 * (u * v)))), Vec(float.(tuple(0.125 + -0.125v + 0.125w + -0.125 * (v * w), -0.125 + -0.125u + -0.125w + -0.125 * (u * w), 0.125 + 0.125u + -0.125v + -0.125 * (u * v)))), Vec(float.(tuple(0.125 + 0.125v + 0.125w + 0.125 * (v * w), 0.125 + 0.125u + 0.125w + 0.125 * (u * w), 0.125 + 0.125u + 0.125v + 0.125 * (u * v)))), Vec(float.(tuple(-0.125 + -0.125v + -0.125w + -0.125 * (v * w), 0.125 + -0.125u + 0.125w + -0.125 * (u * w), 0.125 + -0.125u + 0.125v + -0.125 * (u * v))))) + end + @inline function eval_dbasis!(::Lagrange{Hexahedron, 1}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-0.125 + 0.125v + 0.125w + -0.125 * (v * w), -0.125 + 0.125u + 0.125w + -0.125 * (u * w), -0.125 + 0.125u + 0.125v + -0.125 * (u * v)))), Vec(float.(tuple(0.125 + -0.125v + -0.125w + 0.125 * (v * w), -0.125 + -0.125u + 0.125w + 0.125 * (u * w), -0.125 + -0.125u + 0.125v + 0.125 * (u * v)))), Vec(float.(tuple(0.125 + 0.125v + -0.125w + -0.125 * (v * w), 0.125 + 0.125u + -0.125w + -0.125 * (u * w), -0.125 + -0.125u + -0.125v + -0.125 * (u * v)))), Vec(float.(tuple(-0.125 + -0.125v + 0.125w + 0.125 * (v * w), 0.125 + -0.125u + -0.125w + 0.125 * (u * w), -0.125 + 0.125u + -0.125v + 0.125 * (u * v)))), Vec(float.(tuple(-0.125 + 0.125v + -0.125w + 0.125 * (v * w), -0.125 + 0.125u + -0.125w + 0.125 * (u * w), 0.125 + -0.125u + -0.125v + 0.125 * (u * v)))), Vec(float.(tuple(0.125 + -0.125v + 0.125w + -0.125 * (v * w), -0.125 + -0.125u + -0.125w + -0.125 * (u * w), 0.125 + 0.125u + -0.125v + -0.125 * (u * v)))), Vec(float.(tuple(0.125 + 0.125v + 0.125w + 0.125 * (v * w), 0.125 + 0.125u + 0.125w + 0.125 * (u * w), 0.125 + 0.125u + 0.125v + 0.125 * (u * v)))), Vec(float.(tuple(-0.125 + -0.125v + -0.125w + -0.125 * (v * w), 0.125 + -0.125u + 0.125w + -0.125 * (u * w), 0.125 + -0.125u + 0.125v + -0.125 * (u * v))))) end # ────────────────────────────────────────────────────────────────────────────── -# Hex20: 20-node serendipity hexahedral element +# Lagrange{Hexahedron, 2}: 20-node serendipity hexahedral element +# (Old name: Hex20) # ────────────────────────────────────────────────────────────────────────────── - struct Hex20Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Hex20Basis}) - return (3, 20) - end - function Base.size(::Type{Hex20Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 20 + function get_reference_element_coordinates(::Type{Lagrange{Hexahedron, 2}}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 1.0))) end - Base.@pure function Base.length(::Type{Hex20Basis}) - return 20 - end - function get_reference_element_coordinates(::Type{Hex20Basis}) - return Vec{3, Float64}[[-1.0, -1.0, -1.0], [1.0, -1.0, -1.0], [1.0, 1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0], [0.0, -1.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [-1.0, 0.0, -1.0], [-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, -1.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Hexahedron, 2}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 1.0))) end - @inline function eval_basis!(::Type{Hex20Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Hexahedron, 2}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (-0.25 + 0.125u + 0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + -0.125u + 0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + -0.125u + -0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + 0.125u + -0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + 0.125u + 0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + -0.125u + 0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + -0.125u + -0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + 0.125u + -0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), 0.25 + -0.25v + -0.25w + -0.25 * u ^ 2 + 0.25 * (v * w) + 0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * v * w), 0.25 + 0.25u + -0.25w + -0.25 * v ^ 2 + -0.25 * (u * w) + -0.25 * (v ^ 2 * u) + 0.25 * (v ^ 2 * w) + 0.25 * (u * v ^ 2 * w), 0.25 + 0.25v + -0.25w + -0.25 * u ^ 2 + -0.25 * (v * w) + -0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * v * w), 0.25 + -0.25u + -0.25w + -0.25 * v ^ 2 + 0.25 * (u * w) + 0.25 * (v ^ 2 * u) + 0.25 * (v ^ 2 * w) + -0.25 * (u * v ^ 2 * w), 0.25 + -0.25u + -0.25v + -0.25 * w ^ 2 + 0.25 * (u * v) + 0.25 * (w ^ 2 * u) + 0.25 * (w ^ 2 * v) + -0.25 * (u * v * w ^ 2), 0.25 + 0.25u + -0.25v + -0.25 * w ^ 2 + -0.25 * (u * v) + -0.25 * (w ^ 2 * u) + 0.25 * (w ^ 2 * v) + 0.25 * (u * v * w ^ 2), 0.25 + 0.25u + 0.25v + -0.25 * w ^ 2 + 0.25 * (u * v) + -0.25 * (w ^ 2 * u) + -0.25 * (w ^ 2 * v) + -0.25 * (u * v * w ^ 2), 0.25 + -0.25u + 0.25v + -0.25 * w ^ 2 + -0.25 * (u * v) + 0.25 * (w ^ 2 * u) + -0.25 * (w ^ 2 * v) + 0.25 * (u * v * w ^ 2), 0.25 + -0.25v + 0.25w + -0.25 * u ^ 2 + -0.25 * (v * w) + 0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * v * w), 0.25 + 0.25u + 0.25w + -0.25 * v ^ 2 + 0.25 * (u * w) + -0.25 * (v ^ 2 * u) + -0.25 * (v ^ 2 * w) + -0.25 * (u * v ^ 2 * w), 0.25 + 0.25v + 0.25w + -0.25 * u ^ 2 + 0.25 * (v * w) + -0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * v * w), 0.25 + -0.25u + 0.25w + -0.25 * v ^ 2 + -0.25 * (u * w) + 0.25 * (v ^ 2 * u) + -0.25 * (v ^ 2 * w) + 0.25 * (u * v ^ 2 * w)) end - @inline function eval_dbasis!(::Type{Hex20Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Hexahedron, 2}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (-0.25 + 0.125u + 0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + -0.125u + 0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + -0.125u + -0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + 0.125u + -0.125v + 0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + -0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + -0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + 0.125u + 0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + -0.125u + 0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + -0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), -0.25 + -0.125u + -0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + 0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + 0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2), -0.25 + 0.125u + -0.125v + -0.125w + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (u ^ 2 * v) + 0.125 * (u ^ 2 * w) + -0.125 * (v ^ 2 * u) + 0.125 * (v ^ 2 * w) + -0.125 * (w ^ 2 * u) + 0.125 * (w ^ 2 * v) + -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2), 0.25 + -0.25v + -0.25w + -0.25 * u ^ 2 + 0.25 * (v * w) + 0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * v * w), 0.25 + 0.25u + -0.25w + -0.25 * v ^ 2 + -0.25 * (u * w) + -0.25 * (v ^ 2 * u) + 0.25 * (v ^ 2 * w) + 0.25 * (u * v ^ 2 * w), 0.25 + 0.25v + -0.25w + -0.25 * u ^ 2 + -0.25 * (v * w) + -0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * v * w), 0.25 + -0.25u + -0.25w + -0.25 * v ^ 2 + 0.25 * (u * w) + 0.25 * (v ^ 2 * u) + 0.25 * (v ^ 2 * w) + -0.25 * (u * v ^ 2 * w), 0.25 + -0.25u + -0.25v + -0.25 * w ^ 2 + 0.25 * (u * v) + 0.25 * (w ^ 2 * u) + 0.25 * (w ^ 2 * v) + -0.25 * (u * v * w ^ 2), 0.25 + 0.25u + -0.25v + -0.25 * w ^ 2 + -0.25 * (u * v) + -0.25 * (w ^ 2 * u) + 0.25 * (w ^ 2 * v) + 0.25 * (u * v * w ^ 2), 0.25 + 0.25u + 0.25v + -0.25 * w ^ 2 + 0.25 * (u * v) + -0.25 * (w ^ 2 * u) + -0.25 * (w ^ 2 * v) + -0.25 * (u * v * w ^ 2), 0.25 + -0.25u + 0.25v + -0.25 * w ^ 2 + -0.25 * (u * v) + 0.25 * (w ^ 2 * u) + -0.25 * (w ^ 2 * v) + 0.25 * (u * v * w ^ 2), 0.25 + -0.25v + 0.25w + -0.25 * u ^ 2 + -0.25 * (v * w) + 0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * v * w), 0.25 + 0.25u + 0.25w + -0.25 * v ^ 2 + 0.25 * (u * w) + -0.25 * (v ^ 2 * u) + -0.25 * (v ^ 2 * w) + -0.25 * (u * v ^ 2 * w), 0.25 + 0.25v + 0.25w + -0.25 * u ^ 2 + 0.25 * (v * w) + -0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * v * w), 0.25 + -0.25u + 0.25w + -0.25 * v ^ 2 + -0.25 * (u * w) + 0.25 * (v ^ 2 * u) + -0.25 * (v ^ 2 * w) + 0.25 * (u * v ^ 2 * w)) + end + @inline function eval_dbasis!(::Type{Lagrange{Hexahedron, 2}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + -0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + 0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + 0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + -0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w), -0.25 + 0.25w + 0.25 * u ^ 2 + -0.25 * (u ^ 2 * w), -0.25 + 0.25v + 0.25 * u ^ 2 + -0.25 * (u ^ 2 * v)))), Vec(float.(tuple(0.25 + -0.25w + -0.25 * v ^ 2 + 0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w), -0.25 + -0.25u + 0.25 * v ^ 2 + 0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w), 0.25 + -0.25w + -0.25 * u ^ 2 + 0.25 * (u ^ 2 * w), -0.25 + -0.25v + 0.25 * u ^ 2 + 0.25 * (u ^ 2 * v)))), Vec(float.(tuple(-0.25 + 0.25w + 0.25 * v ^ 2 + -0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w), -0.25 + 0.25u + 0.25 * v ^ 2 + -0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 + 0.25v + 0.25 * w ^ 2 + -0.25 * (v * w ^ 2), -0.25 + 0.25u + 0.25 * w ^ 2 + -0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * ((2 * w ^ (2 - 1)) * v) + -0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 + -0.25v + -0.25 * w ^ 2 + 0.25 * (v * w ^ 2), -0.25 + -0.25u + 0.25 * w ^ 2 + 0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 + 0.25v + -0.25 * w ^ 2 + -0.25 * (v * w ^ 2), 0.25 + 0.25u + -0.25 * w ^ 2 + -0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + -0.25 * ((2 * w ^ (2 - 1)) * u) + -0.25 * ((2 * w ^ (2 - 1)) * v) + -0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 + -0.25v + 0.25 * w ^ 2 + 0.25 * (v * w ^ 2), 0.25 + -0.25u + -0.25 * w ^ 2 + 0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + 0.25 * ((2 * w ^ (2 - 1)) * u) + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w), -0.25 + -0.25w + 0.25 * u ^ 2 + 0.25 * (u ^ 2 * w), 0.25 + -0.25v + -0.25 * u ^ 2 + 0.25 * (u ^ 2 * v)))), Vec(float.(tuple(0.25 + 0.25w + -0.25 * v ^ 2 + -0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + -0.25 * ((2 * v ^ (2 - 1)) * u) + -0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w), 0.25 + 0.25u + -0.25 * v ^ 2 + -0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w), 0.25 + 0.25w + -0.25 * u ^ 2 + -0.25 * (u ^ 2 * w), 0.25 + 0.25v + -0.25 * u ^ 2 + -0.25 * (u ^ 2 * v)))), Vec(float.(tuple(-0.25 + -0.25w + 0.25 * v ^ 2 + 0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + 0.25 * ((2 * v ^ (2 - 1)) * u) + -0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w), 0.25 + -0.25u + -0.25 * v ^ 2 + 0.25 * (u * v ^ 2))))) + end + @inline function eval_dbasis!(::Lagrange{Hexahedron, 2}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + -0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + -0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + 0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + -0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), 0.125 + 0.125 * (2 * w ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + 0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + -0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + -0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), 0.125 + 0.125 * (2 * v ^ (2 - 1)) + -0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + -0.125 * w ^ 2 + -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + -0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + 0.125 * v ^ 2 + 0.125 * w ^ 2 + 0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + 0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 + 0.125 * (2 * u ^ (2 - 1)) + 0.125 * ((2 * u ^ (2 - 1)) * v) + 0.125 * ((2 * u ^ (2 - 1)) * w) + -0.125 * v ^ 2 + -0.125 * w ^ 2 + -0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2), -0.125 + 0.125 * (2 * v ^ (2 - 1)) + 0.125 * u ^ 2 + -0.125 * ((2 * v ^ (2 - 1)) * u) + 0.125 * ((2 * v ^ (2 - 1)) * w) + 0.125 * w ^ 2 + -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2), -0.125 + 0.125 * (2 * w ^ (2 - 1)) + 0.125 * u ^ 2 + 0.125 * v ^ 2 + -0.125 * ((2 * w ^ (2 - 1)) * u) + 0.125 * ((2 * w ^ (2 - 1)) * v) + -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w), -0.25 + 0.25w + 0.25 * u ^ 2 + -0.25 * (u ^ 2 * w), -0.25 + 0.25v + 0.25 * u ^ 2 + -0.25 * (u ^ 2 * v)))), Vec(float.(tuple(0.25 + -0.25w + -0.25 * v ^ 2 + 0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w), -0.25 + -0.25u + 0.25 * v ^ 2 + 0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w), 0.25 + -0.25w + -0.25 * u ^ 2 + 0.25 * (u ^ 2 * w), -0.25 + -0.25v + 0.25 * u ^ 2 + 0.25 * (u ^ 2 * v)))), Vec(float.(tuple(-0.25 + 0.25w + 0.25 * v ^ 2 + -0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w), -0.25 + 0.25u + 0.25 * v ^ 2 + -0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 + 0.25v + 0.25 * w ^ 2 + -0.25 * (v * w ^ 2), -0.25 + 0.25u + 0.25 * w ^ 2 + -0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * ((2 * w ^ (2 - 1)) * v) + -0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 + -0.25v + -0.25 * w ^ 2 + 0.25 * (v * w ^ 2), -0.25 + -0.25u + 0.25 * w ^ 2 + 0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 + 0.25v + -0.25 * w ^ 2 + -0.25 * (v * w ^ 2), 0.25 + 0.25u + -0.25 * w ^ 2 + -0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + -0.25 * ((2 * w ^ (2 - 1)) * u) + -0.25 * ((2 * w ^ (2 - 1)) * v) + -0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 + -0.25v + 0.25 * w ^ 2 + 0.25 * (v * w ^ 2), 0.25 + -0.25u + -0.25 * w ^ 2 + 0.25 * (u * w ^ 2), -0.25 * (2 * w ^ (2 - 1)) + 0.25 * ((2 * w ^ (2 - 1)) * u) + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (u * v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w), -0.25 + -0.25w + 0.25 * u ^ 2 + 0.25 * (u ^ 2 * w), 0.25 + -0.25v + -0.25 * u ^ 2 + 0.25 * (u ^ 2 * v)))), Vec(float.(tuple(0.25 + 0.25w + -0.25 * v ^ 2 + -0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + -0.25 * ((2 * v ^ (2 - 1)) * u) + -0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w), 0.25 + 0.25u + -0.25 * v ^ 2 + -0.25 * (u * v ^ 2)))), Vec(float.(tuple(-0.25 * (2 * u ^ (2 - 1)) + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w), 0.25 + 0.25w + -0.25 * u ^ 2 + -0.25 * (u ^ 2 * w), 0.25 + 0.25v + -0.25 * u ^ 2 + -0.25 * (u ^ 2 * v)))), Vec(float.(tuple(-0.25 + -0.25w + 0.25 * v ^ 2 + 0.25 * (v ^ 2 * w), -0.25 * (2 * v ^ (2 - 1)) + 0.25 * ((2 * v ^ (2 - 1)) * u) + -0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w), 0.25 + -0.25u + -0.25 * v ^ 2 + 0.25 * (u * v ^ 2))))) end # ────────────────────────────────────────────────────────────────────────────── -# Hex27: 27-node triquadratic hexahedral element +# Lagrange{Hexahedron, 2}: 27-node triquadratic hexahedral element +# (Old name: Hex27) # ────────────────────────────────────────────────────────────────────────────── - struct Hex27Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Hex27Basis}) - return (3, 27) - end - function Base.size(::Type{Hex27Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 27 + function get_reference_element_coordinates(::Type{Lagrange{Hexahedron, 2}}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.0))) end - Base.@pure function Base.length(::Type{Hex27Basis}) - return 27 - end - function get_reference_element_coordinates(::Type{Hex27Basis}) - return Vec{3, Float64}[[-1.0, -1.0, -1.0], [1.0, -1.0, -1.0], [1.0, 1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0], [0.0, -1.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [-1.0, 0.0, -1.0], [-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, -1.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0], [0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0]] + function get_reference_element_coordinates(::Lagrange{Hexahedron, 2}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, -1.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.0))) end - @inline function eval_basis!(::Type{Hex27Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Hexahedron, 2}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (-0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (v * w) + -0.25 * (v ^ 2 * w) + -0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * w) + -0.25 * (u ^ 2 * w) + 0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * v ^ 2 * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (v * w) + -0.25 * (v ^ 2 * w) + 0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + 0.25 * (u ^ 2 * v * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * w) + -0.25 * (u ^ 2 * w) + -0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * v ^ 2 * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * v) + -0.25 * (u ^ 2 * v) + -0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v * w ^ 2) + 0.25 * (u ^ 2 * v * w ^ 2) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * v) + -0.25 * (u ^ 2 * v) + 0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v * w ^ 2) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * v) + 0.25 * (u ^ 2 * v) + 0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v * w ^ 2) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * v) + 0.25 * (u ^ 2 * v) + -0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v * w ^ 2) + -0.25 * (u ^ 2 * v * w ^ 2) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (v * w) + 0.25 * (v ^ 2 * w) + -0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + 0.25 * (u ^ 2 * v * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * w) + 0.25 * (u ^ 2 * w) + 0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * v ^ 2 * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (v * w) + 0.25 * (v ^ 2 * w) + 0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * w) + 0.25 * (u ^ 2 * w) + -0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * v ^ 2 * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5w + 0.5 * w ^ 2 + 0.5 * (u ^ 2 * w) + 0.5 * (v ^ 2 * w) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + -0.5 * (u ^ 2 * v ^ 2 * w) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5w + 0.5 * w ^ 2 + -0.5 * (u ^ 2 * w) + -0.5 * (v ^ 2 * w) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5v + 0.5 * v ^ 2 + 0.5 * (u ^ 2 * v) + 0.5 * (w ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + -0.5 * (u ^ 2 * v * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5u + 0.5 * u ^ 2 + -0.5 * (v ^ 2 * u) + -0.5 * (w ^ 2 * u) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u * v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5v + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * v) + -0.5 * (w ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5u + 0.5 * u ^ 2 + 0.5 * (v ^ 2 * u) + 0.5 * (w ^ 2 * u) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (u * v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 1 + -1.0 * u ^ 2 + -1.0 * v ^ 2 + -1.0 * w ^ 2 + u ^ 2 * v ^ 2 + u ^ 2 * w ^ 2 + v ^ 2 * w ^ 2 + -1.0 * (u ^ 2 * v ^ 2 * w ^ 2)) end - @inline function eval_dbasis!(::Type{Hex27Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Hexahedron, 2}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (-0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + -0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + -0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + -0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + 0.125 * (u * v ^ 2 * w) + 0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + 0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), -0.125 * (u * v * w) + 0.125 * (u ^ 2 * v * w) + -0.125 * (u * v ^ 2 * w) + -0.125 * (u * v * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w) + 0.125 * (u ^ 2 * v * w ^ 2) + -0.125 * (u * v ^ 2 * w ^ 2) + 0.125 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (v * w) + -0.25 * (v ^ 2 * w) + -0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * w) + -0.25 * (u ^ 2 * w) + 0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * v ^ 2 * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (v * w) + -0.25 * (v ^ 2 * w) + 0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + 0.25 * (u ^ 2 * v * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * w) + -0.25 * (u ^ 2 * w) + -0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * v ^ 2 * w) + 0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * v) + -0.25 * (u ^ 2 * v) + -0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v * w ^ 2) + 0.25 * (u ^ 2 * v * w ^ 2) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * v) + -0.25 * (u ^ 2 * v) + 0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v * w ^ 2) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * v) + 0.25 * (u ^ 2 * v) + 0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v * w ^ 2) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * v) + 0.25 * (u ^ 2 * v) + -0.25 * (v ^ 2 * u) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v * w ^ 2) + -0.25 * (u ^ 2 * v * w ^ 2) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (v * w) + 0.25 * (v ^ 2 * w) + -0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + 0.25 * (u ^ 2 * v * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (u * w) + 0.25 * (u ^ 2 * w) + 0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * v ^ 2 * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), 0.25 * (v * w) + 0.25 * (v ^ 2 * w) + 0.25 * (w ^ 2 * v) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + -0.25 * (u ^ 2 * v * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.25 * (u * w) + 0.25 * (u ^ 2 * w) + -0.25 * (w ^ 2 * u) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * v ^ 2 * w) + -0.25 * (u ^ 2 * v ^ 2 * w) + 0.25 * (u * v ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5w + 0.5 * w ^ 2 + 0.5 * (u ^ 2 * w) + 0.5 * (v ^ 2 * w) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + -0.5 * (u ^ 2 * v ^ 2 * w) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5w + 0.5 * w ^ 2 + -0.5 * (u ^ 2 * w) + -0.5 * (v ^ 2 * w) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5v + 0.5 * v ^ 2 + 0.5 * (u ^ 2 * v) + 0.5 * (w ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + -0.5 * (u ^ 2 * v * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5u + 0.5 * u ^ 2 + -0.5 * (v ^ 2 * u) + -0.5 * (w ^ 2 * u) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u * v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 0.5v + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * v) + -0.5 * (w ^ 2 * v) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), -0.5u + 0.5 * u ^ 2 + 0.5 * (v ^ 2 * u) + 0.5 * (w ^ 2 * u) + -0.5 * (u ^ 2 * v ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + -0.5 * (u * v ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * w ^ 2), 1 + -1.0 * u ^ 2 + -1.0 * v ^ 2 + -1.0 * w ^ 2 + u ^ 2 * v ^ 2 + u ^ 2 * w ^ 2 + v ^ 2 * w ^ 2 + -1.0 * (u ^ 2 * v ^ 2 * w ^ 2)) + end + @inline function eval_dbasis!(::Type{Lagrange{Hexahedron, 2}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * ((2 * u ^ (2 - 1)) * v * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25w + -0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25v + -0.25 * v ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25w + -0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.25 * (v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25 * (u * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 * ((2 * u ^ (2 - 1)) * v * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25w + -0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25v + -0.25 * v ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25w + -0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.25 * (v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25 * (u * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.25 * (v * w ^ 2) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.25 * (u * w ^ 2) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25 * (u * v * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.25 * (v * w ^ 2) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.25 * (u * w ^ 2) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25 * (u * v * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.25 * (v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.25 * (u * w ^ 2) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25 * (u * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.25 * (v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.25 * (u * w ^ 2) + -0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25 * (u * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 * ((2 * u ^ (2 - 1)) * v * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25w + 0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25v + 0.25 * v ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25w + 0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.25 * (v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25 * (u * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2) + -0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * ((2 * u ^ (2 - 1)) * v * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25w + 0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25v + 0.25 * v ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25w + 0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.25 * (v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25 * (u * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2) + -0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * w) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 * ((2 * v ^ (2 - 1)) * w) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 + 0.5 * (2 * w ^ (2 - 1)) + 0.5 * u ^ 2 + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u ^ 2 * v ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * w) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 * ((2 * v ^ (2 - 1)) * w) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 + 0.5 * (2 * w ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * v ^ 2 + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 + 0.5 * (2 * v ^ (2 - 1)) + 0.5 * u ^ 2 + 0.5 * w ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 * ((2 * w ^ (2 - 1)) * v) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1)) + -0.5 * v ^ 2 + -0.5 * w ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.5 * (v ^ 2 * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 * ((2 * v ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.5 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 * ((2 * w ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 + 0.5 * (2 * v ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * w ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 * ((2 * w ^ (2 - 1)) * v) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1)) + 0.5 * v ^ 2 + 0.5 * w ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 * ((2 * v ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 * ((2 * w ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)) + (2 * u ^ (2 - 1)) * v ^ 2 + (2 * u ^ (2 - 1)) * w ^ 2 + -1.0 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -1.0 * (2 * v ^ (2 - 1)) + u ^ 2 * (2 * v ^ (2 - 1)) + (2 * v ^ (2 - 1)) * w ^ 2 + -1.0 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -1.0 * (2 * w ^ (2 - 1)) + u ^ 2 * (2 * w ^ (2 - 1)) + v ^ 2 * (2 * w ^ (2 - 1)) + -1.0 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1))))))) + end + @inline function eval_dbasis!(::Lagrange{Hexahedron, 2}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + -0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + -0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + -0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + -0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + -0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + -0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + -0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + -0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + 0.125 * (v ^ 2 * w) + 0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.125 * (u * w) + 0.125 * (u ^ 2 * w) + 0.125 * (u * (2 * v ^ (2 - 1)) * w) + 0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + 0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.125 * (u * v) + 0.125 * (u ^ 2 * v) + 0.125 * (u * v ^ 2) + 0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.125 * (v * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w) + -0.125 * (v ^ 2 * w) + -0.125 * (v * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.125 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.125 * (v ^ 2 * w ^ 2) + 0.125 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.125 * (u * w) + 0.125 * (u ^ 2 * w) + -0.125 * (u * (2 * v ^ (2 - 1)) * w) + -0.125 * (u * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.125 * (u ^ 2 * w ^ 2) + -0.125 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.125 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.125 * (u * v) + 0.125 * (u ^ 2 * v) + -0.125 * (u * v ^ 2) + -0.125 * (u * v * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2) + 0.125 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.125 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.125 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * ((2 * u ^ (2 - 1)) * v * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25w + -0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25v + -0.25 * v ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25w + -0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.25 * (v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25 * (u * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 * ((2 * u ^ (2 - 1)) * v * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25w + -0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.25 * (u ^ 2 * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25v + -0.25 * v ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v) + 0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25w + -0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.25 * (v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25 * (u * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2) + 0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.25 * (v * w ^ 2) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25u + -0.25 * u ^ 2 + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.25 * (u * w ^ 2) + 0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25 * (u * v * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + -0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.25 * (v * w ^ 2) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25u + -0.25 * u ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.25 * (u * w ^ 2) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25 * (u * v * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + 0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.25 * (v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.25 * (u * w ^ 2) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25 * (u * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25v + 0.25 * ((2 * u ^ (2 - 1)) * v) + -0.25 * v ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.25 * (v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * ((2 * v ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.25 * (u * w ^ 2) + -0.25 * (u ^ 2 * w ^ 2) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25 * (u * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25 * ((2 * u ^ (2 - 1)) * v * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25w + 0.25 * ((2 * v ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25v + 0.25 * v ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.25w + 0.25 * ((2 * u ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.25 * (v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.25 * (u * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25u + 0.25 * u ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u * v ^ 2) + -0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25 * ((2 * u ^ (2 - 1)) * v * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25w + 0.25 * ((2 * v ^ (2 - 1)) * w) + 0.25 * w ^ 2 + 0.25 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.25v + 0.25 * v ^ 2 + 0.25 * ((2 * w ^ (2 - 1)) * v) + 0.25 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v) + -0.25 * (u ^ 2 * v ^ 2) + -0.25 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.25w + 0.25 * ((2 * u ^ (2 - 1)) * w) + -0.25 * w ^ 2 + 0.25 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.25 * (v ^ 2 * w) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.25 * (v ^ 2 * w ^ 2) + -0.25 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.25 * (u * (2 * v ^ (2 - 1)) * w) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.25 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + -0.25 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.25u + 0.25 * u ^ 2 + -0.25 * ((2 * w ^ (2 - 1)) * u) + 0.25 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.25 * (u * v ^ 2) + -0.25 * (u ^ 2 * v ^ 2) + 0.25 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + -0.25 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * w) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 * ((2 * v ^ (2 - 1)) * w) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 + 0.5 * (2 * w ^ (2 - 1)) + 0.5 * u ^ 2 + 0.5 * v ^ 2 + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u ^ 2 * v ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * w) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 * ((2 * v ^ (2 - 1)) * w) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 + 0.5 * (2 * w ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * v ^ 2 + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 + 0.5 * (2 * v ^ (2 - 1)) + 0.5 * u ^ 2 + 0.5 * w ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + -0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 * ((2 * w ^ (2 - 1)) * v) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0.5 + 0.5 * (2 * u ^ (2 - 1)) + -0.5 * v ^ 2 + -0.5 * w ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + 0.5 * (v ^ 2 * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -0.5 * ((2 * v ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + 0.5 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 * ((2 * w ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 * ((2 * u ^ (2 - 1)) * v) + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 + 0.5 * (2 * v ^ (2 - 1)) + -0.5 * u ^ 2 + -0.5 * w ^ 2 + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * ((2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -0.5 * ((2 * w ^ (2 - 1)) * v) + -0.5 * (v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-0.5 + 0.5 * (2 * u ^ (2 - 1)) + 0.5 * v ^ 2 + 0.5 * w ^ 2 + -0.5 * ((2 * u ^ (2 - 1)) * v ^ 2) + -0.5 * ((2 * u ^ (2 - 1)) * w ^ 2) + -0.5 * (v ^ 2 * w ^ 2) + 0.5 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), 0.5 * ((2 * v ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * v ^ (2 - 1))) + -0.5 * (u * (2 * v ^ (2 - 1)) * w ^ 2) + 0.5 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), 0.5 * ((2 * w ^ (2 - 1)) * u) + -0.5 * (u ^ 2 * (2 * w ^ (2 - 1))) + -0.5 * (u * v ^ 2 * (2 * w ^ (2 - 1))) + 0.5 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 * (2 * u ^ (2 - 1)) + (2 * u ^ (2 - 1)) * v ^ 2 + (2 * u ^ (2 - 1)) * w ^ 2 + -1.0 * ((2 * u ^ (2 - 1)) * v ^ 2 * w ^ 2), -1.0 * (2 * v ^ (2 - 1)) + u ^ 2 * (2 * v ^ (2 - 1)) + (2 * v ^ (2 - 1)) * w ^ 2 + -1.0 * (u ^ 2 * (2 * v ^ (2 - 1)) * w ^ 2), -1.0 * (2 * w ^ (2 - 1)) + u ^ 2 * (2 * w ^ (2 - 1)) + v ^ 2 * (2 * w ^ (2 - 1)) + -1.0 * (u ^ 2 * v ^ 2 * (2 * w ^ (2 - 1))))))) end # ────────────────────────────────────────────────────────────────────────────── -# Pyr5: 5-node linear pyramid element +# Lagrange{Pyramid, 1}: 5-node linear pyramid element +# (Old name: Pyr5) # ────────────────────────────────────────────────────────────────────────────── - struct Pyr5Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Pyr5Basis}) - return (3, 5) - end - function Base.size(::Type{Pyr5Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 5 + function get_reference_element_coordinates(::Type{Lagrange{Pyramid, 1}}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0))) end - Base.@pure function Base.length(::Type{Pyr5Basis}) - return 5 - end - function get_reference_element_coordinates(::Type{Pyr5Basis}) - return Vec{3, Float64}[[-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Pyramid, 1}) + return (Vec{3, Float64}(tuple(-1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, -1.0, 0.0)), Vec{3, Float64}(tuple(1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(-1.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0))) end - @inline function eval_basis!(::Type{Pyr5Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Pyramid, 1}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (0.25 + -0.25u + -0.25v + -0.25w + 0.25 * (u * v), 0.25 + 0.25u + -0.25v + -0.25w + -0.25 * (u * v), 0.25 + 0.25u + 0.25v + -0.25w + 0.25 * (u * v), 0.25 + -0.25u + 0.25v + -0.25w + -0.25 * (u * v), +w) end - @inline function eval_dbasis!(::Type{Pyr5Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Pyramid, 1}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (0.25 + -0.25u + -0.25v + -0.25w + 0.25 * (u * v), 0.25 + 0.25u + -0.25v + -0.25w + -0.25 * (u * v), 0.25 + 0.25u + 0.25v + -0.25w + 0.25 * (u * v), 0.25 + -0.25u + 0.25v + -0.25w + -0.25 * (u * v), +w) + end + @inline function eval_dbasis!(::Type{Lagrange{Pyramid, 1}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-0.25 + 0.25v, -0.25 + 0.25u, -0.25))), Vec(float.(tuple(0.25 + -0.25v, -0.25 + -0.25u, -0.25))), Vec(float.(tuple(0.25 + 0.25v, 0.25 + 0.25u, -0.25))), Vec(float.(tuple(-0.25 + -0.25v, 0.25 + -0.25u, -0.25))), Vec(float.(tuple(0, 0, 1)))) + end + @inline function eval_dbasis!(::Lagrange{Pyramid, 1}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-0.25 + 0.25v, -0.25 + 0.25u, -0.25))), Vec(float.(tuple(0.25 + -0.25v, -0.25 + -0.25u, -0.25))), Vec(float.(tuple(0.25 + 0.25v, 0.25 + 0.25u, -0.25))), Vec(float.(tuple(-0.25 + -0.25v, 0.25 + -0.25u, -0.25))), Vec(float.(tuple(0, 0, 1)))) end # ────────────────────────────────────────────────────────────────────────────── -# Wedge6: 6-node linear wedge element (triangular prism) +# Lagrange{Wedge, 1}: 6-node linear wedge element (triangular prism) +# (Old name: Wedge6) # ────────────────────────────────────────────────────────────────────────────── - struct Wedge6Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Wedge6Basis}) - return (3, 6) - end - function Base.size(::Type{Wedge6Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 6 + function get_reference_element_coordinates(::Type{Lagrange{Wedge, 1}}) + return (Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0))) end - Base.@pure function Base.length(::Type{Wedge6Basis}) - return 6 - end - function get_reference_element_coordinates(::Type{Wedge6Basis}) - return Vec{3, Float64}[[0.0, 0.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0]] + function get_reference_element_coordinates(::Lagrange{Wedge, 1}) + return (Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0))) end - @inline function eval_basis!(::Type{Wedge6Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Wedge, 1}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (0.5 + -0.5u + -0.5v + -0.5w + 0.5 * (u * w) + 0.5 * (v * w), 0.5u + -0.5 * (u * w), 0.5v + -0.5 * (v * w), 0.5 + -0.5u + -0.5v + 0.5w + -0.5 * (u * w) + -0.5 * (v * w), 0.5u + 0.5 * (u * w), 0.5v + 0.5 * (v * w)) end - @inline function eval_dbasis!(::Type{Wedge6Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Wedge, 1}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (0.5 + -0.5u + -0.5v + -0.5w + 0.5 * (u * w) + 0.5 * (v * w), 0.5u + -0.5 * (u * w), 0.5v + -0.5 * (v * w), 0.5 + -0.5u + -0.5v + 0.5w + -0.5 * (u * w) + -0.5 * (v * w), 0.5u + 0.5 * (u * w), 0.5v + 0.5 * (v * w)) + end + @inline function eval_dbasis!(::Type{Lagrange{Wedge, 1}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-0.5 + 0.5w, -0.5 + 0.5w, -0.5 + 0.5u + 0.5v))), Vec(float.(tuple(0.5 + -0.5w, 0, -0.5u))), Vec(float.(tuple(0, 0.5 + -0.5w, -0.5v))), Vec(float.(tuple(-0.5 + -0.5w, -0.5 + -0.5w, 0.5 + -0.5u + -0.5v))), Vec(float.(tuple(0.5 + 0.5w, 0, 0.5u))), Vec(float.(tuple(0, 0.5 + 0.5w, 0.5v)))) + end + @inline function eval_dbasis!(::Lagrange{Wedge, 1}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-0.5 + 0.5w, -0.5 + 0.5w, -0.5 + 0.5u + 0.5v))), Vec(float.(tuple(0.5 + -0.5w, 0, -0.5u))), Vec(float.(tuple(0, 0.5 + -0.5w, -0.5v))), Vec(float.(tuple(-0.5 + -0.5w, -0.5 + -0.5w, 0.5 + -0.5u + -0.5v))), Vec(float.(tuple(0.5 + 0.5w, 0, 0.5u))), Vec(float.(tuple(0, 0.5 + 0.5w, 0.5v)))) end # ────────────────────────────────────────────────────────────────────────────── -# Wedge15: 15-node quadratic wedge element +# Lagrange{Wedge, 2}: 15-node quadratic wedge element +# (Old name: Wedge15) # ────────────────────────────────────────────────────────────────────────────── - struct Wedge15Basis <: AbstractBasis{3} - end - Base.@pure function Base.size(::Type{Wedge15Basis}) - return (3, 15) - end - function Base.size(::Type{Wedge15Basis}, j::Int) - j == 1 && return 3 - j == 2 && return 15 + function get_reference_element_coordinates(::Type{Lagrange{Wedge, 2}}) + return (Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.0, -1.0)), Vec{3, Float64}(tuple(0.5, 0.5, -1.0)), Vec{3, Float64}(tuple(0.0, 0.5, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.5, 0.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.5, 1.0)), Vec{3, Float64}(tuple(0.0, 0.5, 1.0))) end - Base.@pure function Base.length(::Type{Wedge15Basis}) - return 15 - end - function get_reference_element_coordinates(::Type{Wedge15Basis}) - return Vec{3, Float64}[[0.0, 0.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [0.5, 0.0, -1.0], [0.5, 0.5, -1.0], [0.0, 0.5, -1.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.5, 0.0, 1.0], [0.5, 0.5, 1.0], [0.0, 0.5, 1.0]] + function get_reference_element_coordinates(::Lagrange{Wedge, 2}) + return (Vec{3, Float64}(tuple(0.0, 0.0, -1.0)), Vec{3, Float64}(tuple(1.0, 0.0, -1.0)), Vec{3, Float64}(tuple(0.0, 1.0, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 1.0)), Vec{3, Float64}(tuple(1.0, 0.0, 1.0)), Vec{3, Float64}(tuple(0.0, 1.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.0, -1.0)), Vec{3, Float64}(tuple(0.5, 0.5, -1.0)), Vec{3, Float64}(tuple(0.0, 0.5, -1.0)), Vec{3, Float64}(tuple(0.0, 0.0, 0.0)), Vec{3, Float64}(tuple(1.0, 0.0, 0.0)), Vec{3, Float64}(tuple(0.0, 1.0, 0.0)), Vec{3, Float64}(tuple(0.5, 0.0, 1.0)), Vec{3, Float64}(tuple(0.5, 0.5, 1.0)), Vec{3, Float64}(tuple(0.0, 0.5, 1.0))) end - @inline function eval_basis!(::Type{Wedge15Basis}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{Wedge, 2}}, ::Type{T}, xi::Vec) where T (u, v, w) = xi @inbounds return (-1.0u + -1.0v + -0.5w + u ^ 2 + v ^ 2 + 0.5 * w ^ 2 + 2.0 * (u * v) + 1.5 * (u * w) + 1.5 * (v * w) + -1.0 * (u ^ 2 * w) + -1.0 * (v ^ 2 * w) + -2.0 * (u * v * w) + -0.5 * (u * w ^ 2) + -0.5 * (v * w ^ 2), -1.0u + u ^ 2 + 0.5 * (u * w) + -1.0 * (u ^ 2 * w) + 0.5 * (u * w ^ 2), -1.0v + v ^ 2 + 0.5 * (v * w) + -1.0 * (v ^ 2 * w) + 0.5 * (v * w ^ 2), -1.0u + -1.0v + 0.5w + u ^ 2 + v ^ 2 + 0.5 * w ^ 2 + 2.0 * (u * v) + -1.5 * (u * w) + -1.5 * (v * w) + u ^ 2 * w + v ^ 2 * w + 2.0 * (u * v * w) + -0.5 * (u * w ^ 2) + -0.5 * (v * w ^ 2), -1.0u + u ^ 2 + -0.5 * (u * w) + u ^ 2 * w + 0.5 * (u * w ^ 2), -1.0v + v ^ 2 + -0.5 * (v * w) + v ^ 2 * w + 0.5 * (v * w ^ 2), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v) + -2.0 * (u * w) + 2.0 * (u ^ 2 * w) + 2.0 * (u * v * w), 2.0 * (u * v) + -2.0 * (u * v * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v) + -2.0 * (v * w) + 2.0 * (v ^ 2 * w) + 2.0 * (u * v * w), 1 + -1.0u + -1.0v + -1.0 * w ^ 2 + u * w ^ 2 + v * w ^ 2, u + -1.0 * (u * w ^ 2), v + -1.0 * (v * w ^ 2), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v) + 2.0 * (u * w) + -2.0 * (u ^ 2 * w) + -2.0 * (u * v * w), 2.0 * (u * v) + 2.0 * (u * v * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v) + 2.0 * (v * w) + -2.0 * (v ^ 2 * w) + -2.0 * (u * v * w)) end - @inline function eval_dbasis!(::Type{Wedge15Basis}, xi::Vec) + @inline function eval_basis!(::Lagrange{Wedge, 2}, ::Type{T}, xi::Vec) where T + (u, v, w) = xi + @inbounds return (-1.0u + -1.0v + -0.5w + u ^ 2 + v ^ 2 + 0.5 * w ^ 2 + 2.0 * (u * v) + 1.5 * (u * w) + 1.5 * (v * w) + -1.0 * (u ^ 2 * w) + -1.0 * (v ^ 2 * w) + -2.0 * (u * v * w) + -0.5 * (u * w ^ 2) + -0.5 * (v * w ^ 2), -1.0u + u ^ 2 + 0.5 * (u * w) + -1.0 * (u ^ 2 * w) + 0.5 * (u * w ^ 2), -1.0v + v ^ 2 + 0.5 * (v * w) + -1.0 * (v ^ 2 * w) + 0.5 * (v * w ^ 2), -1.0u + -1.0v + 0.5w + u ^ 2 + v ^ 2 + 0.5 * w ^ 2 + 2.0 * (u * v) + -1.5 * (u * w) + -1.5 * (v * w) + u ^ 2 * w + v ^ 2 * w + 2.0 * (u * v * w) + -0.5 * (u * w ^ 2) + -0.5 * (v * w ^ 2), -1.0u + u ^ 2 + -0.5 * (u * w) + u ^ 2 * w + 0.5 * (u * w ^ 2), -1.0v + v ^ 2 + -0.5 * (v * w) + v ^ 2 * w + 0.5 * (v * w ^ 2), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v) + -2.0 * (u * w) + 2.0 * (u ^ 2 * w) + 2.0 * (u * v * w), 2.0 * (u * v) + -2.0 * (u * v * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v) + -2.0 * (v * w) + 2.0 * (v ^ 2 * w) + 2.0 * (u * v * w), 1 + -1.0u + -1.0v + -1.0 * w ^ 2 + u * w ^ 2 + v * w ^ 2, u + -1.0 * (u * w ^ 2), v + -1.0 * (v * w ^ 2), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v) + 2.0 * (u * w) + -2.0 * (u ^ 2 * w) + -2.0 * (u * v * w), 2.0 * (u * v) + 2.0 * (u * v * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v) + 2.0 * (v * w) + -2.0 * (v ^ 2 * w) + -2.0 * (u * v * w)) + end + @inline function eval_dbasis!(::Type{Lagrange{Wedge, 2}}, xi::Vec) + (u, v, w) = xi + @inbounds return (Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 2.0v + 1.5w + -1.0 * ((2 * u ^ (2 - 1)) * w) + -2.0 * (v * w) + -0.5 * w ^ 2, -1.0 + 2 * v ^ (2 - 1) + 2.0u + 1.5w + -1.0 * ((2 * v ^ (2 - 1)) * w) + -2.0 * (u * w) + -0.5 * w ^ 2, -0.5 + 0.5 * (2 * w ^ (2 - 1)) + 1.5u + 1.5v + -1.0 * u ^ 2 + -1.0 * v ^ 2 + -2.0 * (u * v) + -0.5 * (u * (2 * w ^ (2 - 1))) + -0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 0.5w + -1.0 * ((2 * u ^ (2 - 1)) * w) + 0.5 * w ^ 2, 0, 0.5u + -1.0 * u ^ 2 + 0.5 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, -1.0 + 2 * v ^ (2 - 1) + 0.5w + -1.0 * ((2 * v ^ (2 - 1)) * w) + 0.5 * w ^ 2, 0.5v + -1.0 * v ^ 2 + 0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 2.0v + -1.5w + (2 * u ^ (2 - 1)) * w + 2.0 * (v * w) + -0.5 * w ^ 2, -1.0 + 2 * v ^ (2 - 1) + 2.0u + -1.5w + (2 * v ^ (2 - 1)) * w + 2.0 * (u * w) + -0.5 * w ^ 2, 0.5 + 0.5 * (2 * w ^ (2 - 1)) + -1.5u + -1.5v + u ^ 2 + v ^ 2 + 2.0 * (u * v) + -0.5 * (u * (2 * w ^ (2 - 1))) + -0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + -0.5w + (2 * u ^ (2 - 1)) * w + 0.5 * w ^ 2, 0, -0.5u + u ^ 2 + 0.5 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, -1.0 + 2 * v ^ (2 - 1) + -0.5w + (2 * v ^ (2 - 1)) * w + 0.5 * w ^ 2, -0.5v + v ^ 2 + 0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(2.0 + -2.0 * (2 * u ^ (2 - 1)) + -2.0v + -2.0w + 2.0 * ((2 * u ^ (2 - 1)) * w) + 2.0 * (v * w), -2.0u + 2.0 * (u * w), -2.0u + 2.0 * u ^ 2 + 2.0 * (u * v)))), Vec(float.(tuple(2.0v + -2.0 * (v * w), 2.0u + -2.0 * (u * w), -2.0 * (u * v)))), Vec(float.(tuple(-2.0v + 2.0 * (v * w), 2.0 + -2.0 * (2 * v ^ (2 - 1)) + -2.0u + -2.0w + 2.0 * ((2 * v ^ (2 - 1)) * w) + 2.0 * (u * w), -2.0v + 2.0 * v ^ 2 + 2.0 * (u * v)))), Vec(float.(tuple(-1.0 + w ^ 2, -1.0 + w ^ 2, -1.0 * (2 * w ^ (2 - 1)) + u * (2 * w ^ (2 - 1)) + v * (2 * w ^ (2 - 1))))), Vec(float.(tuple(1 + -1.0 * w ^ 2, 0, -1.0 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, 1 + -1.0 * w ^ 2, -1.0 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(2.0 + -2.0 * (2 * u ^ (2 - 1)) + -2.0v + 2.0w + -2.0 * ((2 * u ^ (2 - 1)) * w) + -2.0 * (v * w), -2.0u + -2.0 * (u * w), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v)))), Vec(float.(tuple(2.0v + 2.0 * (v * w), 2.0u + 2.0 * (u * w), 2.0 * (u * v)))), Vec(float.(tuple(-2.0v + -2.0 * (v * w), 2.0 + -2.0 * (2 * v ^ (2 - 1)) + -2.0u + 2.0w + -2.0 * ((2 * v ^ (2 - 1)) * w) + -2.0 * (u * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v))))) + end + @inline function eval_dbasis!(::Lagrange{Wedge, 2}, xi::Vec) (u, v, w) = xi @inbounds return (Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 2.0v + 1.5w + -1.0 * ((2 * u ^ (2 - 1)) * w) + -2.0 * (v * w) + -0.5 * w ^ 2, -1.0 + 2 * v ^ (2 - 1) + 2.0u + 1.5w + -1.0 * ((2 * v ^ (2 - 1)) * w) + -2.0 * (u * w) + -0.5 * w ^ 2, -0.5 + 0.5 * (2 * w ^ (2 - 1)) + 1.5u + 1.5v + -1.0 * u ^ 2 + -1.0 * v ^ 2 + -2.0 * (u * v) + -0.5 * (u * (2 * w ^ (2 - 1))) + -0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 0.5w + -1.0 * ((2 * u ^ (2 - 1)) * w) + 0.5 * w ^ 2, 0, 0.5u + -1.0 * u ^ 2 + 0.5 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, -1.0 + 2 * v ^ (2 - 1) + 0.5w + -1.0 * ((2 * v ^ (2 - 1)) * w) + 0.5 * w ^ 2, 0.5v + -1.0 * v ^ 2 + 0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + 2.0v + -1.5w + (2 * u ^ (2 - 1)) * w + 2.0 * (v * w) + -0.5 * w ^ 2, -1.0 + 2 * v ^ (2 - 1) + 2.0u + -1.5w + (2 * v ^ (2 - 1)) * w + 2.0 * (u * w) + -0.5 * w ^ 2, 0.5 + 0.5 * (2 * w ^ (2 - 1)) + -1.5u + -1.5v + u ^ 2 + v ^ 2 + 2.0 * (u * v) + -0.5 * (u * (2 * w ^ (2 - 1))) + -0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(-1.0 + 2 * u ^ (2 - 1) + -0.5w + (2 * u ^ (2 - 1)) * w + 0.5 * w ^ 2, 0, -0.5u + u ^ 2 + 0.5 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, -1.0 + 2 * v ^ (2 - 1) + -0.5w + (2 * v ^ (2 - 1)) * w + 0.5 * w ^ 2, -0.5v + v ^ 2 + 0.5 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(2.0 + -2.0 * (2 * u ^ (2 - 1)) + -2.0v + -2.0w + 2.0 * ((2 * u ^ (2 - 1)) * w) + 2.0 * (v * w), -2.0u + 2.0 * (u * w), -2.0u + 2.0 * u ^ 2 + 2.0 * (u * v)))), Vec(float.(tuple(2.0v + -2.0 * (v * w), 2.0u + -2.0 * (u * w), -2.0 * (u * v)))), Vec(float.(tuple(-2.0v + 2.0 * (v * w), 2.0 + -2.0 * (2 * v ^ (2 - 1)) + -2.0u + -2.0w + 2.0 * ((2 * v ^ (2 - 1)) * w) + 2.0 * (u * w), -2.0v + 2.0 * v ^ 2 + 2.0 * (u * v)))), Vec(float.(tuple(-1.0 + w ^ 2, -1.0 + w ^ 2, -1.0 * (2 * w ^ (2 - 1)) + u * (2 * w ^ (2 - 1)) + v * (2 * w ^ (2 - 1))))), Vec(float.(tuple(1 + -1.0 * w ^ 2, 0, -1.0 * (u * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(0, 1 + -1.0 * w ^ 2, -1.0 * (v * (2 * w ^ (2 - 1)))))), Vec(float.(tuple(2.0 + -2.0 * (2 * u ^ (2 - 1)) + -2.0v + 2.0w + -2.0 * ((2 * u ^ (2 - 1)) * w) + -2.0 * (v * w), -2.0u + -2.0 * (u * w), 2.0u + -2.0 * u ^ 2 + -2.0 * (u * v)))), Vec(float.(tuple(2.0v + 2.0 * (v * w), 2.0u + 2.0 * (u * w), 2.0 * (u * v)))), Vec(float.(tuple(-2.0v + -2.0 * (v * w), 2.0 + -2.0 * (2 * v ^ (2 - 1)) + -2.0u + 2.0w + -2.0 * ((2 * v ^ (2 - 1)) * w) + -2.0 * (u * w), 2.0v + -2.0 * v ^ 2 + -2.0 * (u * v))))) end