Finite Elements¶

CalcpadCE implements transparent finite-element method solvers in worksheet form — every stiffness matrix, load vector and post-processing step is plain code rather than a black-box routine.

The deep-beam plane-stress model discretises a 2D elastic continuum with rectangular 8-DOF elements and handles elastic supports through a distributed-stiffness function. The rectangular slab uses a 16-DOF Bogner–Fox–Schmit (BFS) plate element on a simply-supported plan, while the multi-span flat slab extends the same BFS element to irregular column grids with point supports. The analytically integrated variant swaps numerical Gauss integration for closed-form element matrices, cutting solve time on dense meshes by an order of magnitude.

Mesh density, geometry and material constants are parametric, so each worksheet doubles as a benchmark for convergence studies and for cross-checking commercial FEA results.

Deep Beam FEA¶

Plane-stress FEA of a deep beam using rectangular 8-DOF elements, with a partial uniform load and elastic supports modelled by a distributed-stiffness function. Recovers the \(σ_x\), \(σ_y\) and \(τ_{xy}\) stress fields and the deflected shape on a parametric mesh.

Code:

'<h3>Input data</h3>
'<table><tr><td style="padding:0;">
'Length -'l = 4'm, Height -'h = 2'm
'Thickness -'t = 0.1'm
'Distributed load -'q = 100'kN/m
'<p><b>Loads</b></p>
'Load length -'b = 0.8'm
'Load position -'a = (l - b)/2'm
'Load function -'q(x) = q*(x ≥ a ∧ x ≤ a + b)
'</td><td><img src="./deep-beam-fea.png" style="width:140; height:160pt; margin-left:20pt;"></td></tr></table>
'<p><b>Supports</b></p>
'Support length -'d = 0.4'm
'Support elastic stiffness -'k_s = 50000'MN/m²
'Support function -'s(x) = k_s*t*(x ≤ d ∨ x ≥ l - d)
'<p><b>Material properties</b></p>
'Modulus of elasticity -'E = 20000'MPa
'Poisson`s ratio -'ν = 0.2
'<h3>Finite element mesh</h3>
'We will use rectangular finite element with'n = 8'DOFs
'Number of elements along <var>l</var> and <var>h</var> -'n_l = 20', 'n_h = 10
'Total number of elements -'n_e = n_l*n_h
'Total number of joints -'n_j = (n_l + 1)*(n_h + 1)
'Element dimensions -'l_1 = l/n_l', 'h_1 = h/n_h
#hide
x_j = vector(n_j)','y_j = vector(n_j)
x = 0', 'y = 0
#for j = 1 : n_j
    x_j.j = x
    y_j.j = y
    y = y + h_1
    #if y > h
        y = 0
        x = x + l_1
    #end if
#loop
e_j = matrix(n_e; 4)
#for i_l = 1 : n_l
    #for i_h = 1 : n_h
        e = i_h + n_h*(i_l - 1)
        j = e + i_l - 1
        e_j.(e; 1) = j
        e_j.(e; 2) = j + 1
        e_j.(e; 3) = j + n_h + 1
        e_j.(e; 4) = j + n_h + 2
    #loop
#loop
#show
'Joint coordinates
x_j'm
y_j'm
'Joint numbers for elements
transp(e_j)
'Coordinates of element centers
j_e(e) = row(e_j; e)
x_c(e) = sum(extract(x_j; j_e(e)))/4',' _
y_c(e) = sum(extract(y_j; j_e(e)))/4
#hide
x_c = vector(n_e)', 'y_c = vector(n_e)
$Repeat{x_c.e = x_c(e) @ e = 1 : n_e}
$Repeat{y_c.e = y_c(e) @ e = 1 : n_e}
#show
x_c'm
y_c'm
'Elements along the bottom (supported) edge
e_S = find_lt(y_c; h_1; 1)
'Elements along the top (loaded) edge
e_L = find_gt(y_c; h - h_1; 1)
#val
#hide
W = 400
k = W/l
δ = 20
H = h/l*W
r = 0.02*k
#show
'<svg viewbox="'-δ/2' '-3*δ' 'W + 2*δ' 'H + 5*δ'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Segoe UI; font-size:8px; width:'W + δ'pt; height:'H + 5*δ'pt">
'<style>.joint{fill:orangeRed;} .support{stroke:red; stroke-width:1; fill:lightpink;} .element{stroke:seaGreen; stroke-width:1; fill:lime; fill-opacity:0.1; stroke-opacity:0.5;} .struct{fill:yellow; fill-opacity:0.25;stroke:goldenRod; stroke-opacity:0.5;}</style>
'<rect x="'0'" y="'0'" width="'W'" height="'H'" class="struct"/>
'<rect x="'0'" y="'H'" width="'d*k'" height="'1.5*δ'" class="struct"/>
'<rect x="'W - d*k'" y="'H'" width="'d*k'" height="'1.5*δ'" class="struct"/>
'<rect x="'a*k'" y="'-1.5*δ'" width="'b*k'" height="'1.5*δ'" style="fill:cyan; fill-opacity:0.15; stroke:darkCyan; stroke-opacity:0.5;"/>
'<text x="'(a + b/2)*k + r'" y="'-1.8*δ'" text-anchor="middle" style="font-size:11px;">q='q' kN/m</text>
#for e = 1 : n_e
    #hide
    x = x_c.e*k
    y = y_c.e*k
    #show
    '<text x="'x'" y="'H - y'" text-anchor="middle">'e'</text>
    '<rect x="'x - l_1/2*k'" y="'H - y - h_1/2*k'" width="'l_1*k'" height="'h_1*k'" class="element" />
#loop
#for j = 1 : n_j
    #hide
    x = x_j.j*k
    y = H - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'r'" class="joint" />
    '<text x="'x + r'" y="'y - r'" text-anchor="start">'j'</text>
#loop
'</svg>
#noc
'<h3>Finite element formation</h3>
'<p><b>Shape functions</b></p>
'They are defined in relative coordinates to the center of the element:
'<img src="./plane-element.png" alt="plane-elemen.png" style="float:right; width:150pt;" />
ξ = 2*(x - x_c)/l_1'∈ (-1; 1),'η = 2*(y - y_c)/h_1'∈ (-1; 1)
#equ
'Base functions -'Φ_1(ξ) = (1 - ξ)/2', 'Φ_2(ξ) = (1 + ξ)/2
'First derivatives -'Φ′_1(ξ) = -1/2', 'Φ′_2(ξ) = 1/2
'Shape functions for joints at elements` corners
N_1(ξ; η) = Φ_1(ξ)*Φ_1(η)', 'N_2(ξ; η) = Φ_1(ξ)*Φ_2(η)
N_3(ξ; η) = Φ_2(ξ)*Φ_1(η)', 'N_4(ξ; η) = Φ_2(ξ)*Φ_2(η)
'<p><b>Constitutive matrix</b> (stress-strain relationship)</p>
D = E*t/(1 - ν^2)*[1; ν; 0|ν; 1; 0|0; 0; (1 - ν)/2]
'<p><b>Strain-displacement matrix</b></p>
B_1(j; η) = (1/l_1)*take(j; -Φ_1(η); 0; -Φ_2(η); 0; Φ_1(η); 0; Φ_2(η); 0)
B_2(j; ξ) = (1/h_1)*take(j; 0; -Φ_1(ξ); 0; Φ_1(ξ); 0; -Φ_2(ξ); 0; Φ_2(ξ))
B_3(j; ξ; η) = take(j; _
-Φ_1(ξ)/h_1; -Φ_1(η)/l_1; Φ_1(ξ)/h_1; -Φ_2(η)/l_1; _
-Φ_2(ξ)/h_1; Φ_1(η)/l_1; Φ_2(ξ)/h_1; Φ_2(η)/l_1)
B(j; ξ; η) = [B_1(j; η); B_2(j; ξ); B_3(j; ξ; η)]
'The elements of the stiffness matrix will be calculated by using direct integration
#noc
K_e,ij = l_1*h_1/4*$Area{$Area{B_i(ξ; η)^T*D*B_j(ξ; η) @ ξ = -1 : 1} @ η = -1 : 1}
#equ
'<p><b>Element stiffness matrix</b> (above the main diagonal only)</p>
#hide
K_e = utriang(8)
Precision = 10^-4
e = 1
#show
BTDB_e(i; j; ξ; η) = transp(B(i; ξ; η))*D*B(j; ξ; η)
K_e(i; j) = l_1*h_1/4*$Area{$Area{BTDB_e(i; j; ξ; η) @ ξ = -1 : 1} @ η = -1 : 1}
$Repeat{$Repeat{K_e.(i; j) = K_e(i; j) @ j = i : n} @ i = 1 : n}
K_e'MN/m
'Boundary conditions
'Supports
'Number of elements along the supported edge -'n_S = len(e_S)
'Element′s joint springs stiffness factors
#noc
K_s,j(e) = l_1/2*$Area{N_j(ξ; -1)*s(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}', 'j = 1', '3
#equ
'Results for element 1
K_s,1(e) = l_1/2*$Area{N_1(ξ; -1)*s(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}','K_s,1(1)'MN/m
K_s,3(e) = l_1/2*$Area{N_3(ξ; -1)*s(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}','K_s,3(1)'MN/m
'Number of elements along the loaded edge -'n_L = len(e_L)
'Element load vector
#noc
F_j(e) = l_1/2*$Area{N_j(ξ; 1)*q(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}', 'j = 2', '4
#val
'Results for element 'e = e_L.(n_l/2)
#equ
F_2(e) = l_1/2*$Area{N_2(ξ; 1)*q(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}','F_2(e)'kN
F_4(e) = l_1/2*$Area{N_4(ξ; 1)*q(x_c.e + ξ*l_1/2) @ ξ = -1 : 1}','F_4(e)'kN
'<h3>Solution</h3>
'Global stiffness matrix
#hide
k_1 = n/4
n = k_1*n_j
K = symmetric(n)
'Addition of element stiffness matrix coefficients
#for e = 1 : n_e
    #for i = 1 : 4
        #for j = i : 4
            j_1 = e_j.(e; i)', 'i_1 = k_1*(j_1 - 1) + 1
            j_2 = e_j.(e; j)', 'i_2 = k_1*(j_2 - 1) + 1
            k_ij = submatrix(K_e; k_1*(i - 1) + 1; k_1*i; k_1*(j - 1) + 1; k_1*j)
            #if j_1 > j_2
                add(transp(k_ij); K; i_2; i_1)
            #else
                add(k_ij; K; i_1; i_2)
            #end if
        #loop
    #loop
#loop
'Addition of supports
#for i = 1 : n_S
    e = e_S.i
    K_se,1 = K_s,1(e)
    K_se,3 = K_s,3(e)
    j_1 = e_j.(e; 1)', 'i_1 = k_1*(j_1 - 1) + 2
    j_3 = e_j.(e; 3)', 'i_3 = k_1*(j_3 - 1) + 2
    K.(i_1; i_1) = K.(i_1; i_1) + K_se,1
    K.(i_3; i_3) = K.(i_3; i_3) + K_se,3
#loop
j_M = find(x_j; l/2; 1)
n_M = len(j_M)
#for i = 1 : n_M
    j = j_M.i', 'i_1 = k_1*(j - 1) + 1
    K.(i_1; i_1) = K.(i_1; i_1) + k_s
#loop
#show
K
'Global load vector
#hide
F = vector(n)
#for i = 1 : n_L
    e = e_L.i
    F_e,2 = F_2(e)
    F_e,4 = F_4(e)
    j_2 = e_j.(e; 2)', 'i_2 = k_1*(j_2 - 1) + 2
    j_4 = e_j.(e; 4)', 'i_4 = k_1*(j_4 - 1) + 2
    F.i_2 = F.i_2 + F_e,2
    F.i_4 = F.i_4 + F_e,4
#loop
#show
F
sum(F)'kN
'Solution of the system of equations
Z = clsolve(K; F)'mm
'<h3>Results</h3>
#hide
'Joint displacements
u = matrix(n_l + 1; n_h + 1)
u(j) = Z.(k_1*(j - 1) + 1)
$Repeat{$Repeat{u.(i; j) = round(u((i - 1)*(n_h + 1) + j)*1000000)/1000 @ j = 1 : n_h + 1} @ i = 1 : n_l + 1}
u(x; y) = spline(1 + x/l_1; 1 + y/h_1; u)
v = matrix(n_l + 1; n_h + 1)
v(j) = Z.(k_1*(j - 1) + 2)
$Repeat{$Repeat{v.(i; j) = round(v((i - 1)*(n_h + 1) + j)*1000000)/1000 @ j = 1 : n_h + 1} @ i = 1 : n_l + 1}
v(x; y) = spline(1 + x/l_1; 1 + y/h_1; v)
PlotWidth = 100*l
PlotHeight = 100*h
PlotStep = 10
#show
'Horizontal joint displacements, ·10<sup>-3</sup>mm
transp(u)
$Map{u(x; y) @ x = 0 : l & y = 0 : h}
'Vertical joint displacements, ·10<sup>-3</sup>mm
transp(v)
$Map{v(x; y) @ x = 0 : l & y = 0 : h}
'Calculation of internal forces
'Displacements for joint - 'Z_j(j) = slice(Z; k_1*(j - 1) + 1; k_1*j)
'Displacements for element - 'Z_e(e) = [Z_j(e_j.(e; 1)); Z_j(e_j.(e; 2)); Z_j(e_j.(e; 3)); Z_j(e_j.(e; 4))]
#hide
B(ξ; η) = [ _
B_1(1; η); B_1(2; η); B_1(3; η); B_1(4; η); _
B_1(5; η); B_1(6; η); B_1(7; η); B_1(8; η)| _
B_2(1; ξ); B_2(2; ξ); B_2(3; ξ); B_2(4; ξ); _
B_2(5; ξ); B_2(6; ξ); B_2(7; ξ); B_2(8; ξ)| _
B_3(1; ξ; η); B_3(2; ξ; η); B_3(3; ξ; η); B_3(4; ξ; η); _
B_3(5; ξ; η); B_3(6; ξ; η); B_3(7; ξ; η); B_3(8; ξ; η)]
#show
'Membrane forces in element -'N_e(e; x; y) = -D*B(2*x/l_1; 2*y/h_1)*Z_e(e)
#hide
N_j = matrix(3; n_j)
c_j = vector(n_j)
#for e = 1 : n_e
    Z_е = Z_e(e)
    j_1 = e_j.(e; 1)
    x_1 = x_c.e
    y_1 = y_c.e
    #for i = 1 : 4
        j = e_j.(e; i)
        c_j.j = c_j.j + 1
        x = x_j.j - x_1
        y = y_j.j - y_1
        N_e = N_e(e; x; y)
        add(N_e; N_j; 1; j)
    #loop
#loop
#for j = 1 : n_j
    N_j.(1; j) = N_j.(1; j)/c_j.j
    N_j.(2; j) = N_j.(2; j)/c_j.j
    N_j.(3; j) = N_j.(3; j)/c_j.j
#loop
Nx = matrix(n_l + 1; n_h + 1)
Ny = matrix(n_l + 1; n_h + 1)
Nxy = matrix(n_l + 1; n_h + 1)
#for i = 1 : n_l + 1
    #for k = 1 : n_h + 1
        j = (i - 1)*(n_h + 1) + k
        Nx.(i; k) = N_j.(1; j)
        Ny.(i; k) = N_j.(2; j)
        Nxy.(i; k) = N_j.(3; j)
    #loop
#loop
N_x(x; y) = spline(1 + x/l_1; 1 + y/h_1; Nx)
N_y(x; y) = spline(1 + x/l_1; 1 + y/h_1; Ny)
N_xy(x; y) = spline(1 + x/l_1; 1 + y/h_1; Nxy)
i = round(n_l/2)
j = round(n_h/2)
e = e_S.(n_l/2 + 1)
j = e_j.(e; 1)
x = -l_1/2
y = -h_1/2
#show
#val
'Results for element 'e' and joint 'j':
#equ
Z_e = Z_e(e)'mm
N_e = N_e(e; x; y)'kN/m
'Average membrane forces at joints, kN/m
N_j
'Membrane forces for the structure
'Normal membrane forces - <var>N</var><sub>x</sub>, kN/m
transp(Nx)
'<!--'k = round(abs(N_x(l/2; h))/20)*50'-->
$Map{N_x(x; y) @ x = 0 : l & y = 0 : h}
'Plot for <var>N</var><sub>x</sub>, kN/m at <var>x</var> = <var>l</var>/2
$Plot{N_x(l/2; y)|y & -k|0 & k|0 @ y = 0 : h}
'Bottom value -'N_x(l/2; 0)'kN/m
'Top value -'N_x(l/2; h)'kN/m
'Normal membrane forces - <var>N</var><sub>y</sub>, kN/m
transp(Ny)
$Map{N_y(x; y) @ x = 0 : l & y = 0 : h}
'<!--'k = round(abs(N_y(l/2; h))/20)*50'-->
$Plot{N_y(l/2; y)|y & -k|0 & k|0 @ y = 0 : h}
'Top value -'N_y(l/2; h)'kN/m
'Shear membrane forces - <var>N</var><sub>xy</sub>, kN/m
transp(Nxy)
$Map{N_xy(x; y) @ x = 0 : l & y = 0 : h}
'Max. value at 3/4 of span -'N_xy(3*l/4; h/2)'kN/m
'Principal membrane forces, kN/m
N_m(x; y) = 0.5*(N_x(x; y) + N_y(x; y))
ΔN(x; y) = 0.5*sqr((N_x(x; y) - N_y(x; y))^2 + 4*N_xy(x; y)^2)
N_max(x; y) = N_m(x; y) + ΔN(x; y)'kN/m
$Map{N_max(x; y) @ x = 0 : l & y = 0 : h}
N_min(x; y) = N_m(x; y) - ΔN(x; y)'kN/m
$Map{N_min(x; y) @ x = 0 : l & y = 0 : h}

Flat Slab FEA¶

Multi-span flat-slab analysis with the 16-DOF Bogner–Fox–Schmit plate element using numerical Gauss integration of the stiffness matrix. Handles arbitrary column grids and returns deflection together with the \(M_x\) and \(M_y\) bending-moment fields.

Code:

'Using numerical formulation of Bogner-Fox-Schmit (BFS)  plate element
'<h4>Input data</h4>
'Span lengths
a = hp([3.6; 4.2; 4.2; 3.6])'m
b = hp([3; 3.6; 3])'m
'Number of axes -'n_sa = len(a) + 1', 'n_sb = len(b) + 1
'<img src="./flat-slab.png" width="360">
#hide
x_s = vector_hp(n_sa)','y_s = vector_hp(n_sb)
$Repeat{x_s.(k + 1) = x_s.k + a.k @ k = 1 : n_sa - 1}
$Repeat{y_s.(k + 1) = y_s.k + b.k @ k = 1 : n_sb - 1}
#show
'Axis coordinates -'x_s'm, 'y_s'm
'Slab dimensions -'l_a = x_s.n_sa'm, 'l_b = y_s.n_sb'm
'Thickness -'t = 0.2'm
'Load -'q = 10'kN/m²
'Modulus of elasticity -'E = 35000'MPa
'Poisson`s ratio -'ν = 0.2
'<h4>Finite element mesh</h4>
'We will use Bogner-Fox-Schmit rectangular finite element with'n_DOFs = 16
'Element dimensions -'a_1 = 0.6'm, 'b_1 = 0.6'm
'Number of elements and joints along <var>a</var> and <var>b</var> -'
n_a = ceiling(a/a_1)', 'n_ea = sum(n_a)', 'n_ja = n_ea + 1
n_b = ceiling(b/b_1)', 'n_eb = sum(n_b)', 'n_jb = n_eb + 1
#if count(a ≡ a_1*n_a; 0; 1) + count(b ≡ b_1*n_b; 0; 1) > 0
    '<p class="err">Error! All span lengts must be multiples of element size.</p>
    #break
#end if
'Total number of elements -'n_e = n_ea*n_eb
'Total number of joints -'n_j = n_ja*n_jb
'Supported joints count -'n_s = n_sa*n_sb
#hide
x_j = vector_hp(n_j)','y_j = vector_hp(n_j)
x = 0', 'y = 0', 'j = 0
#for j_a = 1 : n_ja
    #for j_b = 1 : n_jb
        j = j + 1
        x_j.j = x
        y_j.j = y
        y = y + b_1
        #if y > l_b
            y = 0
            i_b = 1
            x = x + a_1
        #end if
    #loop
#loop
e_j = matrix_hp(n_e; 4)
#for e_a = 1 : n_ea
    #for e_b = 1 : n_eb
        e = e_b + n_eb*(e_a - 1)
        j = e + e_a - 1
        e_j.(e; 1) = j
        e_j.(e; 2) = j + n_eb + 1
        e_j.(e; 3) = j + n_eb + 2
        e_j.(e; 4) = j + 1
    #loop
#loop
s_j = vector(n_s)', 'j_s = 0', 'j = 1', 's_a = 0
#for i_a = 1 : n_sa
    s_b = 0
    #for i_b = 1 : n_sb
        j = j + s_b
        j_s = j_s + 1
        s_j.j_s = j
        #if i_b < n_sb
            s_b = n_b.i_b
        #end if
    #loop
    #if i_a < n_sa
        s_a = n_a.i_a
        j = j + (s_a - 1)*n_jb + 1
    #end if
#loop
#show
'Joint coordinates
x_j'm
y_j'm
'Numbers of joints at elements'' corners
transp(e_j)
'Supported joints
s_j
'Joints for element e -'j_e(e) = row(e_j; e)
#val
#hide
w = 400
k = w/l_a
d = 10
h = l_b/l_a*w
r = 0.04*k
#show
'<svg viewbox="'-d' '-d' 'w + 2*d' 'h + 2*d'" xmlns="http://www.w3.org/2000/svg" version="1.1" style=" font-family: Segoe UI; font-size:5px; width:'w'pt; height:'h'pt">
'<style>.joint{fill:orangeRed;} .support{stroke:red; stroke-width:1; fill:lightpink;} .element{stroke:seaGreen; stroke-width:1; fill:lime; fill-opacity:0.1; stroke-opacity:0.5}</style>
'<rect x="'0'" y="'0'" width="'w'" height="'h'" style="fill:yellow; fill-opacity:0.2" />
#for e = 1 : n_e
    #hide
    j_e = j_e(e)
    j1 = j_e.1
    x1 = x_j.j1*k
    y1 = y_j.j1*k
    j2 = j_e.3
    x2 = x_j.j2*k
    y2 = y_j.j2*k
    x = (x1 + x2)/2
    y = (y1 + y2)/2
    #show
    '<text x="'x'" y="'h - y'" text-anchor="middle">'e'</text>
    '<rect x="'x1'" y="'h - y2'" width="'x2 - x1'" height="'y2 - y1'" class="element" />
#loop
#for i = 1 : n_s
    j = s_j.i
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'2*r'" class="support"/>
#loop
#for j = 1 : n_j
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'r'" class="joint" />
    '<text x="'x + 2*r'" y="'y - r'" text-anchor="start">'j'</text>
#loop
'</svg>
#equ
'<h4>Finite element formulation</h4>
'<p><b>Shape functions</b></p>
#def FI$(x$; l$)
    '<table><tr><td>Base functions
    Φ_1l$(x$) = 1 - x$^2*(3 - 2*x$)
    Φ_2l$(x$) = x$*l$_1*(1 - x$*(2 - x$))
    Φ_3l$(x$) = x$^2*(3 - 2*x$)
    Φ_4l$(x$) = x$^2*l$_1*(-1 + x$)
    '</td><td>First derivatives
    Φ′_1l$(x$) = -6*(x$/l$_1)*(1 - x$)
    Φ′_2l$(x$) = 1 - x$*(4 - 3*x$)
    Φ′_3l$(x$) = 6*(x$/l$_1)*(1 - x$)
    Φ′_4l$(x$) = -x$*(2 - 3*x$)
    '</td><td>Second derivatives
    Φ″_1l$(x$) = -(6/l$_1^2)*(1 - 2*x$)
    Φ″_2l$(x$) = -(2/l$_1)*(2 - 3*x$)
    Φ″_3l$(x$) = 6/l$_1^2*(1 - 2*x$)
    Φ″_4l$(x$) = -(2/l$_1)*(1 - 3*x$)
    '</td></tr></table>
#end def
'Along dimension <var>a</var>
FI$(ξ; a)
'Along dimension <var>b</var>
FI$(η; b)
'<table><tr><td>For vertical displacements <var>w</var>
N_1,w(ξ; η) = Φ_1a(ξ)*Φ_1b(η)
N_2,w(ξ; η) = Φ_3a(ξ)*Φ_1b(η)
N_3,w(ξ; η) = Φ_3a(ξ)*Φ_3b(η)
N_4,w(ξ; η) = Φ_1a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ₓ
N_1,θₓ(ξ; η) = Φ_2a(ξ)*Φ_1b(η)
N_2,θₓ(ξ; η) = Φ_4a(ξ)*Φ_1b(η)
N_3,θₓ(ξ; η) = Φ_4a(ξ)*Φ_3b(η)
N_4,θₓ(ξ; η) = Φ_2a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ᵧ
N_1,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_2b(η)
N_2,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_2b(η)
N_3,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_4b(η)
N_4,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_4b(η)
'</td></tr><tr><td>
'For twist ψ
N_1,ψ(ξ; η) = Φ_2a(ξ)*Φ_2b(η)
N_2,ψ(ξ; η) = Φ_4a(ξ)*Φ_2b(η)
N_3,ψ(ξ; η) = Φ_4a(ξ)*Φ_4b(η)
N_4,ψ(ξ; η) = Φ_2a(ξ)*Φ_4b(η)
'</td><td colspan="2" style="width:160pt; text-align:center;">
'<img src="./plate-element.png" alt="plane-elemen.png" style="width:160pt;" />
'</td></tr></table>
'<!--'PlotSVG = 1''PlotWidth = 160''PlotHeight = PlotWidth''PlotPalette = 7'-->
'<table><tr><td>
'<p><var>N</var><sub>1,w</sub> shape function plot</p>
$Map{N_1,w(ξ; η) @ ξ = 0 : 1 & η = 0 : 1}
'</td><td>
'<p><var>N</var><sub>1,θₓ</sub> shape function plot</p>
$Map{N_1,θₓ(ξ; η) @ ξ = 0 : 1 & η = 0 : 1}
'</td></tr><tr><td>
'<!--'PlotHeight = 50'-->
$Plot{N_1,w(ξ; 0) @ ξ = 0 : 1}
'</td><td>
$Plot{N_1,θₓ(ξ; 0) @ ξ = 0 : 1}
'</td></tr></table>
'Shape functions vector
N(i; ξ; η) = take(i; _
N_1,w(ξ; η); N_1,θₓ(ξ; η); N_1,θᵧ(ξ; η); N_1,ψ(ξ; η); _
N_2,w(ξ; η); N_2,θₓ(ξ; η); N_2,θᵧ(ξ; η); N_2,ψ(ξ; η); _
N_3,w(ξ; η); N_3,θₓ(ξ; η); N_3,θᵧ(ξ; η); N_3,ψ(ξ; η); _
N_4,w(ξ; η); N_4,θₓ(ξ; η); N_4,θᵧ(ξ; η); N_4,ψ(ξ; η))
'<p><b>Constitutive matrix</b> (stress - strain relationship)</p>
D = E*t^3/(12*(1 - ν^2))*hp([1; ν; 0|ν; 1; 0|0; 0; (1 - ν)/2])'kNm
'<p><b>Strain-displacement matrix</b></p>
B_1(j; ξ; η) = take(j; _
Φ″_1a(ξ)*Φ_1b(η); Φ″_2a(ξ)*Φ_1b(η); Φ″_1a(ξ)*Φ_2b(η); Φ″_2a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_1b(η); Φ″_4a(ξ)*Φ_1b(η); Φ″_3a(ξ)*Φ_2b(η); Φ″_4a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_3b(η); Φ″_4a(ξ)*Φ_3b(η); Φ″_3a(ξ)*Φ_4b(η); Φ″_4a(ξ)*Φ_4b(η); _
Φ″_1a(ξ)*Φ_3b(η); Φ″_2a(ξ)*Φ_3b(η); Φ″_1a(ξ)*Φ_4b(η); Φ″_2a(ξ)*Φ_4b(η))
B_2(j; ξ; η) = take(j; _
Φ_1a(ξ)*Φ″_1b(η); Φ_2a(ξ)*Φ″_1b(η); Φ_1a(ξ)*Φ″_2b(η); Φ_2a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_1b(η); Φ_4a(ξ)*Φ″_1b(η); Φ_3a(ξ)*Φ″_2b(η); Φ_4a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_3b(η); Φ_4a(ξ)*Φ″_3b(η); Φ_3a(ξ)*Φ″_4b(η); Φ_4a(ξ)*Φ″_4b(η); _
Φ_1a(ξ)*Φ″_3b(η); Φ_2a(ξ)*Φ″_3b(η); Φ_1a(ξ)*Φ″_4b(η); Φ_2a(ξ)*Φ″_4b(η))
B_3(j; ξ; η) = 2*take(j; _
Φ′_1a(ξ)*Φ′_1b(η); Φ′_2a(ξ)*Φ′_1b(η); Φ′_1a(ξ)*Φ′_2b(η); Φ′_2a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_1b(η); Φ′_4a(ξ)*Φ′_1b(η); Φ′_3a(ξ)*Φ′_2b(η); Φ′_4a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_3b(η); Φ′_4a(ξ)*Φ′_3b(η); Φ′_3a(ξ)*Φ′_4b(η); Φ′_4a(ξ)*Φ′_4b(η); _
Φ′_1a(ξ)*Φ′_3b(η); Φ′_2a(ξ)*Φ′_3b(η); Φ′_1a(ξ)*Φ′_4b(η); Φ′_2a(ξ)*Φ′_4b(η))
B(j; ξ; η) = hp([B_1(j; ξ; η); B_2(j; ξ; η); B_3(j; ξ; η)])
'The coefficients of the stiffness matrix will be calculated by using the equation
#noc
K_e,ij = a_1*b_1*$Integral{$Integral{B_i(ξ; η)^T*D*B_j(ξ; η) @ ξ = 0 : 1} @ η = 0 : 1}
#equ
'<p><b>Element stiffness matrix</b></p> (above the main diagonal only)
#hide
n = n_DOFs
Precision = 10^-4
K_e = utriang_hp(n)
#show
BTDB_e(i; j; ξ; η) = transp(B(i; ξ; η))*D*B(j; ξ; η)
K_e(i; j) = a_1*b_1*$Integral{$Integral{BTDB_e(i; j; ξ; η) @ η = 0 : 1} @ ξ = 0 : 1}
$Repeat{$Repeat{K_e.(i; j) = K_e(i; j) @ j = i : n} @ i = 1 : n}
K_e
'Element load vector
#noc
F_e,i = a_1*b_1*$Integral{$Integral{N_i(ξ; η)^T*q @ ξ = 0 : 1} @ η = 0 : 1}
#equ
#hide
F_e = vector_hp(n)
#for i = 1 : n
    F_e.i = a_1*b_1*$Integral{$Integral{N(i; ξ; η)*q @ ξ = 0 : 1} @ η = 0 : 1}
#loop
#show
F_e'kN
'<h4>Solution</h4>
'Global stiffness matrix
#hide
k_1 = n/4
n = k_1*n_j
K = symmetric_hp(n)
'Addition of element stiffness matrix coefficients
#for e = 1 : n_e
    #for i = 1 : 4
        #for j = i : 4
            j_1 = e_j.(e; i)
            j_2 = e_j.(e; j)
            i_1 = k_1*(j_1 - 1) + 1
            i_2 = k_1*(j_2 - 1) + 1
            k_ij = submatrix(K_e; k_1*(i - 1) + 1; k_1*i; k_1*(j - 1) + 1; k_1*j)
            #if j_1 > j_2
                add(transp(k_ij); K; i_2; i_1)
            #else
                add(k_ij; K; i_1; i_2)
            #end if
        #loop
    #loop
#loop
k_s = 10^20
'Addition of supports
#for i = 1 : n_s
    j = k_1*(s_j.i - 1) + 1
    K.(j; j) = K.(j; j) + k_s
#loop
#show
K
'Global load vector
#hide
F = vector_hp(n)
#for e = 1 : n_e
    #for i = 1 : 4
        j = e_j.(e; i)
        #for k = 1 : 4
            l = k_1*(j - 1) + k
            F.l = F.l + F_e.(k_1*(i - 1) + k)
        #loop
    #loop
#loop
#show
F'kN
'Solution of the system of equations<!--'Tol = 10^-2'-->
Z = slsolve(K; F)'mm
'<h4>Results</h4>
'Joint displacements
#hide
W_z = matrix_hp(n_ja; n_jb)
w(j) = Z.(4*j - 3)
$Repeat{$Repeat{W_z.(i; j) = round(w((i - 1)*n_jb + j)*1000)/1000 @ j = 1 : n_jb} @ i = 1 : n_ja}
w(x; y) = spline(1 + x/a_1; 1 + y/b_1; W_z)
PlotWidth = 450', 'PlotHeight = PlotWidth*l_b/l_a
PlotPalette = 2', 'PlotShadows = 1
PlotLightDir = 7', 'PlotStep = 2
#show
transp(W_z)'mm
$Map{-w(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments
Z_j(j) = slice(Z; k_1*(j - 1) + 1; k_1*j)
Z_e(e) = hp([Z_j(e_j.(e; 1)); Z_j(e_j.(e; 2)); Z_j(e_j.(e; 3)); Z_j(e_j.(e; 4))])
#hide
B(ξ; η) = hp([ _
B_1(1; ξ; η); B_1(2; ξ; η); B_1(3; ξ; η); B_1(4; ξ; η); _
B_1(5; ξ; η); B_1(6; ξ; η); B_1(7; ξ; η); B_1(8; ξ; η); _
B_1(9; ξ; η); B_1(10; ξ; η); B_1(11; ξ; η); B_1(12; ξ; η); _
B_1(13; ξ; η); B_1(14; ξ; η); B_1(15; ξ; η); B_1(16; ξ; η)| _
B_2(1; ξ; η); B_2(2; ξ; η); B_2(3; ξ; η); B_2(4; ξ; η); _
B_2(5; ξ; η); B_2(6; ξ; η); B_2(7; ξ; η); B_2(8; ξ; η); _
B_2(9; ξ; η); B_2(10; ξ; η); B_2(11; ξ; η); B_2(12; ξ; η); _
B_2(13; ξ; η); B_2(14; ξ; η); B_2(15; ξ; η); B_2(16; ξ; η)| _
B_3(1; ξ; η); B_3(2; ξ; η); B_3(3; ξ; η); B_3(4; ξ; η); _
B_3(5; ξ; η); B_3(6; ξ; η); B_3(7; ξ; η); B_3(8; ξ; η); _
B_3(9; ξ; η); B_3(10; ξ; η); B_3(11; ξ; η); B_3(12; ξ; η); _
B_3(13; ξ; η); B_3(14; ξ; η); B_3(15; ξ; η); B_3(16; ξ; η)])
M_j = matrix_hp(3; n_j)
c_j = vector_hp(n_j)
#for e = 1 : n_e
    Z_е = Z_e(e)
    j_1 = e_j.(e; 1)
    x_1 = x_j.j_1
    y_1 = y_j.j_1
    M_e(x; y) = -D*B(x/a_1; y/b_1)*Z_е
    #for i = 1 : 4
        j = e_j.(e; i)
        c_j.j = c_j.j + 1
        x = x_j.j - x_1
        y = y_j.j - y_1
        M_e = M_e(x; y)
        add(M_e; M_j; 1; j)
    #loop
#loop
#for j = 1 : n_j
     M_j.(1; j) = M_j.(1; j)/c_j.j
     M_j.(2; j) = M_j.(2; j)/c_j.j
     M_j.(3; j) = M_j.(3; j)/c_j.j
#loop
Mx = matrix_hp(n_ja; n_jb)
My = matrix_hp(n_ja; n_jb)
Mxy = matrix_hp(n_ja; n_jb)
#for i = 1 : n_ja
    #for k = 1 : n_jb
        j = (i - 1)*n_jb + k
        Mx.(i; k) = M_j.(1; j)
        My.(i; k) = M_j.(2; j)
        Mxy.(i; k) = M_j.(3; j)
    #loop
#loop
M_x(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mx)
M_y(x; y) = spline(1 + x/a_1; 1 + y/b_1; My)
M_xy(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mxy)
PlotShadows = 0
#show
'Average bending moments at joints, kNm/m
M_j
'Bending moments for the plate
'Bending moments - <var>M</var><sub>x</sub>
transp(Mx)'kNm/m
$Map{M_x(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments <var>M</var><sub>y</sub>
transp(My)'kNm/m
$Map{M_y(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments <var>M</var><sub>xy</sub>
transp(Mxy)'kNm/m
$Map{M_xy(x; y) @ x = 0 : l_a & y = 0 : l_b}

Flat Slab FEA Optimized¶

Same Bogner–Fox–Schmit flat-slab model, but with closed-form analytically integrated element matrices instead of Gauss quadrature. An order of magnitude faster on dense meshes — well suited for parametric studies and mesh-refinement sweeps.

Code:

'Using analytical formulation of Bogner-Fox-Schmit (BFS)  plate element
'<h4>Input data</h4>
'Span lengths
a = hp([3.6; 4.2; 4.2; 3.6])'m
b = hp([3; 3.6; 3])'m
'Number of axes -'n_sa = len(a) + 1', 'n_sb = len(b) + 1
'<img src="./flat-slab.png" width="360">
#hide
x_s = vector_hp(n_sa)','y_s = vector_hp(n_sb)
$Repeat{x_s.(k + 1) = x_s.k + a.k @ k = 1 : n_sa - 1}
$Repeat{y_s.(k + 1) = y_s.k + b.k @ k = 1 : n_sb - 1}
#show
'Axis coordinates -'x_s'm, 'y_s'm
'Slab dimensions -'l_a = x_s.n_sa'm, 'l_b = y_s.n_sb'm
'Thickness -'t = 0.2'm
'Load -'q = 10'kN/m²
'Modulus of elasticity -'E = 35000'MPa
'Poisson`s ratio -'ν = 0.2
'<h4>Finite element mesh</h4>
'We will use BFS rectangular finite element with'n_DOFs = 16
'Element dimensions -'a_1 = 0.6'm, 'b_1 = 0.6'm
'Number of elements and joints along <var>a</var> and <var>b</var> -'
n_a = ceiling(a/a_1)', 'n_ea = sum(n_a)', 'n_ja = n_ea + 1
n_b = ceiling(b/b_1)', 'n_eb = sum(n_b)', 'n_jb = n_eb + 1
#if count(a ≡ a_1*n_a; 0; 1) + count(b ≡ b_1*n_b; 0; 1) > 0
    '<p class="err">Error! All span lengts must be multiples of element size.</p>
    #break
#end if
'Total number of elements -'n_e = n_ea*n_eb
'Total number of joints -'n_j = n_ja*n_jb
'Supported joints count -'n_s = n_sa*n_sb
#hide
x_j = vector_hp(n_j)','y_j = vector_hp(n_j)
x = 0', 'y = 0', 'j = 0
#for j_a = 1 : n_ja
    #for j_b = 1 : n_jb
        j = j + 1
        x_j.j = x
        y_j.j = y
        y = y + b_1
        #if y > l_b
            y = 0
            i_b = 1
            x = x + a_1
        #end if
    #loop
#loop
e_j = matrix_hp(n_e; 4)
#for e_a = 1 : n_ea
    #for e_b = 1 : n_eb
        e = e_b + n_eb*(e_a - 1)
        j = e + e_a - 1
        e_j.(e; 1) = j
        e_j.(e; 2) = j + n_eb + 1
        e_j.(e; 3) = j + n_eb + 2
        e_j.(e; 4) = j + 1
    #loop
#loop
s_j = vector(n_s)', 'j_s = 0', 'j = 1', 's_a = 0
#for i_a = 1 : n_sa
    s_b = 0
    #for i_b = 1 : n_sb
        j = j + s_b
        j_s = j_s + 1
        s_j.j_s = j
        #if i_b < n_sb
            s_b = n_b.i_b
        #end if
    #loop
    #if i_a < n_sa
        s_a = n_a.i_a
        j = j + (s_a - 1)*n_jb + 1
    #end if
#loop
#show
'Joint coordinates
x_j'm
y_j'm
'Numbers of joints at elements'' corners
transp(e_j)
'Supported joints
s_j
'Joints for element e -'j_e(e) = row(e_j; e)
#val
#hide
w = 400
k = w/l_a
d = 10
h = l_b/l_a*w
r = 0.04*k
#show
'<svg viewbox="'-d' '-d' 'w + 2*d' 'h + 2*d'" xmlns="http://www.w3.org/2000/svg" version="1.1" style=" font-family: Segoe UI; font-size:5px; width:'w'pt; height:'h'pt">
'<style>.joint{fill:orangeRed;} .support{stroke:red; stroke-width:1; fill:lightpink;} .element{stroke:seaGreen; stroke-width:1; fill:lime; fill-opacity:0.1; stroke-opacity:0.5}</style>
'<rect x="'0'" y="'0'" width="'w'" height="'h'" style="fill:yellow; fill-opacity:0.2" />
#for e = 1 : n_e
    #hide
    j_e = j_e(e)
    j1 = j_e.1
    x1 = x_j.j1*k
    y1 = y_j.j1*k
    j2 = j_e.3
    x2 = x_j.j2*k
    y2 = y_j.j2*k
    x = (x1 + x2)/2
    y = (y1 + y2)/2
    #show
    '<text x="'x'" y="'h - y'" text-anchor="middle">'e'</text>
    '<rect x="'x1'" y="'h - y2'" width="'x2 - x1'" height="'y2 - y1'" class="element" />
#loop
#for i = 1 : n_s
    j = s_j.i
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'2*r'" class="support"/>
#loop
#for j = 1 : n_j
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'r'" class="joint" />
    '<text x="'x + 2*r'" y="'y - r'" text-anchor="start">'j'</text>
#loop
'</svg>
#equ
'<h4>Finite element formulation</h4>
'<p><b>Shape functions</b></p>
#def FI$(x$; l$)
    '<table><tr><td>Base functions
    Φ_1l$(x$) = 1 - x$^2*(3 - 2*x$)
    Φ_2l$(x$) = x$*l$_1*(1 - x$*(2 - x$))
    Φ_3l$(x$) = x$^2*(3 - 2*x$)
    Φ_4l$(x$) = x$^2*l$_1*(-1 + x$)
    '</td><td>First derivatives
    Φ′_1l$(x$) = -6*(x$/l$_1)*(1 - x$)
    Φ′_2l$(x$) = 1 - x$*(4 - 3*x$)
    Φ′_3l$(x$) = 6*(x$/l$_1)*(1 - x$)
    Φ′_4l$(x$) = -x$*(2 - 3*x$)
    '</td><td>Second derivatives
    Φ″_1l$(x$) = -(6/l$_1^2)*(1 - 2*x$)
    Φ″_2l$(x$) = -(2/l$_1)*(2 - 3*x$)
    Φ″_3l$(x$) = 6/l$_1^2*(1 - 2*x$)
    Φ″_4l$(x$) = -(2/l$_1)*(1 - 3*x$)
    '</td></tr></table>
#end def
'Along dimension <var>a</var>
FI$(ξ; a)
'Along dimension <var>b</var>
FI$(η; b)
'<table><tr><td>For vertical displacements <var>w</var>
N_1,w(ξ; η) = Φ_1a(ξ)*Φ_1b(η)
N_2,w(ξ; η) = Φ_3a(ξ)*Φ_1b(η)
N_3,w(ξ; η) = Φ_3a(ξ)*Φ_3b(η)
N_4,w(ξ; η) = Φ_1a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ₓ
N_1,θₓ(ξ; η) = Φ_2a(ξ)*Φ_1b(η)
N_2,θₓ(ξ; η) = Φ_4a(ξ)*Φ_1b(η)
N_3,θₓ(ξ; η) = Φ_4a(ξ)*Φ_3b(η)
N_4,θₓ(ξ; η) = Φ_2a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ᵧ
N_1,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_2b(η)
N_2,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_2b(η)
N_3,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_4b(η)
N_4,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_4b(η)
'</td></tr><tr><td>
'For twist ψ
N_1,ψ(ξ; η) = Φ_2a(ξ)*Φ_2b(η)
N_2,ψ(ξ; η) = Φ_4a(ξ)*Φ_2b(η)
N_3,ψ(ξ; η) = Φ_4a(ξ)*Φ_4b(η)
N_4,ψ(ξ; η) = Φ_2a(ξ)*Φ_4b(η)
'</td><td colspan="2" style="width:160pt; text-align:center;">
'<img src="./plate-element.png" alt="plane-elemen.png" style="width:160pt;" />
'</td></tr></table>
'<!--'PlotSVG = 1''PlotWidth = 160''PlotHeight = PlotWidth''PlotPalette = 7'-->
'<table><tr><td>
'<p><var>N</var><sub>1,w</sub> shape function plot</p>
$Map{N_1,w(ξ; η) @ ξ = 0 : 1 & η = 0 : 1}
'</td><td>
'<p><var>N</var><sub>1,θₓ</sub> shape function plot</p>
$Map{N_1,θₓ(ξ; η) @ ξ = 0 : 1 & η = 0 : 1}
'</td></tr><tr><td>
'<!--'PlotHeight = 40'-->
$Plot{N_1,w(ξ; 0) @ ξ = 0 : 1}
'</td><td>
$Plot{N_1,θₓ(ξ; 0) @ ξ = 0 : 1}
'</td></tr></table>
'<p><b>Constitutive matrix</b> (stress - strain relationship)</p>
D_1,1 = E*t^3/(12*(1 - ν^2))'kNm
D = D_1,1*hp([1; ν; 0|ν; 1; 0|0; 0; (1 - ν)/2])'kNm
#hide
'<p><b>Strain-displacement matrix</b></p>
B_1(j; ξ; η) = take(j; _
Φ″_1a(ξ)*Φ_1b(η); Φ″_2a(ξ)*Φ_1b(η); Φ″_1a(ξ)*Φ_2b(η); Φ″_2a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_1b(η); Φ″_4a(ξ)*Φ_1b(η); Φ″_3a(ξ)*Φ_2b(η); Φ″_4a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_3b(η); Φ″_4a(ξ)*Φ_3b(η); Φ″_3a(ξ)*Φ_4b(η); Φ″_4a(ξ)*Φ_4b(η); _
Φ″_1a(ξ)*Φ_3b(η); Φ″_2a(ξ)*Φ_3b(η); Φ″_1a(ξ)*Φ_4b(η); Φ″_2a(ξ)*Φ_4b(η))
B_2(j; ξ; η) = take(j; _
Φ_1a(ξ)*Φ″_1b(η); Φ_2a(ξ)*Φ″_1b(η); Φ_1a(ξ)*Φ″_2b(η); Φ_2a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_1b(η); Φ_4a(ξ)*Φ″_1b(η); Φ_3a(ξ)*Φ″_2b(η); Φ_4a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_3b(η); Φ_4a(ξ)*Φ″_3b(η); Φ_3a(ξ)*Φ″_4b(η); Φ_4a(ξ)*Φ″_4b(η); _
Φ_1a(ξ)*Φ″_3b(η); Φ_2a(ξ)*Φ″_3b(η); Φ_1a(ξ)*Φ″_4b(η); Φ_2a(ξ)*Φ″_4b(η))
B_3(j; ξ; η) = 2*take(j; _
Φ′_1a(ξ)*Φ′_1b(η); Φ′_2a(ξ)*Φ′_1b(η); Φ′_1a(ξ)*Φ′_2b(η); Φ′_2a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_1b(η); Φ′_4a(ξ)*Φ′_1b(η); Φ′_3a(ξ)*Φ′_2b(η); Φ′_4a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_3b(η); Φ′_4a(ξ)*Φ′_3b(η); Φ′_3a(ξ)*Φ′_4b(η); Φ′_4a(ξ)*Φ′_4b(η); _
Φ′_1a(ξ)*Φ′_3b(η); Φ′_2a(ξ)*Φ′_3b(η); Φ′_1a(ξ)*Φ′_4b(η); Φ′_2a(ξ)*Φ′_4b(η))
n = n_DOFs
K_e = utriang_hp(n)
#show
'<div class="fold">
'<h5>Element stiffness matrix calculation</h5>
α = b_1/a_1
#hide
α² = α^2
A_1 = a_1*b_1
#show
K_e.(1; 1) = 156/35*(α/a_1^2 + 1/(α*b_1^2)) + 72/(25*A_1)
K_e.(5; 5) = K_e.(1; 1)', ' _
K_e.(9; 9) = K_e.(1; 1)', ' _
K_e.(13; 13) = K_e.(1; 1)
K_e.(1; 2) = 2/35*(39*α/a_1 + 11/(α²*b_1)) + (30*ν + 6)/(25*b_1)
K_e.(5; 6) = -K_e.(1; 2)', ' _
K_e.(9; 10) = -K_e.(1; 2)', ' _
K_e.(13; 14) = K_e.(1; 2)
K_e.(1; 3) = 2/35*(11*α²/a_1 + 39/(α*b_1)) + (30*ν + 6)/(25*a_1)
K_e.(5; 7) = K_e.(1; 3)', ' _
K_e.(9; 11) = -K_e.(1; 3)', ' _
K_e.(13; 15) = -K_e.(1; 3)
K_e.(1; 4) = 11/35*(α² + 1/α²) + (10*ν + 1)/50
K_e.(5; 8) = -K_e.(1; 4)', ' _
K_e.(9; 12) = K_e.(1; 4)', ' _
K_e.(13; 16) = -K_e.(1; 4)
K_e.(1; 5) = 2/35*(27/(α*b_1^2) - 78*α/a_1^2) - 72/(25*A_1)
K_e.(9; 13) = K_e.(1; 5)
K_e.(1; 6) = 13/35*(6*α/a_1 - 1/(α²*b_1)) + 6/(25*b_1)
K_e.(2; 5) = -K_e.(1; 6)', ' _
K_e.(9; 14) = -K_e.(1; 6)', ' _
K_e.(10; 13) = K_e.(1; 6)
K_e.(1; 7) = 1/35*(27/(α*b_1) - 22*α²/a_1) - (30*ν + 6)/(25*a_1)
K_e.(3; 5) = K_e.(1; 7)', ' _
K_e.(9; 15) = -K_e.(1; 7)', ' _
K_e.(11; 13) = -K_e.(1; 7)
K_e.(2; 4) = 2/35*(11*α*b_1/3 + a_1/α²) + 2*a_1*(5*ν + 1)/75
K_e.(6; 8) = K_e.(2; 4)', ' _
K_e.(10; 12) = -K_e.(2; 4)', ' _
K_e.(14; 16) = -K_e.(2; 4)
K_e.(2; 6) = 1/35*(26*α - 3/α^3) - 2/(25*α)
K_e.(10; 14) = K_e.(2; 6)
K_e.(2; 8) = 1/35*(11*α*b_1/3 - 3*a_1/(2*α²)) - a_1*(5*ν + 1)/150
K_e.(4; 6) = K_e.(2; 8)', ' _
K_e.(10; 16) = -K_e.(2; 8)', ' _
K_e.(12; 14) = -K_e.(2; 8)
K_e.(2; 10) = 3/35*(1/α^3 + 3*α) + 2/(25*α)
K_e.(6; 14) = K_e.(2; 10)
K_e.(2; 12) = -(1/70)*(3*a_1/α² + 13*α*b_1/3) - a_1/150
K_e.(4; 10) = -K_e.(2; 12)', ' _
K_e.(6; 16) = K_e.(2; 12)', ' _
K_e.(8; 14) = -K_e.(2; 12)
K_e.(2; 14) = 2/35*(9*α - 2/α^3) - 8/(25*α)
K_e.(6; 10) = K_e.(2; 14)
K_e.(2; 15) = 1/35*(11/α² - 13*α²/2) + (5*ν + 1)/50
K_e.(3; 14) = -K_e.(2; 15)', ' _
K_e.(6; 11) = -K_e.(2; 15)', ' _
K_e.(7; 10) = K_e.(2; 15)
K_e.(1; 14) = 1/35*(27*α/a_1 - 22/(α²*b_1)) - (30*ν + 6)/(25*b_1)
K_e.(2; 13) = K_e.(1; 14)', ' _
K_e.(5; 10) = -K_e.(1; 14)', ' _
K_e.(6; 9) = -K_e.(1; 14)
K_e.(1; 15) = 13/35*(6/(α*b_1) - α²/a_1) + 6/(25*a_1)
K_e.(3; 13) = -K_e.(1; 15)', ' _
K_e.(5; 11) = K_e.(1; 15)', ' _
K_e.(7; 9) = -K_e.(1; 15)
K_e.(1; 16) = 1/35*(11/α² - 13*α²/2) + (5*ν + 1)/50
K_e.(4; 13) = -K_e.(1; 16)', ' _
K_e.(8; 9) = K_e.(1; 16)', ' _
K_e.(5; 12) = -K_e.(1; 16)
K_e.(2; 2) = 4/35*(13*α + 1/α^3) + 8/(25*α)
K_e.(6; 6) = K_e.(2; 2)', ' _
K_e.(10; 10) = K_e.(2; 2)', ' _
K_e.(14; 14) = K_e.(2; 2)
K_e.(2; 3) = 11/35*(α² + 1/α²) + (60*ν + 1)/50
K_e.(6; 7) = -K_e.(2; 3)', ' _
K_e.(10; 11) = K_e.(2; 3)', ' _
K_e.(14; 15) = -K_e.(2; 3)
K_e.(3; 16) = -(1/35)*(3*α²*b_1/2 - 11*a_1/(3*α)) - b_1*(5*ν + 1)/150
K_e.(4; 15) = K_e.(3; 16)', ' _
K_e.(7; 12) = -K_e.(3; 16)', ' _
K_e.(8; 11) = -K_e.(3; 16)
K_e.(4; 4) = 4/105*(α*b_1^2 + a_1^2/α) + 8*A_1/225
K_e.(8; 8) = K_e.(4; 4)', ' _
K_e.(12; 12) = K_e.(4; 4)', ' _
K_e.(16; 16) = K_e.(4; 4)
K_e.(4; 8) = 1/35*(2*α*b_1^2/3 - a_1^2/α) - 2*A_1/225
K_e.(12; 16) = K_e.(4; 8)
K_e.(4; 12) = -(1/70)*(a_1^2/α + α*b_1^2) + A_1/450
K_e.(8; 16) = K_e.(4; 12)
K_e.(4; 16) = 1/35*(2*a_1^2/(3*α) - α*b_1^2) - 2*A_1/225
K_e.(8; 12) = K_e.(4; 16)
K_e.(1; 8) = (1/35)*(11*α² - 13/(2*α²)) + (5*ν + 1)/50
K_e.(2; 7) = -K_e.(1; 8)', ' _
K_e.(3; 6) = K_e.(1; 8)', ' _
K_e.(4; 5) = -K_e.(1; 8)', ' _
K_e.(9; 16) = K_e.(1; 8)
K_e.(10; 15) = -K_e.(1; 8)', ' _
K_e.(11; 14) = K_e.(1; 8)', ' _
K_e.(12; 13) = -K_e.(1; 8)
K_e.(1; 9) = -(54/35)*(α/a_1^2 + 1/(α*b_1^2)) + 72/(25*A_1)
K_e.(5; 13) = K_e.(1; 9)
K_e.(1; 10) = 1/35*(27*α/a_1 + 13/(α²*b_1)) - 6/(25*b_1)
K_e.(2; 9) = -K_e.(1; 10)', ' _
K_e.(5; 14) = -K_e.(1; 10)', ' _
K_e.(6; 13) = K_e.(1; 10)
K_e.(1; 11) = 1/35*(13*α²/a_1 + 27/(α*b_1)) - 6/(25*a_1)
K_e.(3; 9) = -K_e.(1; 11)', ' _
K_e.(5; 15) = K_e.(1; 11)', ' _
K_e.(7; 13) = -K_e.(1; 11)
K_e.(1; 12) = -(13/70)*(α² + 1/α²) + 1/50
K_e.(2; 11) = -K_e.(1; 12)', ' _
K_e.(3; 10) = -K_e.(1; 12)', ' _
K_e.(4; 9) = K_e.(1; 12)', ' _
K_e.(5; 16) = -K_e.(1; 12)
K_e.(6; 15) = K_e.(1; 12)', ' _
K_e.(7; 14) = K_e.(1; 12)', ' _
K_e.(8; 13) = -K_e.(1; 12)
K_e.(1; 13) = 2/35*(27*α/a_1^2 - 78/(α*b_1^2)) - 72/(25*A_1)
K_e.(5; 9) = K_e.(1; 13)
K_e.(2; 16) = 1/35*(2*a_1/α² - 13*α*b_1/3) + 2*a_1/75
K_e.(4; 14) = -K_e.(2; 16)', ' _
K_e.(6; 12) = K_e.(2; 16)', ' _
K_e.(8; 10) = -K_e.(2; 16)
K_e.(3; 3) = 4/35*(α^3 + 13/α) + 8*α/25
K_e.(7; 7) = K_e.(3; 3)', ' _
K_e.(11; 11) = K_e.(3; 3)', ' _
K_e.(15; 15) = K_e.(3; 3)
K_e.(3; 4) = 2/35*(α²*b_1 + 11*a_1/(3*α)) + 2*b_1*(5*ν + 1)/75
K_e.(7; 8) = -K_e.(3; 4)', ' _
K_e.(11; 12) = -K_e.(3; 4)', ' _
K_e.(15; 16) = K_e.(3; 4)
K_e.(3; 7) = 2/35*(-2*α^3 + 9/α) - 8*α/25
K_e.(11; 15) = K_e.(3; 7)
K_e.(3; 8) = 1/35*(2*α²*b_1 - 13*a_1/(3*α)) + 2*b_1/75
K_e.(4; 7) = -K_e.(3; 8)', ' _
K_e.(11; 16) = -K_e.(3; 8)', ' _
K_e.(12; 15) = K_e.(3; 8)
K_e.(3; 11) = 3/35*(α^3 + 3/α) + 2*α/25
K_e.(7; 15) = K_e.(3; 11)
K_e.(3; 12) = -(1/70)*(3*α²*b_1 + 13*a_1/(3*α)) - b_1/150
K_e.(4; 11) = -K_e.(3; 12)', ' _
K_e.(7; 16) = -K_e.(3; 12)', ' _
K_e.(8; 15) = K_e.(3; 12)
K_e.(3; 15) = 1/35*(26/α - 3*α^3) - 2*α/25
K_e.(7; 11) = K_e.(3; 15)
'</div>
'Element stiffness matrix coefficients (above the main diagonal only)
K_e = D_1,1*K_e
'<p><strong>Element load vector</strong></p>
F_e = q*A_1/24*hp([6; a_1; b_1; A_1/6; 6; -a_1; b_1; -A_1/6; 6; -a_1; -b_1; A_1/6; 6; a_1; -b_1; -A_1/6])'kN
'<h4>Solution</h4>
'Global stiffness matrix
#hide
k_1 = n/4
k_1,1 = submatrix(K_e; 1; 4; 1; 4)
k_1,2 = submatrix(K_e; 1; 4; 5; 8)
k_1,3 = submatrix(K_e; 1; 4; 9; 12)
k_1,4 = submatrix(K_e; 1; 4; 13; 16)
k_2,2 = submatrix(K_e; 5; 8; 5; 8)
k_2,3 = submatrix(K_e; 5; 8; 9; 12)
k_2,4 = submatrix(K_e; 5; 8; 13; 16)
k_3,3 = submatrix(K_e; 9; 12; 9; 12)
k_3,4 = submatrix(K_e; 9; 12; 13; 16)
k_4,4 = submatrix(K_e; 13; 16; 13; 16)
n = k_1*n_j
K = symmetric_hp(n)
'Addition of element stiffness matrix coefficients
#for e = 1 : n_e
    #for i = 1 : 4
        #for j = i : 4
            j_1 = e_j.(e; i)
            j_2 = e_j.(e; j)
            i_1 = k_1*(j_1 - 1) + 1
            i_2 = k_1*(j_2 - 1) + 1
            #if i ≡ 1
                #if j ≡ 1
                    k_ij = k_1,1
                #else if j ≡ 2
                    k_ij = k_1,2
                #else if j ≡ 3
                    k_ij = k_1,3
                #else
                    k_ij = k_1,4
                #end if
            #else if i ≡ 2
                #if j ≡ 2
                    k_ij = k_2,2
                #else if j ≡ 3
                    k_ij = k_2,3
                #else
                    k_ij = k_2,4
                #end if
            #else if i ≡ 3
                #if j ≡ 3
                    k_ij = k_3,3
                #else
                    k_ij = k_3,4
                #end if
            #else
                k_ij = k_4,4
            #end if
            #if j_1 > j_2
                add(transp(k_ij); K; i_2; i_1)
            #else
                add(k_ij; K; i_1; i_2)
            #end if
        #loop
    #loop
#loop
k_s = 10^20
'Addition of supports
#for i = 1 : n_s
    j = k_1*(s_j.i - 1) + 1
    K.(j; j) = K.(j; j) + k_s
#loop
#show
K
'Global load vector
#hide
F = vector_hp(n)
#for e = 1 : n_e
    #for i = 1 : 4
        j = e_j.(e; i)
        #for k = 1 : 4
            l = k_1*(j - 1) + k
            F.l = F.l + F_e.(k_1*(i - 1) + k)
        #loop
    #loop
#loop
#show
F'kN
'Solution of the system of equations<!--'Tol = 10^-2'-->
Z = slsolve(K; F)'mm
'<h4>Results</h4>
'Joint displacements
#hide
W_z = matrix_hp(n_ja; n_jb)
w(j) = Z.(4*j - 3)
$Repeat{$Repeat{W_z.(i; j) = round(w((i - 1)*n_jb + j)*1000)/1000 @ j = 1 : n_jb} @ i = 1 : n_ja}
w(x; y) = spline(1 + x/a_1; 1 + y/b_1; W_z)
PlotWidth = 450', 'PlotHeight = PlotWidth*l_b/l_a
PlotPalette = 2', 'PlotShadows = 1
PlotLightDir = 7', 'PlotStep = 5
#show
transp(W_z)'mm
$Map{-w(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments
Z_j(j) = slice(Z; k_1*(j - 1) + 1; k_1*j)
Z_e(e) = hp([Z_j(e_j.(e; 1)); Z_j(e_j.(e; 2)); Z_j(e_j.(e; 3)); Z_j(e_j.(e; 4))])
#hide
B_0,0 = matrix_hp(3; n_DOFs)
B_1,0 = matrix_hp(3; n_DOFs)
B_1,1 = matrix_hp(3; n_DOFs)
B_0,1 = matrix_hp(3; n_DOFs)
#for j = 1 : n_DOFs
    B_0,0.(1; j) = B_1(j; 0; 0)
    B_0,0.(2; j) = B_2(j; 0; 0)
    B_0,0.(3; j) = B_3(j; 0; 0)
    B_1,0.(1; j) = B_1(j; 1; 0)
    B_1,0.(2; j) = B_2(j; 1; 0)
    B_1,0.(3; j) = B_3(j; 1; 0)
    B_1,1.(1; j) = B_1(j; 1; 1)
    B_1,1.(2; j) = B_2(j; 1; 1)
    B_1,1.(3; j) = B_3(j; 1; 1)
    B_0,1.(1; j) = B_1(j; 0; 1)
    B_0,1.(2; j) = B_2(j; 0; 1)
    B_0,1.(3; j) = B_3(j; 0; 1)
#loop
M_j = matrix_hp(3; n_j)
c_j = vector_hp(n_j)
j = vector_hp(4)
#for e = 1 : n_e
    #for i = 1 : 4
        j.i = e_j.(e; i)
        c_j.j.i = c_j.j.i + 1
    #loop
    Z_е = Z_e(e)
    add(-D*B_0,0*Z_е; M_j; 1; j.1)
    add(-D*B_1,0*Z_е; M_j; 1; j.2)
    add(-D*B_1,1*Z_е; M_j; 1; j.3)
    add(-D*B_0,1*Z_е; M_j; 1; j.4)
#loop
#for j = 1 : n_j
     M_j.(1; j) = M_j.(1; j)/c_j.j
     M_j.(2; j) = M_j.(2; j)/c_j.j
     M_j.(3; j) = M_j.(3; j)/c_j.j
#loop
Mx = matrix_hp(n_ja; n_jb)
My = matrix_hp(n_ja; n_jb)
Mxy = matrix_hp(n_ja; n_jb)
#for i = 1 : n_ja
    #for k = 1 : n_jb
        j = (i - 1)*n_jb + k
        Mx.(i; k) = M_j.(1; j)
        My.(i; k) = M_j.(2; j)
        Mxy.(i; k) = M_j.(3; j)
    #loop
#loop
M_x(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mx)
M_y(x; y) = spline(1 + x/a_1; 1 + y/b_1; My)
M_xy(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mxy)
PlotShadows = 0
#show
'Average bending moments at joints, kNm/m
M_j
'Bending moments for the plate
'Bending moments - <var>M</var><sub>x</sub>
transp(Mx)'kNm/m
$Map{M_x(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments <var>M</var><sub>y</sub>
transp(My)'kNm/m
$Map{M_y(x; y) @ x = 0 : l_a & y = 0 : l_b}
'Bending moments <var>M</var><sub>xy</sub>
transp(Mxy)'kNm/m
$Map{M_xy(x; y) @ x = 0 : l_a & y = 0 : l_b}

Rectangular Slab FEA¶

Single-panel simply-supported rectangular slab solved with the 16-DOF Bogner–Fox–Schmit plate element under uniform load. Returns deflection plus the bending moments \(M_x\) and \(M_y\) on a user-defined mesh — a compact testbed for convergence studies.

Code:

'<h4>Input data</h4>
'<table><tr><td>
'Slab dimensions -'a = 6'm,'b = 4'm
'Thickness -'t = 0.1'm
'Load -'q = 10'kN/m²
'Modulus of elasticity -'E = 35000'MPa
'Poisson`s ratio -'ν = 0.15
'</td><td><img src="../../Images/mechanics/elastic/slab.png" style="height:110pt; margin-left:40pt;"></td></tr></table>
'<h4>Finite element mesh</h4>
'We will use rectangular finite element with 'n = 16' DOFs
'Number of elements along <var>a</var> and <var>b</var> -'n_a = 6', 'n_b = 4
'Total number of elements -'n_e = n_a*n_b
'Total number of joints -'n_j = (n_a + 1)*(n_b + 1)
'Element dimensions -'a_1 = a/n_a', 'b_1 = b/n_b
'Supported joints count -'n_s = 2*(n_a + n_b)
#hide
x_j = vector(n_j)','y_j = vector(n_j)
x = 0', 'y = 0
#for j = 1 : n_j
    x_j.j = x
    y_j.j = y
    y = y + b_1
    #if y > b
        y = 0
        x = x + a_1
    #end if
#loop
e_j = matrix(n_e; 4)
#for i_a = 1 : n_a
    #for i_b = 1 : n_b
        e = i_b + n_b*(i_a - 1)
        j = e + i_a - 1
        e_j.(e; 1) = j
        e_j.(e; 2) = j + n_b + 1
        e_j.(e; 3) = j + n_b + 2
        e_j.(e; 4) = j + 1
    #loop
#loop
s_j = vector(n_s)
i_s = 0
#for i = 1 : n_a + 1
    j = (n_b + 1)*i - n_b
    i_s = i_s + 1
    s_j.i_s = j
#loop
#for i = 1 : n_a + 1
    j = (n_b + 1)*i
    i_s = i_s + 1
    s_j.i_s = j
#loop
#for i = 2 : n_b
    j = i
    i_s = i_s + 1
    s_j.i_s = j
#loop
#for i = 2 : n_b
    j = n_a*(n_b + 1) + i
    i_s = i_s + 1
    s_j.i_s = j
#loop
#show
'Joint coordinates
x_j'm
y_j'm
'Numbers of elements joints
transp(e_j)
'Supported joints
s_j
'Coordinates of elements centers
j_e(e) = row(e_j; e)
x_c(e) = sum(extract(x_j; j_e(e)))/4', ' _
y_c(e) = sum(extract(y_j; j_e(e)))/4
#val
#hide
w = 400
k = w/a
d = 20
h = b/a*w
r = 0.04*k
#show
'<svg viewbox="'-d' '-d' 'w + 2*d' 'h + 2*d'" xmlns="http://www.w3.org/2000/svg" version="1.1" style=" font-family: Segoe UI; font-size:10px; width:'w'pt; height:'h'pt">
'<style>.joint{fill:orangeRed;} .support{stroke:red; stroke-width:1; fill:lightpink;} .element{stroke:seaGreen; stroke-width:1; fill:lime; fill-opacity:0.1; stroke-opacity:0.5}</style>
'<rect x="'0'" y="'0'" width="'w'" height="'h'" style="fill:yellow; fill-opacity:0.2" />
#for e = 1 : n_e
    #hide
    x = x_c(e)*k
    y = y_c(e)*k
    #show
    '<text x="'x'" y="'h - y'" text-anchor="middle">'e'</text>
    '<rect x="'x - a_1/2*k'" y="'h - y - b_1/2*k'" width="'a_1*k'" height="'b_1*k'" class="element" />
#loop
#for i = 1 : n_s
    j = s_j.i
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    #if y_j.j ≡ 0 ∨ y_j.j ≡ b
        '<line x1="'x - 4*r'" y1="'y'" x2="'x + 4*r'" y2="'y'" class="support"/>
    #end if
    #if x_j.j ≡ 0 ∨ x_j.j ≡ a
        '<line x1="'x'" y1="'y - 4*r'" x2="'x'" y2="'y + 4*r'" class="support"/>
    #end if
    '<circle cx="'x'" cy="'y'" r="'2*r'" class="support"/>
#loop
#for j = 1 : n_j
    #hide
    x = x_j.j*k
    y = h - y_j.j*k
    #show
    '<circle cx="'x'" cy="'y'" r="'r'" class="joint" />
    '<text x="'x + 2*r'" y="'y - r'" text-anchor="start">'j'</text>
#loop
'</svg>
#equ
'<h4>Finite element formulation</h4>
'<p><b>Shape functions</b></p>
#def FI$(x$; l$)
    '<table><tr><td>Base functions
    Φ_1l$(x$) = 1 - x$^2*(3 - 2*x$)
    Φ_2l$(x$) = x$*l$_1*(1 - x$*(2 - x$))
    Φ_3l$(x$) = x$^2*(3 - 2*x$)
    Φ_4l$(x$) = x$^2*l$_1*(-1 + x$)
    '</td><td>First derivatives
    Φ′_1l$(x$) = -6*(x$/l$_1)*(1 - x$)
    Φ′_2l$(x$) = 1 - x$*(4 - 3*x$)
    Φ′_3l$(x$) = 6*(x$/l$_1)*(1 - x$)
    Φ′_4l$(x$) = -x$*(2 - 3*x$)
    '</td><td>Second derivatives
    Φ″_1l$(x$) = -(6/l$_1^2)*(1 - 2*x$)
    Φ″_2l$(x$) = -(2/l$_1)*(2 - 3*x$)
    Φ″_3l$(x$) = 6/l$_1^2*(1 - 2*x$)
    Φ″_4l$(x$) = -(2/l$_1)*(1 - 3*x$)
    '</td></tr></table>
#end def
'Along dimension <var>a</var>
FI$(ξ; a)
'Along dimension <var>b</var>
FI$(η; b)
'<table><tr><td>For vertical displacements <var>w</var>
N_1,w(ξ; η) = Φ_1a(ξ)*Φ_1b(η)
N_2,w(ξ; η) = Φ_3a(ξ)*Φ_1b(η)
N_3,w(ξ; η) = Φ_3a(ξ)*Φ_3b(η)
N_4,w(ξ; η) = Φ_1a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ₓ
N_1,θₓ(ξ; η) = Φ_2a(ξ)*Φ_1b(η)
N_2,θₓ(ξ; η) = Φ_4a(ξ)*Φ_1b(η)
N_3,θₓ(ξ; η) = Φ_4a(ξ)*Φ_3b(η)
N_4,θₓ(ξ; η) = Φ_2a(ξ)*Φ_3b(η)
'</td><td>For rotations <var>θ</var>ᵧ
N_1,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_2b(η)
N_2,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_2b(η)
N_3,θᵧ(ξ; η) = Φ_3a(ξ)*Φ_4b(η)
N_4,θᵧ(ξ; η) = Φ_1a(ξ)*Φ_4b(η)
'</td></tr></table>
'<img src="./plate-element.png" alt="plane-elemen.png" style="float:right; width:200pt;" />
'For twist ψ
N_1,ψ(ξ; η) = Φ_2a(ξ)*Φ_2b(η)
N_2,ψ(ξ; η) = Φ_4a(ξ)*Φ_2b(η)
N_3,ψ(ξ; η) = Φ_4a(ξ)*Φ_4b(η)
N_4,ψ(ξ; η) = Φ_2a(ξ)*Φ_4b(η)
'Shape functions vector
N(i; ξ; η) = take(i; _
N_1,w(ξ; η); N_1,θₓ(ξ; η); N_1,θᵧ(ξ; η); N_1,ψ(ξ; η); _
N_2,w(ξ; η); N_2,θₓ(ξ; η); N_2,θᵧ(ξ; η); N_2,ψ(ξ; η); _
N_3,w(ξ; η); N_3,θₓ(ξ; η); N_3,θᵧ(ξ; η); N_3,ψ(ξ; η); _
N_4,w(ξ; η); N_4,θₓ(ξ; η); N_4,θᵧ(ξ; η); N_4,ψ(ξ; η))
'<p><b>Constitutive matrix</b> (stress - strain relationship)</p>
D = E*t^3/(12*(1 - ν^2))*[1; ν; 0|ν; 1; 0|0; 0; (1 - ν)/2]
'<p><b>Strain-displacement matrix</b></p>
B_1(j; ξ; η) = take(j; _
Φ″_1a(ξ)*Φ_1b(η); Φ″_2a(ξ)*Φ_1b(η); Φ″_1a(ξ)*Φ_2b(η); Φ″_2a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_1b(η); Φ″_4a(ξ)*Φ_1b(η); Φ″_3a(ξ)*Φ_2b(η); Φ″_4a(ξ)*Φ_2b(η); _
Φ″_3a(ξ)*Φ_3b(η); Φ″_4a(ξ)*Φ_3b(η); Φ″_3a(ξ)*Φ_4b(η); Φ″_4a(ξ)*Φ_4b(η); _
Φ″_1a(ξ)*Φ_3b(η); Φ″_2a(ξ)*Φ_3b(η); Φ″_1a(ξ)*Φ_4b(η); Φ″_2a(ξ)*Φ_4b(η))
B_2(j; ξ; η) = take(j; _
Φ_1a(ξ)*Φ″_1b(η); Φ_2a(ξ)*Φ″_1b(η); Φ_1a(ξ)*Φ″_2b(η); Φ_2a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_1b(η); Φ_4a(ξ)*Φ″_1b(η); Φ_3a(ξ)*Φ″_2b(η); Φ_4a(ξ)*Φ″_2b(η); _
Φ_3a(ξ)*Φ″_3b(η); Φ_4a(ξ)*Φ″_3b(η); Φ_3a(ξ)*Φ″_4b(η); Φ_4a(ξ)*Φ″_4b(η); _
Φ_1a(ξ)*Φ″_3b(η); Φ_2a(ξ)*Φ″_3b(η); Φ_1a(ξ)*Φ″_4b(η); Φ_2a(ξ)*Φ″_4b(η))
B_3(j; ξ; η) = 2*take(j; _
Φ′_1a(ξ)*Φ′_1b(η); Φ′_2a(ξ)*Φ′_1b(η); Φ′_1a(ξ)*Φ′_2b(η); Φ′_2a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_1b(η); Φ′_4a(ξ)*Φ′_1b(η); Φ′_3a(ξ)*Φ′_2b(η); Φ′_4a(ξ)*Φ′_2b(η); _
Φ′_3a(ξ)*Φ′_3b(η); Φ′_4a(ξ)*Φ′_3b(η); Φ′_3a(ξ)*Φ′_4b(η); Φ′_4a(ξ)*Φ′_4b(η); _
Φ′_1a(ξ)*Φ′_3b(η); Φ′_2a(ξ)*Φ′_3b(η); Φ′_1a(ξ)*Φ′_4b(η); Φ′_2a(ξ)*Φ′_4b(η))
B(j; ξ; η) = [B_1(j; ξ; η); B_2(j; ξ; η); B_3(j; ξ; η)]
x_1(e) = x_j.e_j.(e; 1)','y_1(e) = y_j.e_j.(e; 1)
'The elements of the stiffness matrix will be calculated by using the equation
#noc
K_e,ij = a_1*b_1*$Area{$Area{B_i(ξ; η)^T*D*B_j(ξ; η) @ ξ = 0 : 1} @ η = 0 : 1}
#equ
'<p><b>Element stiffness matrix</b></p> (above the main diagonal only)
#hide
K_e = utriang(n)
Precision = 10^-4
e = 1
#show
BTDB_e(i; j; ξ; η) = transp(B(i; ξ; η))*D*B(j; ξ; η)
K_e(i; j) = a_1*b_1*$Area{$Area{BTDB_e(i; j; ξ; η) @ η = 0 : 1} @ ξ = 0 : 1}
$Repeat{$Repeat{K_e.(i; j) = K_e(i; j) @ j = i : n} @ i = 1 : n}
K_e
'Element load vector
#noc
F_e,i = a_1*b_1*$Area{$Area{N_i(ξ; η)^T*q @ ξ = 0 : 1} @ η = 0 : 1}
#equ
#hide
F_e = vector(n)
#for i = 1 : n
    F_e.i = a_1*b_1*$Area{$Area{N(i; ξ; η)*q @ ξ = 0 : 1} @ η = 0 : 1}
#loop
#show
F_e'kN
'<h4>Solution</h4>
'Global stiffness matrix
#hide
k_1 = n/4
n = k_1*n_j
K = symmetric(n)
'Addition of element stiffness matrix coefficients
#for e = 1 : n_e
    #for i = 1 : 4
        #for j = i : 4
            j_1 = e_j.(e; i)
            j_2 = e_j.(e; j)
            i_1 = k_1*(j_1 - 1) + 1
            i_2 = k_1*(j_2 - 1) + 1
            k_ij = submatrix(K_e; k_1*(i - 1) + 1; k_1*i; k_1*(j - 1) + 1; k_1*j)
            #if j_1 > j_2
                add(transp(k_ij); K; i_2; i_1)
            #else
                add(k_ij; K; i_1; i_2)
            #end if
        #loop
    #loop
#loop
k_s = 10^20
'Addition of supports
#for i = 1 : n_s
    j = k_1*(s_j.i - 1) + 1
    K.(j; j) = K.(j; j) + k_s
    j = j + 1
    #if y_j.s_j.i ≡ 0 ∨ y_j.s_j.i ≡ b
        K.(j; j) = K.(j; j) + k_s
    #end if
    j = j + 1
    #if x_j.s_j.i ≡ 0 ∨ x_j.s_j.i ≡ a
        K.(j; j) = K.(j; j) + k_s
    #end if
#loop
#show
K
'Global load vector
#hide
F = vector(n)
#for e = 1 : n_e
    #for i = 1 : 4
        j = e_j.(e; i)
        #for k = 1 : 4
            l = k_1*(j - 1) + k
            F.l = F.l + F_e.(k_1*(i - 1) + k)
        #loop
    #loop
#loop
#show
F'kN
'Solution of the system of equations
Z = clsolve(K; F)'mm
'<h4>Results</h4>
'Joint displacements
#hide
W_z = matrix(n_a + 1; n_b + 1)
w(j) = Z.(4*j - 3)
$Repeat{$Repeat{W_z.(i; j) = round(w((i - 1)*(n_b + 1) + j)*1000)/1000 @ j = 1 : n_b + 1} @ i = 1 : n_a + 1}
w(x; y) = spline(1 + x/a_1; 1 + y/b_1; W_z)
PlotWidth = 50*a
PlotHeight = 50*b
PlotStep = 10
#show
transp(W_z)'mm
$Map{-w(x; y) @ x = 0 : a & y = 0 : b}
w(a/2; b/2)'mm
'Bending moments
Z_j(j) = slice(Z; k_1*(j - 1) + 1; k_1*j)
Z_e(e) = [Z_j(e_j.(e; 1)); Z_j(e_j.(e; 2)); Z_j(e_j.(e; 3)); Z_j(e_j.(e; 4))]
#hide
B(ξ; η) = [ _
B_1(1; ξ; η); B_1(2; ξ; η); B_1(3; ξ; η); B_1(4; ξ; η); _
B_1(5; ξ; η); B_1(6; ξ; η); B_1(7; ξ; η); B_1(8; ξ; η); _
B_1(9; ξ; η); B_1(10; ξ; η); B_1(11; ξ; η); B_1(12; ξ; η); _
B_1(13; ξ; η); B_1(14; ξ; η); B_1(15; ξ; η); B_1(16; ξ; η)| _
B_2(1; ξ; η); B_2(2; ξ; η); B_2(3; ξ; η); B_2(4; ξ; η); _
B_2(5; ξ; η); B_2(6; ξ; η); B_2(7; ξ; η); B_2(8; ξ; η); _
B_2(9; ξ; η); B_2(10; ξ; η); B_2(11; ξ; η); B_2(12; ξ; η); _
B_2(13; ξ; η); B_2(14; ξ; η); B_2(15; ξ; η); B_2(16; ξ; η)| _
B_3(1; ξ; η); B_3(2; ξ; η); B_3(3; ξ; η); B_3(4; ξ; η); _
B_3(5; ξ; η); B_3(6; ξ; η); B_3(7; ξ; η); B_3(8; ξ; η); _
B_3(9; ξ; η); B_3(10; ξ; η); B_3(11; ξ; η); B_3(12; ξ; η); _
B_3(13; ξ; η); B_3(14; ξ; η); B_3(15; ξ; η); B_3(16; ξ; η)]
M_j = matrix(3; n_j)
c_j = vector(n_j)
#for e = 1 : n_e
    Z_е = Z_e(e)
    j_1 = e_j.(e; 1)
    x_1 = x_j.j_1
    y_1 = y_j.j_1
    M_e(x; y) = -D*B(x/a_1; y/b_1)*Z_е
    #for i = 1 : 4
        j = e_j.(e; i)
        c_j.j = c_j.j + 1
        x = x_j.j - x_1
        y = y_j.j - y_1
        M_e = M_e(x; y)
        add(M_e; M_j; 1; j)
    #loop
#loop
#for j = 1 : n_j
     M_j.(1; j) = M_j.(1; j)/c_j.j
     M_j.(2; j) = M_j.(2; j)/c_j.j
     M_j.(3; j) = M_j.(3; j)/c_j.j
#loop
Mx = matrix(n_a + 1; n_b + 1)
My = matrix(n_a + 1; n_b + 1)
Mxy = matrix(n_a + 1; n_b + 1)
#for i = 1 : n_a + 1
    #for k = 1 : n_b + 1
        j = (i - 1)*(n_b + 1) + k
        Mx.(i; k) = M_j.(1; j)
        My.(i; k) = M_j.(2; j)
        Mxy.(i; k) = M_j.(3; j)
    #loop
#loop
M_x(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mx)
M_y(x; y) = spline(1 + x/a_1; 1 + y/b_1; My)
M_xy(x; y) = spline(1 + x/a_1; 1 + y/b_1; Mxy)
i = round(n_a/2)
j = round(n_b/2)
e = (i - 1)*n_a + j + 1
j = e_j.(e; 1)
x = 0
y = 0
#show
#val
'Results for element 'e' and joint 'j':
#equ
Z_e = Z_e(e)'mm
M_e(x; y) = -D*B(x/a_1; y/b_1)*Z_e(e)
M_e = M_e(0; 0)'kNm/m
'Average bending moments at joints, kNm/m
M_j
'Bending moments for the plate
'Bending moments - <var>M</var><sub>x</sub>
transp(Mx)'kNm/m
$Map{M_x(x; y) @ x = 0 : a & y = 0 : b}
'Maximal value -'M_x(a/2; b/2)'kNm/m
'Bending moments <var>M</var><sub>y</sub>
transp(My)'kNm/m
$Map{M_y(x; y) @ x = 0 : a & y = 0 : b}
'Maximal value -'M_y(a/2; b/2)'kNm/m
'Bending moments <var>M</var><sub>xy</sub>
transp(Mxy)'kNm/m
$Map{M_xy(x; y) @ x = 0 : a & y = 0 : b}
'Maximal value -'M_xy(0; 0)'kNm/m

Spotted an error? Edit these examples.

`N`_1,w shape function plot	`N`_1,θₓ shape function plot

`N`_1,w shape function plot	`N`_1,θₓ shape function plot