Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 17:22:10 +00:00
Code Issues Packages Projects Releases Wiki Activity

feat: Copy FEMBasis.jl files to src/basis/ (Phase 1 start)

Browse Source 
- 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
This commit is contained in:
Jukka Aho
2025-11-08 09:02:40 +02:00
parent 724ed52923
commit 7b6fcfdaf5
16 changed files with 1042 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/basis/abstract.jl
+73
View File
@@ -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
src/basis/create_basis.jl
+119
View File
@@ -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...))
src/basis/lagrange_hexahedrons.jl
+81
View File
@@ -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),
)
src/basis/lagrange_pyramids.jl
+42
View File
@@ -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)),
],
)
src/basis/lagrange_quadrangles.jl
+47
View File
@@ -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),
)
src/basis/lagrange_segments.jl
+21
View File
@@ -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))
src/basis/lagrange_tetrahedrons.jl
+32
View File
@@ -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),
)
src/basis/lagrange_triangles.jl
+42
View File
@@ -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),
)
src/basis/lagrange_wedges.jl
+64
View File
@@ -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),
],
)
src/basis/math.jl
+255
View File
@@ -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
src/basis/nurbs.jl
+16
View File
@@ -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
src/basis/nurbs_segment.jl
+37
View File
@@ -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
src/basis/nurbs_solid.jl
+59
View File
@@ -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
src/basis/nurbs_surface.jl
+47
View File
@@ -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
src/basis/subs.jl
+54
View File
@@ -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
src/basis/vandermonde.jl
+53
View File
@@ -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