"Finite Element Analysis of Flat Slab<hr/>
'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}