From 7b6fcfdaf5255357b4e89386a71e4fac88363248 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 8 Nov 2025 09:02:40 +0200 Subject: [PATCH] feat: Copy FEMBasis.jl files to src/basis/ (Phase 1 start) - Create src/basis/ directory structure - Copy all FEMBasis.jl source files verbatim: - abstract.jl: AbstractBasis type definition and interface - create_basis.jl: Metaprogramming for basis generation - lagrange_*.jl: All Lagrange element bases (Seg, Quad, Tri, Tet, Hex, Wedge, Pyr) - nurbs*.jl: NURBS basis functions - math.jl: jacobian, grad, interpolate functions - subs.jl, vandermonde.jl: Symbolic/mathematical utilities Strategy: Copy first, integrate later (safest approach) Next: Integrate into src/JuliaFEM.jl module --- src/basis/abstract.jl | 73 +++++++++ src/basis/create_basis.jl | 119 ++++++++++++++ src/basis/lagrange_hexahedrons.jl | 81 +++++++++ src/basis/lagrange_pyramids.jl | 42 +++++ src/basis/lagrange_quadrangles.jl | 47 ++++++ src/basis/lagrange_segments.jl | 21 +++ src/basis/lagrange_tetrahedrons.jl | 32 ++++ src/basis/lagrange_triangles.jl | 42 +++++ src/basis/lagrange_wedges.jl | 64 ++++++++ src/basis/math.jl | 255 +++++++++++++++++++++++++++++ src/basis/nurbs.jl | 16 ++ src/basis/nurbs_segment.jl | 37 +++++ src/basis/nurbs_solid.jl | 59 +++++++ src/basis/nurbs_surface.jl | 47 ++++++ src/basis/subs.jl | 54 ++++++ src/basis/vandermonde.jl | 53 ++++++ 16 files changed, 1042 insertions(+) create mode 100644 src/basis/abstract.jl create mode 100644 src/basis/create_basis.jl create mode 100644 src/basis/lagrange_hexahedrons.jl create mode 100644 src/basis/lagrange_pyramids.jl create mode 100644 src/basis/lagrange_quadrangles.jl create mode 100644 src/basis/lagrange_segments.jl create mode 100644 src/basis/lagrange_tetrahedrons.jl create mode 100644 src/basis/lagrange_triangles.jl create mode 100644 src/basis/lagrange_wedges.jl create mode 100644 src/basis/math.jl create mode 100644 src/basis/nurbs.jl create mode 100644 src/basis/nurbs_segment.jl create mode 100644 src/basis/nurbs_solid.jl create mode 100644 src/basis/nurbs_surface.jl create mode 100644 src/basis/subs.jl create mode 100644 src/basis/vandermonde.jl diff --git a/src/basis/abstract.jl b/src/basis/abstract.jl new file mode 100644 index 0000000..c3cefab --- /dev/null +++ b/src/basis/abstract.jl @@ -0,0 +1,73 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE + +# AbstractBasis type and interface +# Consolidated from FEMBasis.jl package + +using Tensors +using LinearAlgebra +import Calculus + +# Re-export Vec for convenience (from Tensors.jl) +export Vec + +# Type alias for coordinate inputs (tuples or Vec) +const Vecish{N, T} = Union{NTuple{N, T}, Vec{N, T}} + +""" + AbstractBasis{dim} + +Abstract base type for all finite element basis functions. + +# Type parameter +- `dim`: Dimensionality of the reference element (1, 2, or 3) + +# Interface requirements + +Concrete basis types must implement: +- `Base.length(::Type{<:AbstractBasis})` - Number of basis functions +- `Base.size(::Type{<:AbstractBasis})` - (dim, n_basis) +- `get_reference_element_coordinates(::Type{<:AbstractBasis})` - Reference coordinates +- `eval_basis!(::Type{<:AbstractBasis}, N, xi)` - Evaluate basis functions +- `eval_dbasis!(::Type{<:AbstractBasis}, dN, xi)` - Evaluate basis derivatives + +# Example + +```julia +struct Seg2 <: AbstractBasis{1} end +length(Seg2) == 2 +size(Seg2) == (1, 2) +``` +""" +abstract type AbstractBasis{dim} end + +# Forward methods on instances to types +# This allows calling methods on both Seg2 and Seg2() +Base.length(B::T) where {T<:AbstractBasis} = length(T) +Base.size(B::T) where {T<:AbstractBasis} = size(T) +eval_basis!(B::T, N, xi) where {T<:AbstractBasis} = eval_basis!(T, N, xi) +eval_dbasis!(B::T, dN, xi) where {T<:AbstractBasis} = eval_dbasis!(T, dN, xi) + +# Allocating versions (convenience wrappers) +""" + eval_basis(basis::AbstractBasis{dim}, xi) -> Vector{Float64} + +Evaluate basis functions at point `xi`, allocating return vector. + +See also: [`eval_basis!`](@ref) for non-allocating version. +""" +eval_basis(B::AbstractBasis{dim}, xi) where {dim} = eval_basis!(B, zeros(length(B)), xi) + +""" + eval_dbasis(basis::AbstractBasis{dim}, xi) -> Vector{Vec{dim, Float64}} + +Evaluate basis function derivatives at point `xi`, allocating return vector. + +See also: [`eval_dbasis!`](@ref) for non-allocating version. +""" +eval_dbasis(B::AbstractBasis{dim}, xi) where {dim} = eval_dbasis!(B, zeros(Vec{dim}, length(B)), xi) + +# Declare interface functions (will be implemented by basis generator) +function get_reference_element_coordinates end +function eval_basis! end +function eval_dbasis! end diff --git a/src/basis/create_basis.jl b/src/basis/create_basis.jl new file mode 100644 index 0000000..8bd8356 --- /dev/null +++ b/src/basis/create_basis.jl @@ -0,0 +1,119 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +__precompile__(false) + +function get_reference_element_coordinates end +function eval_basis! end +function eval_dbasis! end + +function calculate_interpolation_polynomials(p, V) + basis = [] + first(p.args) == :+ || error("Use only summation between terms of polynomial") + args = p.args[2:end] + n = size(V, 1) + b = zeros(n) + for i in 1:n + fill!(b, 0.0) + b[i] = 1.0 + # TODO: Think about numerical stability with + # inverting Vandermonde matrix? + solution = V \ b + N = Expr(:call, :+) + for (ai, bi) in zip(solution, args) + isapprox(ai, 0.0) && continue + push!(N.args, Calculus.simplify( :( $ai * $bi ) )) + end + push!(basis, N) + end + return basis +end + +function calculate_interpolation_polynomial_derivatives(basis, D) + vars = [:u, :v, :w] + dbasis = Matrix(undef, D, length(basis)) + for (i, N) in enumerate(basis) + partial_derivatives = [] + for j in 1:D + dbasis[j, i] = Calculus.simplify(Calculus.differentiate(N, vars[j])) + end + end + return dbasis +end + +function create_basis(name, description, X::Vector{<:Vecish{D}}, p::Expr) where D + @debug "create basis given antsatz polynomial" name description X p + V = vandermonde_matrix(p, X) + basis = calculate_interpolation_polynomials(p, V) + return create_basis(name, description, X, basis) +end + +function create_basis(name, description, X::Vector{<:Vecish{D}}, basis::Vector) where D + @assert length(X) == length(basis) + @debug "create basis given basis functions" name description X basis + dbasis = calculate_interpolation_polynomial_derivatives(basis, D) + return create_basis(name, description, Vec.(X), basis, dbasis) +end + +function create_basis(name, description, X::Vector{<:Vecish{D, T}}, basis, dbasis) where {D, T} + N = length(X) + @debug "create basis given basis functions and derivatives" name description X basis dbasis + + Q = Expr(:block) + for i=1:N + push!(Q.args, :(N[$i] = $(basis[i]))) + end + + V = Expr(:block) + for i=1:N + push!(V.args, :(dN[$i] = Vec(float.(tuple($(dbasis[:, i]...)))))) + end + + if D == 1 + unpack = :((u,) = xi) + elseif D == 2 + unpack = :((u, v) = xi) + else + unpack = :((u, v, w) = xi) + end + + code = quote + struct $name <: FEMBasis.AbstractBasis{$D} + end + + Base.@pure function Base.size(::Type{$name}) + return ($D, $N) + end + + function Base.size(::Type{$name}, j::Int) + j == 1 && return $D + j == 2 && return $N + end + + Base.@pure function Base.length(::Type{$name}) + return $N + end + + function FEMBasis.get_reference_element_coordinates(::Type{$name}) + return $X + end + + @inline function FEMBasis.eval_basis!(::Type{$name}, N::Vector{<:Number}, xi::Vec) + @assert length(N) == $N + $unpack + @inbounds $Q + return N + end + + @inline function FEMBasis.eval_dbasis!(::Type{$name}, dN::Vector{<:Vec{$D}}, xi::Vec) + @assert length(dN) == $N + $unpack + @inbounds $V + return dN + end + end + return code +end + +create_basis_and_eval(args...) = eval(create_basis(args...)) + diff --git a/src/basis/lagrange_hexahedrons.jl b/src/basis/lagrange_hexahedrons.jl new file mode 100644 index 0000000..053f11e --- /dev/null +++ b/src/basis/lagrange_hexahedrons.jl @@ -0,0 +1,81 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +code = create_basis_and_eval( + :Hex8, + "8 node linear hexahedral element", + [ + (-1.0, -1.0, -1.0), # N1 + ( 1.0, -1.0, -1.0), # N2 + ( 1.0, 1.0, -1.0), # N3 + (-1.0, 1.0, -1.0), # N4 + (-1.0, -1.0, 1.0), # N5 + ( 1.0, -1.0, 1.0), # N6 + ( 1.0, 1.0, 1.0), # N7 + (-1.0, 1.0, 1.0), # N8 + ], + :(1 + u + v + w + u*v + v*w + w*u + u*v*w), + ) + +code = create_basis_and_eval( + :Hex20, + "20 node biquadratic hexahedral element", + [ + (-1.0, -1.0, -1.0), # N1 + ( 1.0, -1.0, -1.0), # N2 + ( 1.0, 1.0, -1.0), # N3 + (-1.0, 1.0, -1.0), # N4 + (-1.0, -1.0, 1.0), # N5 + ( 1.0, -1.0, 1.0), # N6 + ( 1.0, 1.0, 1.0), # N7 + (-1.0, 1.0, 1.0), # N8 + ( 0.0, -1.0, -1.0), # N9 + ( 1.0, 0.0, -1.0), # N10 + ( 0.0, 1.0, -1.0), # N11 + (-1.0, 0.0, -1.0), # N12 + (-1.0, -1.0, 0.0), # N13 + ( 1.0, -1.0, 0.0), # N14 + ( 1.0, 1.0, 0.0), # N15 + (-1.0, 1.0, 0.0), # N16 + ( 0.0, -1.0, 1.0), # N17 + ( 1.0, 0.0, 1.0), # N18 + ( 0.0, 1.0, 1.0), # N19 + (-1.0, 0.0, 1.0), # N20 + ], + :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2), + ) + +code = create_basis_and_eval( + :Hex27, + "27 node quadratic hexahedral element", + [ + (-1.0, -1.0, -1.0), # N1 + ( 1.0, -1.0, -1.0), # N2 + ( 1.0, 1.0, -1.0), # N3 + (-1.0, 1.0, -1.0), # N4 + (-1.0, -1.0, 1.0), # N5 + ( 1.0, -1.0, 1.0), # N6 + ( 1.0, 1.0, 1.0), # N7 + (-1.0, 1.0, 1.0), # N8 + ( 0.0, -1.0, -1.0), # N9 + ( 1.0, 0.0, -1.0), # N10 + ( 0.0, 1.0, -1.0), # N11 + (-1.0, 0.0, -1.0), # N12 + (-1.0, -1.0, 0.0), # N13 + ( 1.0, -1.0, 0.0), # N14 + ( 1.0, 1.0, 0.0), # N15 + (-1.0, 1.0, 0.0), # N16 + ( 0.0, -1.0, 1.0), # N17 + ( 1.0, 0.0, 1.0), # N18 + ( 0.0, 1.0, 1.0), # N19 + (-1.0, 0.0, 1.0), # N20 + ( 0.0, 0.0, -1.0), # N21 + ( 0.0, -1.0, 0.0), # N22 + ( 1.0, 0.0, 0.0), # N23 + ( 0.0, 1.0, 0.0), # N24 + (-1.0, 0.0, 0.0), # N25 + ( 0.0, 0.0, 1.0), # N26 + ( 0.0, 0.0, 0.0), # N27 + ], + :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2 + u^2*v^2 + v^2*w^2 + u^2*w^2 + u^2*v^2*w + u*v^2*w^2 + u^2*v*w^2 + u^2*v^2*w^2), + ) diff --git a/src/basis/lagrange_pyramids.jl b/src/basis/lagrange_pyramids.jl new file mode 100644 index 0000000..516e066 --- /dev/null +++ b/src/basis/lagrange_pyramids.jl @@ -0,0 +1,42 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +# Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015. +code = create_basis_and_eval( + :Pyr5A, + "5 node linear pyramid element", + [ + ( 1.0, 1.0, 0.0), # N1 + ( 1.0, -1.0, 0.0), # N2 + (-1.0, -1.0, 0.0), # N3 + (-1.0, 1.0, 0.0), # N4 + ( 0.0, 0.0, 1.0), # N5 + ], + [ + :(1/4*( (1+u)*(1+v) - w + u*v*w/(1-w) )), + :(1/4*( (1+u)*(1-v) - w + u*v*w/(1-w) )), + :(1/4*( (1-u)*(1-v) - w + u*v*w/(1-w) )), + :(1/4*( (1-u)*(1+v) - w + u*v*w/(1-w) )), + :(1.0*w), + ], + ) + +# source: Code Aster documentation? +code = create_basis_and_eval( + :Pyr5, + "5 node linear pyramid element", + [ + (-1.0, -1.0, -1.0), # N1 + ( 1.0, -1.0, -1.0), # N2 + ( 1.0, 1.0, -1.0), # N3 + (-1.0, 1.0, -1.0), # N4 + ( 0.0, 0.0, 1.0), # N5 + ], + [ + :(1/8 * (1-u) * (1-v) * (1-w)), + :(1/8 * (1+u) * (1-v) * (1-w)), + :(1/8 * (1+u) * (1+v) * (1-w)), + :(1/8 * (1-u) * (1+v) * (1-w)), + :(1/2 * (1+w)), + ], + ) diff --git a/src/basis/lagrange_quadrangles.jl b/src/basis/lagrange_quadrangles.jl new file mode 100644 index 0000000..927f8cc --- /dev/null +++ b/src/basis/lagrange_quadrangles.jl @@ -0,0 +1,47 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +code = create_basis_and_eval( + :Quad4, + "4 node linear quadrangle element", + [ + (-1.0, -1.0), # N1 + ( 1.0, -1.0), # N2 + ( 1.0, 1.0), # N3 + (-1.0, 1.0) # N4 + ], + :(1 + u + v + u*v), + ) + +code = create_basis_and_eval( + :Quad8, + "8 node quadratic quadrangle element (Serendip)", + [ + (-1.0, -1.0), # N1 + ( 1.0, -1.0), # N2 + ( 1.0, 1.0), # N3 + (-1.0, 1.0), # N4 + ( 0.0, -1.0), # N5 + ( 1.0, 0.0), # N6 + ( 0.0, 1.0), # N7 + (-1.0, 0.0) # N8 + ], + :(1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2), + ) + +code = create_basis_and_eval( + :Quad9, + "9 node quadratic quadrangle element", + [ + (-1.0, -1.0), # N1 + ( 1.0, -1.0), # N2 + ( 1.0, 1.0), # N3 + (-1.0, 1.0), # N4 + ( 0.0, -1.0), # N5 + ( 1.0, 0.0), # N6 + ( 0.0, 1.0), # N7 + (-1.0, 0.0), # N8 + ( 0.0, 0.0) # N9 + ], + :(1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2 + u^2*v^2), + ) diff --git a/src/basis/lagrange_segments.jl b/src/basis/lagrange_segments.jl new file mode 100644 index 0000000..6750d06 --- /dev/null +++ b/src/basis/lagrange_segments.jl @@ -0,0 +1,21 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +code = create_basis_and_eval( + :Seg2, + "2 node linear segment/line element", + [ + (-1.0,), + ( 1.0,), + ], + :(1 + u)) + +code = create_basis_and_eval( + :Seg3, + "3 node quadratic segment/line element", + [ + (-1.0,), + ( 1.0,), + ( 0.0,), + ], + :(1 + u + u^2)) diff --git a/src/basis/lagrange_tetrahedrons.jl b/src/basis/lagrange_tetrahedrons.jl new file mode 100644 index 0000000..efe803d --- /dev/null +++ b/src/basis/lagrange_tetrahedrons.jl @@ -0,0 +1,32 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +code = create_basis_and_eval( + :Tet4, + "4 node linear tetrahedral element", + [ + (0.0, 0.0, 0.0), # N1 + (1.0, 0.0, 0.0), # N2 + (0.0, 1.0, 0.0), # N3 + (0.0, 0.0, 1.0), # N4 + ], + :(1 + u + v + w), + ) + +code = create_basis_and_eval( + :Tet10, + "10 node quadratic tetrahedral element", + [ + (0.0, 0.0, 0.0), # N1 + (1.0, 0.0, 0.0), # N2 + (0.0, 1.0, 0.0), # N3 + (0.0, 0.0, 1.0), # N4 + (0.5, 0.0, 0.0), # N5 + (0.5, 0.5, 0.0), # N6 + (0.0, 0.5, 0.0), # N7 + (0.0, 0.0, 0.5), # N8 + (0.5, 0.0, 0.5), # N9 + (0.0, 0.5, 0.5), # N10 + ], + :(1 + u + v + w + u*v + v*w + w*u + u^2 + v^2 + w^2), + ) diff --git a/src/basis/lagrange_triangles.jl b/src/basis/lagrange_triangles.jl new file mode 100644 index 0000000..7a443dc --- /dev/null +++ b/src/basis/lagrange_triangles.jl @@ -0,0 +1,42 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +code = create_basis_and_eval( + :Tri3, + "3 node linear triangle element", + [ + (0.0, 0.0), # N1 + (1.0, 0.0), # N2 + (0.0, 1.0), # N3 + ], + :(1 + u + v), + ) + +code = create_basis_and_eval( + :Tri6, + "6 node quadratic triangle element", + [ + (0.0, 0.0), # N1 + (1.0, 0.0), # N2 + (0.0, 1.0), # N3 + (0.5, 0.0), # N4 + (0.5, 0.5), # N5 + (0.0, 0.5), # N6 + ], + :(1 + u + v + u^2 + u*v + v^2), + ) + +code = create_basis_and_eval( + :Tri7, + "7 node quadratic triangle element (has middle node)", + [ + (0.0, 0.0), # N1 + (1.0, 0.0), # N2 + (0.0, 1.0), # N3 + (0.5, 0.0), # N4 + (0.5, 0.5), # N5 + (0.0, 0.5), # N6 + (1/3, 1/3), # N7 + ], + :(1 + u + v + u^2 + u*v + v^2 + u^2*v^2), + ) diff --git a/src/basis/lagrange_wedges.jl b/src/basis/lagrange_wedges.jl new file mode 100644 index 0000000..c78cafc --- /dev/null +++ b/src/basis/lagrange_wedges.jl @@ -0,0 +1,64 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +# Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015. +create_basis_and_eval( + :Wedge6, + "6 node linear prismatic/wedge element", + [ + (0.0, 0.0, -1.0), # N1 + (1.0, 0.0, -1.0), # N2 + (0.0, 1.0, -1.0), # N3 + (0.0, 0.0, 1.0), # N4 + (1.0, 0.0, 1.0), # N5 + (0.0, 1.0, 1.0), # N6 + ], + [ + :(1/2 * (1-w) * (1-u-v)), # N1 + :(1/2 * (1-w) * u), # N2 + :(1/2 * (1-w) * v), # N3 + :(1/2 * (1+w) * (1-u-v)), # N4 + :(1/2 * (1+w) * u), # N5 + :(1/2 * (1+w) * v), # N6 + ], + ) + +# Basis functions are from ABAQUS theory manual +create_basis_and_eval( + :Wedge15, + "15 node quadratic prismatic/wedge element", + [ + (0.0, 0.0, -1.0), # N1 + (1.0, 0.0, -1.0), # N2 + (0.0, 1.0, -1.0), # N3 + (0.0, 0.0, 1.0), # N4 + (1.0, 0.0, 1.0), # N5 + (0.0, 1.0, 1.0), # N6 + (0.5, 0.0, -1.0), # N7 + (0.5, 0.5, -1.0), # N8 + (0.0, 0.5, -1.0), # N9 + (0.5, 0.0, 1.0), # N10 + (0.5, 0.5, 1.0), # N11 + (0.0, 0.5, 1.0), # N12 + (0.0, 0.0, 0.0), # N13 + (1.0, 0.0, 0.0), # N14 + (0.0, 1.0, 0.0), # N15 + ], + [ + :(u^2*w^2 - u^2*w + 2*u*v*w^2 - 2*u*v*w - 3*u*w^2/2 + 3*u*w/2 + v^2*w^2 - v^2*w - 3*v*w^2/2 + 3*v*w/2 + w^2/2 - w/2), + :(u^2*w^2 - u^2*w - u*w^2/2 + u*w/2), + :(v^2*w^2 - v^2*w - v*w^2/2 + v*w/2), + :(u^2*w^2 + u^2*w + 2*u*v*w^2 + 2*u*v*w - 3*u*w^2/2 - 3*u*w/2 + v^2*w^2 + v^2*w - 3*v*w^2/2 - 3*v*w/2 + w^2/2 + w/2), + :(u^2*w^2 + u^2*w - u*w^2/2 - u*w/2), + :(v^2*w^2 + v^2*w - v*w^2/2 - v*w/2), + :(-2*u^2*w^2 + 2*u^2*w - 2*u*v*w^2 + 2*u*v*w + 2*u*w^2 - 2*u*w), + :(2*u*v*w^2 - 2*u*v*w), + :(-2*u*v*w^2 + 2*u*v*w - 2*v^2*w^2 + 2*v^2*w + 2*v*w^2 - 2*v*w), + :(2*u^2*w^2 - 2*u^2*w - 4*u^2 + 2*u*v*w^2 - 2*u*v*w - 4*u*v - 2*u*w^2 + 2*u*w + 4*u), + :(-2*u*v*w^2 + 2*u*v*w + 4*u*v), + :(2*u*v*w^2 - 2*u*v*w - 4*u*v + 2*v^2*w^2 - 2*v^2*w - 4*v^2 - 2*v*w^2 + 2*v*w + 4*v), + :(-2*u^2*w^2 + 2*u^2 - 4*u*v*w^2 + 4*u*v + 3*u*w^2 - 3*u - 2*v^2*w^2 + 2*v^2 + 3*v*w^2 - 3*v - w^2 + 1), + :(-2*u^2*w^2 + 2*u^2 + u*w^2 - u), + :(-2*v^2*w^2 + 2*v^2 + v*w^2 - v), + ], + ) diff --git a/src/basis/math.jl b/src/basis/math.jl new file mode 100644 index 0000000..88bac36 --- /dev/null +++ b/src/basis/math.jl @@ -0,0 +1,255 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +""" + interpolate(B, T, xi) + +Given basis B, interpolate T at xi. + +# Example +```jldoctest +B = Quad4() +X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]) +T = [1.0, 2.0, 3.0, 4.0] +interpolate(B, T, Vec(0.0, 0.0)) + +# output + +2.5 +``` +""" +function interpolate(B::AbstractBasis{dim}, T::Vector, xi::Vec{dim}) where {dim} + N = eval_basis(B, xi) + return sum(b*t for (b, t) in zip(N, T)) +end + +""" + jacobian(B, X, xi) + +Given basis B, calculate jacobian at xi. + +# Example +```jldoctest +B = Quad4() +X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]) +jacobian(B, X, Vec((0.0, 0.0))) + +# output + +2×2 Tensor{2,2,Float64,4}: + 0.5 0.0 + 0.0 0.5 + +``` +""" +jacobian(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} = jacobian(B, X, xi, eval_dbasis(B, xi)) + +function jacobian(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}, dB::Vector{<:Vec{dim}}) where {dim} + @assert length(X) == length(B) == length(dB) + J = zero(Tensor{2, dim}) + @inbounds for i in 1:length(X) + J += otimes(dB[i], X[i]) # dB[i] ⊗ X[i] + end + return J +end + + + +""" + grad(B, X, xi) + +Given basis B, calculate gradient dB/dX at xi. + +# Example +```jldoctest +B = Quad4() +X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]) +grad(B, X, Vec(0.0, 0.0)) + +# output + +4-element Array{Tensor{1,2,Float64,2},1}: + [-0.5, -0.5] + [0.5, -0.5] + [0.5, 0.5] + [-0.5, 0.5] + +``` +""" +grad(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} = + grad!(B, similar(X), X, xi, eval_dbasis(B, xi)) + +function grad!(B::AbstractBasis{dim}, dN::Vector{<:Vec{dim}}, X::Vector{<:Vec{dim}}, xi::Vec{dim}, dB::Vector{<:Vec{dim}}) where {dim} + @assert length(dN) == length(dB) + J = jacobian(B, X, xi, dB) + @inbounds for i in 1:length(dN) + dN[i] = inv(J) ⋅ dB[i] + end + return dN +end + +""" + grad(B, T, X, xi) + +Calculate gradient of `T` with respect to `X` in point `xi` using basis `B`. + +# Example +```jldoctest +B = Quad4() +X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]) +u = Vec.([(0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0)]) +grad(B, u, X, Vec(0.0, 0.0)) + +# output + +julia> grad(B, u, X, Vec(0.0, 0.0)) +2×2 Tensor{2,2,Float64,4}: + 1.5 0.5 + 1.0 2.0 + +``` +""" +function grad(B::AbstractBasis{dim}, T::Vector{<:Vec{dim}}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} + G = grad(B, X, xi) # <- allocates + dTdX = sum(T[i] ⊗ G[i] for i=1:length(B)) + return dTdX +end +function grad(B::AbstractBasis{dim}, T::Vector{<:Number}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} + G = grad(B, X, xi) # <- allocates + dTdX = sum(T[i] * G[i] for i=1:length(B)) + return dTdX +end + + +""" +Data type for fast FEM. +""" +mutable struct BasisInfo{B<:AbstractBasis,dim, T, M} + N::Vector{T} + dN::Vector{Vec{dim, T}} + grad::Vector{Vec{dim, T}} + J::Tensor{2, dim, T, M} + invJ::Tensor{2, dim, T, M} + detJ::T + basis::Type{B} +end + +Base.length(B::BasisInfo{T}) where T<:AbstractBasis = length(T) +Base.size(B::BasisInfo{T}) where T<:AbstractBasis = size(T) + +""" +Initialization of data type `BasisInfo`. + +# Examples + +```jldoctest + +BasisInfo(Tri3) + +# output + +FEMBasis.BasisInfo{FEMBasis.Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0) + +``` + +""" +function BasisInfo(::Type{B}, T=Float64) where B <: AbstractBasis{dim} where dim + nbasis = length(B) + N = zeros(T, nbasis) + dN = zeros(Vec{dim, T}, nbasis) + grad = zeros(Vec{dim, T}, nbasis) + J = zero(Tensor{2, dim, T}) + invJ = zero(Tensor{2, dim, T}) + detJ = zero(T) + return BasisInfo(N, dN, grad, J, invJ, detJ, B) +end + +""" +Evaluate basis, gradient and so on for some point `xi`. + +# Examples + +```jldoctest + +b = BasisInfo(Quad4) +X = Vec.([(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0)]) +xi = Vec(0.0, 0.0) +eval_basis!(b, X, xi) + +# output + +BasisInfo{Quad4,2,Float64,4}([0.25, 0.25, 0.25, 0.25], Tensors.Tensor{1,2,Float64,2}[[-0.25, -0.25], [0.25, -0.25], [0.25, 0.25], [-0.25, 0.25]], Tensors.Tensor{1,2,Float64,2}[[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25, Quad4) + +``` +""" +function eval_basis!(bi::BasisInfo{B}, + X::Vector{<:Vec{dim}}, xi::Vec{dim}) where B <: AbstractBasis{dim} where dim + # evaluate basis and derivatives + eval_basis!(B, bi.N, xi) + eval_dbasis!(B, bi.dN, xi) + + # calculate Jacobian + bi.J = jacobian(B(), X, xi, bi.dN) + + # calculate determinant of Jacobian + gradient operator + + # TODO, fixup curve + manifold + # @assert dim[1] == dim[2] + bi.invJ = inv(bi.J) + @inbounds for i in 1:length(bi.dN) + bi.grad[i] = bi.invJ ⋅ bi.dN[i] + end + bi.detJ = det(bi.J) + #= + elseif dim1 == 1 # curve + bi.detJ = norm(bi.J) + elseif dim1 == 2 # manifold + bi.detJ = norm(cross(bi.J[1,:], bi.J[2,:])) + end + =# + + return bi +end + +""" + grad!(bi, gradu, u) + +Evalute gradient ∂u/∂X and store result to matrix `gradu`. It is assumed +that `eval_basis!` has been already run to `bi` so it already contains +all necessary matrices evaluated with some `X` and `xi`. + +# Example + +First setup and evaluate basis using `eval_basis!`: +```jldoctest ex1 +B = BasisInfo(Quad4) +X = Vec.([(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0)]) +xi = Vec(0.0, 0.0) +eval_basis!(B, X, xi) + +# output + +BasisInfo{Quad4,2,Float64}([0.25, 0.25, 0.25, 0.25], Tensors.Tensor{1,2,Float64,2}[[-0.25, -0.25], [0.25, -0.25], [0.25, 0.25], [-0.25, 0.25]], Tensors.Tensor{1,2,Float64,2}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.5, 0.5]], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25) + +``` + +Next, calculate gradient of `u`: +```jldoctest ex1 +u = Vec.([(0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0)]) +grad(B, u) + +# output + +2×2 Tensors.Tensor{2,2,Float64,4}: + 1.5 0.5 + 1.0 2.0 + +``` +""" +function grad(bi::BasisInfo{B}, u::Vector{<:Vec{dim}}) where B <: AbstractBasis{dim} where dim + gradu = zero(Tensor{2, dim}) + for k in 1:length(B) + gradu += otimes(u[k], bi.grad[k]) + end + return gradu +end diff --git a/src/basis/nurbs.jl b/src/basis/nurbs.jl new file mode 100644 index 0000000..f104f53 --- /dev/null +++ b/src/basis/nurbs.jl @@ -0,0 +1,16 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +import Base: size, length + +function NURBS(i, p, u, t) + p == 0 && return t[i] <= u <= t[i+1] ? 1.0 : 0.0 + anom = u-t[i] + adenom = t[i+p]-t[i] + a = isapprox(adenom, 0.0) ? 0.0 : anom/adenom + bnom = t[i+p+1]-u + bdenom = t[i+p+1]-t[i+1] + b = isapprox(bdenom, 0.0) ? 0.0 : bnom/bdenom + result = a*NURBS(i,p-1,u,t) + b*NURBS(i+1,p-1,u,t) + return result +end diff --git a/src/basis/nurbs_segment.jl b/src/basis/nurbs_segment.jl new file mode 100644 index 0000000..06ac906 --- /dev/null +++ b/src/basis/nurbs_segment.jl @@ -0,0 +1,37 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +""" NURBS segment. """ +mutable struct NSeg <: AbstractBasis{1} + order :: Int + knots :: Vector{Float64} + weights :: Vector{Float64} +end + +function NSeg() + NSeg(1, + [-1.0, -1.0, 1.0, 1.0], + ones(4)) +end + +function length(basis::NSeg) + nu = length(basis.knots) - basis.order - 1 + return nu +end + +function size(basis::NSeg) + return (1, length(basis)) +end + +function eval_basis!(basis::NSeg, N::Vector, xi::Vec{1}) + pu = basis.order + tu = basis.knots + w = basis.weights + nu = length(tu)-pu-1 + u = xi[1] + for j=1:nu + N[j] = w[j]*NURBS(j,pu,u,tu) + end + N ./= sum(N) + return N +end diff --git a/src/basis/nurbs_solid.jl b/src/basis/nurbs_solid.jl new file mode 100644 index 0000000..6a4730d --- /dev/null +++ b/src/basis/nurbs_solid.jl @@ -0,0 +1,59 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +mutable struct NSolid <: AbstractBasis{3} + order_u :: Int + order_v :: Int + order_w :: Int + knots_u :: Vector{Float64} + knots_v :: Vector{Float64} + knots_w :: Vector{Float64} + weights :: Array{Float64, 3} +end + +function NSolid() + NSolid(1, 1, 1, + [-1.0, -1.0, 1.0, 1.0], + [-1.0, -1.0, 1.0, 1.0], + [-1.0, -1.0, 1.0, 1.0], + ones(2, 2, 2)) +end + +function length(basis::NSolid) + nu = length(basis.knots_u) - basis.order_u - 1 + nv = length(basis.knots_v) - basis.order_v - 1 + nw = length(basis.knots_w) - basis.order_w - 1 + return nu*nv*nw +end + +function size(basis::NSolid) + return (3, length(basis)) +end + +function eval_basis!(basis::NSolid, N::Vector, xi::Vec{3}) + pu = basis.order_u + pv = basis.order_v + pw = basis.order_w + tu = basis.knots_u + tv = basis.knots_v + tw = basis.knots_w + weights = basis.weights + nu = length(tu)-pu-1 + nv = length(tv)-pv-1 + nw = length(tw)-pw-1 + u, v, w = xi + n = 1 + for i=1:nu + for j=1:nv + for k=1:nw + A = NURBS(i,pu,u,tu) + B = NURBS(j,pv,v,tv) + C = NURBS(k,pw,w,tw) + N[n] = weights[i,j,k]*A*B*C + n += 1 + end + end + end + N ./= sum(N) + return N +end diff --git a/src/basis/nurbs_surface.jl b/src/basis/nurbs_surface.jl new file mode 100644 index 0000000..bf20abb --- /dev/null +++ b/src/basis/nurbs_surface.jl @@ -0,0 +1,47 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +mutable struct NSurf <: AbstractBasis{2} + order_u :: Int + order_v :: Int + knots_u :: Vector{Float64} + knots_v :: Vector{Float64} + weights :: Matrix{Float64} +end + +function NSurf() + NSurf(1, 1, + [-1.0, -1.0, 1.0, 1.0], + [-1.0, -1.0, 1.0, 1.0], + ones(2, 2)) +end + +function length(basis::NSurf) + nu = length(basis.knots_u) - basis.order_u - 1 + nv = length(basis.knots_v) - basis.order_v - 1 + return nu*nv +end + +function size(basis::NSurf) + return (2, length(basis)) +end + +function eval_basis!(basis::NSurf, N::Vector, xi::Vec{2}) + pu = basis.order_u + pv = basis.order_v + tu = basis.knots_u + tv = basis.knots_v + w = basis.weights + nu = length(tu)-pu-1 + nv = length(tv)-pv-1 + u, v = xi + n = 1 + for i=1:nu + for j=1:nv + N[n] = w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv) + n += 1 + end + end + N ./= sum(N) + return N +end diff --git a/src/basis/subs.jl b/src/basis/subs.jl new file mode 100644 index 0000000..a4aa6e2 --- /dev/null +++ b/src/basis/subs.jl @@ -0,0 +1,54 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +function subs(p::Number, ::Any) + return p +end + +function subs(p::Symbol, data::Pair{Symbol, T}) where T + k, v = data + if p == k + return v + end + return p +end + +function subs(p::Symbol, data::NTuple{N,Pair{Symbol, T}}) where {N, T} + for (k, v) in data + if p == k + return v + end + end + return p +end + +function subs(p::Expr, d::Pair) + v = copy(p) + for j in 2:length(p.args) + v.args[j] = subs(v.args[j], d) + end + return v +end + +""" + subs(expression, data) + +Given expression and pair(s) of `symbol => value` data, substitute to expression. + +# Examples + +Let us have polynomial `1 + u + v + u*v^2`, and substitute u=1 and v=2: +```julia +expression = :(1 + u + v + u*v^2) +data = (:u => 1.0, :v => 2.0) +subs(expression, data) +8.0 +``` + +""" +function subs(p::Expr, data::NTuple{N,Pair{Symbol, T}}) where {N, T} + for di in data + p = subs(p, di) + end + return Calculus.simplify(p) +end diff --git a/src/basis/vandermonde.jl b/src/basis/vandermonde.jl new file mode 100644 index 0000000..254d1b2 --- /dev/null +++ b/src/basis/vandermonde.jl @@ -0,0 +1,53 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE + +""" + vandermonde_matrix(polynomial, coordinates) + +Given some polynomial and coordinates points (1-3 dimensions), create a Vandermonde +matrix. + +# Example + +To genererate a Vandermonde matrix for a reference quadrangle `[-1.0, 1.0]^2` for +polynomial `p(u,v) = 1 + u + v + u*v`, one writes: + +```julia +polynomial = :(1 + u + v + u*v) +coordinates = [(-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0)] +V = vandermonde_matrix(polynomial, coordinates) + +# output + +[ +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 +] + +``` + +# References +- Wikipedia contributors. (2018, August 1). Vandermonde matrix. In Wikipedia, The Free Encyclopedia. Retrieved 10:00, August 20, 2018, from https://en.wikipedia.org/w/index.php?title=Vandermonde_matrix&oldid=852930962 +""" +function vandermonde_matrix(polynomial::Expr, coordinates::Vector{NTuple{D, T}}) where {D, T<:Number} + N = length(coordinates) + A = zeros(N, N) + first(polynomial.args) == :+ || error("Use only summation between terms of polynomial") + args = polynomial.args[2:end] + for i in 1:N + X = coordinates[i] + if D == 1 + data = (:u => X[1],) + elseif D == 2 + data = (:u => X[1], :v => X[2]) + elseif D == 3 + data = (:u => X[1], :v => X[2], :w => X[3]) + end + for (j, term) in enumerate(args) + A[i,j] = subs(term, data) + end + end + return A +end