Skip to content
FEMBasis contains interpolation routines for finite element function spaces. Given ansatz and coordinates of domain, shape functions are calculated symbolically in a very general way to get efficient code. Shape functions can also be given directly and in that case partial derivatives are calculated automatically.
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
docs
src
test
.appveyor.yml
.gitignore
.travis.yml
LICENSE
README.md
REQUIRE

README.md

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,
    "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),
   )

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])
grad(B, u, X, (0.0, 0.0))

Result is

2×2 Array{Float64,2}:
 1.5  0.5
 1.0  2.0
You can’t perform that action at this time.