Introduction

`FEMBasis.jl` contains interpolation routines for standard finite element function spaces. Given ansatz and coordinates of domain, interpolation functions are calculated symbolically in a very general way to get efficient code. As a concrete example, to generate basis functions for a standard 10-node tetrahedron one can write

``````code = create_basis(
:Tet10,
(
(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),
)``````

The resulting code is

``````begin
mutable struct Tet10
end
function Base.size(::Type{Tet10})
return (3, 10)
end
function Base.size(::Type{Tet10}, j::Int)
j == 1 && return 3
j == 2 && return 10
end
function Base.length(::Type{Tet10})
return 10
end
function FEMBasis.get_reference_element_coordinates(::Type{Tet10})
return ((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))
end
function FEMBasis.eval_basis!(::Type{Tet10}, N::Matrix{T}, xi::Tuple{T, T, T}) where T
(u, v, w) = xi
begin
N[1] = 1.0 + -3.0u + -3.0v + -3.0w + 4.0 * (u * v) + 4.0 * (v * w) + 4.0 * (w * u) + 2.0 * u ^ 2 + 2.0 * v ^ 2 + 2.0 * w ^ 2
N[2] = -u + 2.0 * u ^ 2
N[3] = -v + 2.0 * v ^ 2
N[4] = -w + 2.0 * w ^ 2
N[5] = 4.0u + -4.0 * (u * v) + -4.0 * (w * u) + -4.0 * u ^ 2
N[6] = +(4.0 * (u * v))
N[7] = 4.0v + -4.0 * (u * v) + -4.0 * (v * w) + -4.0 * v ^ 2
N[8] = 4.0w + -4.0 * (v * w) + -4.0 * (w * u) + -4.0 * w ^ 2
N[9] = +(4.0 * (w * u))
N[10] = +(4.0 * (v * w))
end
return N
end
function FEMBasis.eval_dbasis!(::Type{Tet10}, dN::Matrix{T}, xi::Tuple{T, T, T}) where T
(u, v, w) = xi
begin
dN[1, 1] = -3.0 + 4.0v + 4.0w + 2.0 * (2u)
dN[2, 1] = -3.0 + 4.0u + 4.0w + 2.0 * (2v)
dN[3, 1] = -3.0 + 4.0v + 4.0u + 2.0 * (2w)
dN[1, 2] = -1 + 2.0 * (2u)
dN[2, 2] = 0
dN[3, 2] = 0
dN[1, 3] = 0
dN[2, 3] = -1 + 2.0 * (2v)
dN[3, 3] = 0
dN[1, 4] = 0
dN[2, 4] = 0
dN[3, 4] = -1 + 2.0 * (2w)
dN[1, 5] = 4.0 + -4.0v + -4.0w + -4.0 * (2u)
dN[2, 5] = -4.0u
dN[3, 5] = -4.0u
dN[1, 6] = 4.0v
dN[2, 6] = 4.0u
dN[3, 6] = 0
dN[1, 7] = -4.0v
dN[2, 7] = 4.0 + -4.0u + -4.0w + -4.0 * (2v)
dN[3, 7] = -4.0v
dN[1, 8] = -4.0w
dN[2, 8] = -4.0w
dN[3, 8] = 4.0 + -4.0v + -4.0u + -4.0 * (2w)
dN[1, 9] = 4.0w
dN[2, 9] = 0
dN[3, 9] = 4.0u
dN[1, 10] = 0
dN[2, 10] = 4.0w
dN[3, 10] = 4.0v
end
return dN
end
end``````

Also more unusual elements can be defined. For example, pyramid element cannot be descibed with ansatz, but it's still possible to implement by defining shape functions, `Calculus.jl` is taking care of defining partial derivatives of function:

``````code = create_basis(
: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)),
),
)
eval(code)``````

Basis function can have internal variables if needed, e.g. variable dof basis like hierarchical basis functions or NURBS.

It's also possible to do some very common FEM calculations, like calculate Jacobian or gradient of some variable with respect to some coordinates. For example, to calculate displacement gradient du/dX in unit square [0,1]^2, one could write:

``````B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
``````2×2 Array{Float64,2}: