6  Finite elements

Variables and helper function

@variables a b x y l EA EI ν q

function prettyprint(name, factor, value)
    s1 = name.s[2:end-1]
    s2 = latexify(simplify(factor), env=:raw)
    s3 = latexify(simplify.(1 / factor * value, expand=true), env=:raw)
    return L"%$s1 = %$s2 %$s3"
end;

6.1 1D basis functions

φs = hatfunctions([0, 0.25, 0.75, 1.5, 2.0])

fig = Figure()
fplot!(Axis(fig[1, 1], title="Functions"), φs...)
fplot!(Axis(fig[2, 1], title="Derivatives"), derivative.(φs)...)
fig

6.2 Elements

fig = Figure(size=(700, 700))
Axis(fig)
feplot(fig[1, 1], makeelement(:lagrange, IHat, k=1)) |> feconf()
feplot(fig[1, 2], makeelement(:lagrange, IHat, k=3)) |> feconf()
feplot(fig[2, 1], makeelement(:hermite, IHat)) |> feconf()
feplot(fig[2, 2], makeelement(:lagrange, QHat, k=1)) |> feconf()
feplot(fig[3, 1], makeelement(:lagrange, QHat, k=2)) |> feconf()
feplot(fig[3, 2], makeelement(:serendipity, QHat, k=2)) |> feconf()
feplot(fig[4, 1], makeelement(:hermite, QHat, conforming=true)) |> feconf()
feplot(fig[4, 2], makeelement(:hermite, QHat, conforming=false)) |> feconf()
fig

6.3 Shape functions

Linear shape functions

makeelement(:lagrange, IHat, k=1) |> nodalbasis |> components |> fplot

Hermite polynomials

Ns = nodalbasis(makeelement(:hermite, 1 .. 4)) |> components
fig = Figure()
fplot!(Axis(fig[1, 1], title="Functions"), Ns...)
fplot!(Axis(fig[2, 1], title="Derivatives"), derivative.(Ns)...)
fig

Bilinear shape functions

makeelement(:lagrange, QHat, k=1) |> nodalbasis |> components |> fplot3d

Serendipity shape functions

makeelement(:serendipity, (0 .. 2) × (0 .. 1), k=2) |> nodalbasis |> components |> fplot3d

Conforming rectangular plate element shape functions

makeelement(:hermite, QHat) |> nodalbasis |> components |> fplot3d

Quadratic serendipity shape functions on the domain \([0, a] \times [0, b]\)

K = (0 .. a) × (0 .. b)
Ns = nodalbasis(makeelement(:serendipity, K, k=2))

prettyprint(L"\mathbf{N}(x, y)", 1 / (a^2 * b^2), [Nᵢ(x, y) for Nᵢ ∈ Ns])

\(\mathbf{N}(x, y) = \frac{1}{b^{2} a^{2}} \left[ \begin{array}{c} b^{2} a^{2} - 3 a^{2} b y + 2 y^{2} a^{2} - 3 b^{2} a x + 5 a b x y - 2 y^{2} a x + 2 x^{2} b^{2} - 2 x^{2} b y \\ - b^{2} a x - a b x y + 2 y^{2} a x + 2 x^{2} b^{2} - 2 x^{2} b y \\ - 3 a b x y + 2 y^{2} a x + 2 x^{2} b y \\ - a^{2} b y + 2 y^{2} a^{2} - a b x y - 2 y^{2} a x + 2 x^{2} b y \\ 4 b^{2} a x - 4 a b x y - 4 x^{2} b^{2} + 4 x^{2} b y \\ 4 a b x y - 4 y^{2} a x \\ 4 a b x y - 4 x^{2} b y \\ 4 a^{2} b y - 4 y^{2} a^{2} - 4 a b x y + 4 y^{2} a x \\ \end{array} \right]\)

6.4 Finite element stiffness matrices

6.4.1 Truss

Ns = nodalbasis(makeelement(:lagrange, IHat, k=1))
D(u) = 2 / l * u'
aᵉ(u, δu) = EA * l / 2 * integrate(D(u) * D(δu), IHat)
Kᵉ = [aᵉ(Nᵢ, Nⱼ) for Nᵢ ∈ Ns, Nⱼ ∈ Ns]

prettyprint(L"\mathbf{K}^e", EA / l, Kᵉ)

\(\mathbf{K}^e = \frac{\mathtt{EA}}{l} \left[ \begin{array}{cc} 1 & -1 \\ -1 & 1 \\ \end{array} \right]\)

6.4.2 Beam

Ns = [1, l / 2, 1, l / 2] .* nodalbasis(makeelement(:hermite, IHat))
D(w) = 4 / l^2 * w''
aᵉ(w, δw) = EI * l / 2 * integrate(D(w) * D(δw), IHat)
Kᵉ = [aᵉ(Nᵢ, Nⱼ) for Nᵢ ∈ Ns, Nⱼ ∈ Ns]

prettyprint(L"\mathbf{K}^e", EI / l^3, Kᵉ)

\(\mathbf{K}^e = \frac{\mathtt{EI}}{l^{3}} \left[ \begin{array}{cccc} 12 & 6 l & -12 & 6 l \\ 6 l & 4 l^{2} & - 6 l & 2 l^{2} \\ -12 & - 6 l & 12 & - 6 l \\ 6 l & 2 l^{2} & - 6 l & 4 l^{2} \\ \end{array} \right]\)

bᵉ(δw) = q * l / 2 * integrate(δw, IHat)
rᵉ = [bᵉ(Nᵢ) for Nᵢ ∈ Ns]

prettyprint(L"\mathbf{r}^e", q / 12, rᵉ)

\(\mathbf{r}^e = \frac{1}{12} q \left[ \begin{array}{c} 6 l \\ l^{2} \\ 6 l \\ - 1 l^{2} \\ \end{array} \right]\)

6.4.3 2D Poisson equation

Integration on the reference quadrilateral

Ns = nodalbasis(makeelement(:lagrange, QHat, k=1))
D(u) = [2 / a * ∂x(u), 2 / b * ∂y(u)]
aᵉ(u, δu) = a * b / 4 * integrate(D(u) ⋅ D(δu), QHat)
Kᵉ = [aᵉ(Nᵢ, Nⱼ) for Nᵢ ∈ Ns, Nⱼ ∈ Ns]

prettyprint(L"\mathbf{K}^e", 1 / (6 * a * b), Kᵉ)

\(\mathbf{K}^e = \frac{1}{6 a b} \left[ \begin{array}{cccc} 2 \left( a^{2} + b^{2} \right) & a^{2} - 2 b^{2} & - a^{2} - b^{2} & - 2 a^{2} + b^{2} \\ a^{2} - 2 b^{2} & 2 \left( a^{2} + b^{2} \right) & - 2 a^{2} + b^{2} & - a^{2} - b^{2} \\ - a^{2} - b^{2} & - 2 a^{2} + b^{2} & 2 \left( a^{2} + b^{2} \right) & a^{2} - 2 b^{2} \\ - 2 a^{2} + b^{2} & - a^{2} - b^{2} & a^{2} - 2 b^{2} & 2 \left( a^{2} + b^{2} \right) \\ \end{array} \right]\)

Alternative derivation in physical coordinates

K = (0 .. a) × (0 .. b)
Ns = nodalbasis(makeelement(:lagrange, K, k=1))
aᵉ(u, δu) = integrate(∇(u) ⋅ ∇(δu), K)
Kᵉ = [aᵉ(Nᵢ, Nⱼ) for Nᵢ ∈ Ns, Nⱼ ∈ Ns]

prettyprint(L"\mathbf{K}^e", 1 / (6 * a * b), Kᵉ)

\(\mathbf{K}^e = \frac{1}{6 a b} \left[ \begin{array}{cccc} 2 \left( a^{2} + b^{2} \right) & a^{2} - 2 b^{2} & - a^{2} - b^{2} & - 2 a^{2} + b^{2} \\ a^{2} - 2 b^{2} & 2 \left( a^{2} + b^{2} \right) & - 2 a^{2} + b^{2} & - a^{2} - b^{2} \\ - a^{2} - b^{2} & - 2 a^{2} + b^{2} & 2 \left( a^{2} + b^{2} \right) & a^{2} - 2 b^{2} \\ - 2 a^{2} + b^{2} & - a^{2} - b^{2} & a^{2} - 2 b^{2} & 2 \left( a^{2} + b^{2} \right) \\ \end{array} \right]\)

bᵉ(δu) = integrate(δu, K)
rᵉ = [bᵉ(Nᵢ) for Nᵢ ∈ Ns]

prettyprint(L"\mathbf{r}^e", a * b / 4, rᵉ)

\(\mathbf{r}^e = \frac{1}{4} a b \left[ \begin{array}{c} 1 \\ 1 \\ 1 \\ 1 \\ \end{array} \right]\)