"Material and geometric nonliner analysis of prestressed asymmetric Biot truss
'<h4>Input data</h4>
#deg
#hide
δz = 10^-12
Precision = 10^-9
x_J = [0; 3; 9]*m
y_J = [0; 0; 0]*m
#show
x_J','y_J', 'n_J = len(x_J)
'<h4>Elements - [J1; J2]</h4>
#hide
E_J = [1; 2|2; 3]
#show
transp(E_J)', 'n_E = n_rows(E_J)', 'J_1(e) = E_J.(e; 1)', 'J_2(e) = E_J.(e; 2)
'Element endpoint coordinates
x_1(e) = x_J.J_1(e)','y_1(e) = y_J.J_1(e)', 'x_2(e) = x_J.J_2(e)','y_2(e) = y_J.J_2(e)
'Element lengths - 'l(e) = sqrt((x_2(e) - x_1(e))^2 + (y_2(e) - y_1(e))^2)
'Element directions -'c(e) = (x_2(e) - x_1(e))/l(e)','s(e) = (y_2(e) - y_1(e))/l(e)
'Transformation matrix
t(e) = [c(e); s(e)| - s(e); c(e)]
T(e) = add(t(e); add(t(e); matrix(4; 4); 1; 1); 3; 3)
'<h4>Supports - [Joint; cx; cy]</h4>
#hide
c_J = [1; 10^20kN/m; 10^20kN/m|3; 10^20kN/m; 10^20kN/m]
#show
c_J
n_c = n_rows(c_J)
'<h4>Material - steel</h4>
'We assume bi-linear material model 
'Initial modulus of elasticity -'E_0 = 206GPa
'Yield stress -'f_y = 500MPa
'Ultimate tensile stress -'f_u = 600MPa
'Yield strain -'ε_y = f_y/E_0|‰
'Ultimate  strain -'ε_u = 20‰
'Modulus of elasticity after yield -'E_1 = (f_u - f_y)/(ε_u - ε_y)|GPa
'Idealized stress-strain curve -'σ(ε) = if(ε < ε_y; E_0*ε; f_y + E_1*(ε - ε_y))
'<!--'PlotWidth = 200','PlotHeight = 150'-->
$Plot{σ(ε) @ ε = 0‰ : ε_u}
'<h4>Cross section – circular with diameter'Φ = 20mm'</h4>
'Area -'A = π*Φ^2/4|cm^2
'Stiffness
'Initial -'EA_0 = E_0*A', 'EA = EA_0
'After yield -'EA_1 = E_1*A
'<h4>Load - [Joint, Fx, Fy]</h4>
#hide
F_J = join_rows([2; 0kN; -70kN])
#show
F_J', 'n_F = n_rows(F_J)
'Prestressing force -'P = 20kN
'Initial stress -'σ_P = P/A|MPa
'Initial strain -'ε_P = σ_P/E_0
#include svg_drawing.cpd
#hide
Z = vector(2*n_J)*mm
w = max(x_J)
h = max(y_J)
W = 320
k = W/w
H = (h + 3m)*k
r = 3
x1(j) = (x_J.j + Z.(2*j - 1))*k
y1(j) = (h - y_J.j + Z.(2*j))*k
#def svg$ = '<svg viewbox="'-2.5m*k' '-1m*k' '(w + 5m)*k' 'H'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro Cond; font-size:11px; width:'W + 100'pt; height:'H + 100*H/W'pt">
#def thin_style$ = style="stroke:darkcyan; stroke-width:1; fill:none"
#def pin_style$ = style="stroke:red; stroke-width:1; fill:white"
#def thick_style$ = style="stroke:darkcyan; stroke-width:3; fill:none"
#def load_style$ = style="stroke:deepskyblue; stroke-width:1; fill:deepskyblue;"
#show
#val
svg$
#def frame$
	#for i = 1 : n_E
		#hide
		x1 = x1(J_1(i))
		y1 = y1(J_1(i))
		x2 = x1(J_2(i))
		y2 = y1(J_2(i))
		#show
		line$(x1; y1; x2; y2; main_style$)
		dimh$(x1;x2;H-1.5m*k;'(x2-x1)/k'm)
	#loop
	#for i = 1 : n_c
		#hide
		j = c_J.(i; 1)
		x1 = x1(j)
		y1 = y1(j) + r
		δ = w/30*k*sign(x1 - w/2*k)
		x2 = x1 - δ
		y2 = y1 - abs(δ)
		x3 = x1 + δ
		y3 = y1 + abs(δ)
		#show
		line$(x2; y3; x3; y3; thick_style$)
		line$(x2; y3; x1; y1; thin_style$)
		line$(x3; y3; x1; y1; thin_style$)
	#loop
	#for i = 1 : n_F
		#hide
		δ = 0.5m*k
		j = F_J.(i; 1)
		F_x = F_J.(i; 2)/kN
		F_y = F_J.(i; 3)/kN
		#show
		#if F_x ≠ 0
			#hide
			δ_x = 3*δ*sign(F_x)
			δ_y = 0.2*δ
			x1 = x1(j) + r*sign(F_x)
			y1 = y1(j)
			#show
			'<polygon points="'x1 - δ_x','y1' 'x1','y1' 'x1 - δ_x/5','y1 - δ_y' 'x1 - δ_x/5','y1 + δ_y' 'x1','y1'" load_style$ />
			texth$(x1-2*δ_x/3;y1-2*δ_y;Fx='F_x'kN)
		#end if
		#if F_y ≠ 0
			#hide
			δ_x = 0.2*δ
			δ_y = 3*δ*sign(F_y)
			x1 = x1(j)
			y1 = y1(j) + r*sign(F_y)
			#show
			'<polygon points="'x1','y1 + δ_y' 'x1','y1' 'x1 - δ_x','y1 + δ_y/5' 'x1 + δ_x','y1 + δ_y/5' 'x1','y1'" load_style$ />
			textv$(x1-2*δ_x;y1+2*δ_y/3;Fy='F_y'kN)
		#end if
	#loop
	#for j = 1 : n_J
		circle$(x1(j); y1(j); r; pin_style$)
		#if abs(Z.(2*j)) > δz*mm
			textv$(x1(j)-0.2*δ; y1(j)+1.8*δ;'Z.(2*j)'mm)
		#end if
		#if abs(Z.(2*j - 1)) > δz*mm
			texth$(x1(j)+1.5*δ; y1(j)+0.8*δ;'Z.(2*j - 1)'mm)
		#end if
	#loop
#end def
frame$
'</svg>
#equ
'<h4>Finite element formulation</h4>
'<!--'u = 0','v = 0','z = [0; 0]'--> 
'We will formulate a 2D truss element with 3rd order geometric nonlinearity. The equilibrium equations will be derived in local CS, for the deformed state of the structure, assuming large displacements, as follows:
'<img style="height:115pt; width:198pt;" src="./Images/nonlin-truss-element.png" alt="nonlin-truss-element.png">
#noc
'Relative displacements -'u = u_2 - u_1', 'v = v_2 - v_1', 'z = [u; v]
#equ
'Length and elongation of the element in deformed state
l′(e; z) = sqrt((l(e) + z.1)^2 + z.2^2)','Δl(e; z) = l′(e; z) - l(e)', 'ε(e; z) = Δl(e; z)/l(e) + ε_P
'Directions of the displaced axis of the element in the deformed state
c_d(e; z) = (l(e) + z.1)/l′(e; z)', 's_d(e; z) = z.2/l′(e; z)
'Axial force in element -'N(e; z) = σ(ε(e; z))*A
'Horizontal reactive force at the right end -'F_x(e; z) = N(e; z)*c_d(e; z)
'Vertical reactive force at the right end -'F_y(e; z) = N(e; z)*s_d(e; z)
'Partial derivatives of element end reactive forces
F′_x1(e; z) = (1/l(e) - z.2^2/l′(e; z)^3)', ' _
F′_x2(e; z) = (z.2*(l(e) + z.1)/l′(e; z)^3)
F′_y1(e; z) = F′_x2(e; z)', ' _
F′_y2(e; z) = (1/l(e) - (l(e) + z.1)^2/l′(e; z)^3)
'We are not going to linearize the above equations by Taylor series. We will use them further directly in their implicit form.
'<h4>Solution by Newton-Raphson’s method</h4>
'We will express the parameters of the system as a function of the vector of t he unknown displacements:
Z = vector(2*n_J)*mm
'Initial deflection at the intermediate joint -'Z.(n_J + 1) = 100mm
'Extracting displacements for joint <var>j</var> -'z_J(j) = [Z.(2*j - 1); Z.(2*j)]
'Relative displacements for element <var>e</var> -' _
z_e(e) = t(e)*(z_J(J_2(e)) - z_J(J_1(e)))
'Element reactive forces partial vector -' _
F_e(e) = t(e)*[F_x(e; z_e(e)); F_y(e; z_e(e))]
'Partial element Jacobi matrix
J_e(e) = t(e)*[ _
F′_x1(e; z_e(e)); F′_x2(e; z_e(e))| _
F′_y1(e; z_e(e)); F′_y2(e; z_e(e))]
'Target precision -'ε_max = 10^-5
#hide
#for n = 1 : 100
	z_J(j) = [Z.(2*j - 1); Z.(2*j)]
	'Global load vector
	F = vector(2*n_J)*kN
	s = -1
	'Add internal forces
	#for e = 1 : n_E
		#for k = 1 : 2
			j = E_J.(e; k)
			j_2 = 2*j
			F.(j_2 - 1) = F.(j_2 - 1) + take(1; F_e(e))*s
			F.(j_2) = F.(j_2) + take(2; F_e(e))*s
			s = -s
		#loop
	#loop
	'Add support rections
	#for i = 1 : n_c
		j = c_J.(i; 1)
		j_2 = 2*j
		F.(j_2 - 1) = F.(j_2 - 1) + c_J.(i; 2)*Z.(j_2 - 1)
		F.(j_2) = F.(j_2) + c_J.(i; 3)*Z.(j_2)
	#loop
	'Add external loads
	#for i = 1 : n_F
		j = F_J.(i; 1)
		j_2 = 2*j
		F.(j_2 - 1) = F.(j_2 - 1) + F_J.(i; 2)
		F.(j_2) = F.(j_2) + F_J.(i; 3)
	#loop
	'Global Jacobi matrix
	J = symmetric(2*n_J)*kN/m
	'Add element Jacobi matrices
	#for e = 1 : n_E
		ε = ε(e; z_e(e))
		EA = if(ε < ε_y; EA_0; max(EA_0/10; EA_1))
		#for i = 1 : 2
			#for j = 1 : 2
				k = if(i ≡ j; EA; -EA)
				i_1 = 2*E_J.(e; i) - 1
				j_1 = 2*E_J.(e; j) - 1
				add(k*J_e(e); J; i_1; j_1)
			#loop
		#loop
	#loop
	'Add support springs
	#for i = 1 : n_c
		j = c_J.(i; 1)
		C_x = c_J.(i; 2)
		C_y = c_J.(i; 3)
		j_1 = 2*j - 1
		j_2 = 2*j
		J.(j_1; j_1) = J.(j_1; j_1) + C_x
		J.(j_2; j_2) = J.(j_2; j_2) + C_y
	#loop
	Z_0 = Z
	δZ = clsolve(J; F)|mm
	Z = Z_0 - δZ
	ε = norm_e(δZ)/norm_e(Z)
	#if ε < ε_max
		#break
	#end if
#loop
#show
'Iteratively calculate the following equations of the Newton-Raphson’s method:
#noc
Z_0 = Z','δZ = clsolve(J(Z); F(Z))|mm','Z = Z_0 - δZ','ε = norm_e(δZ)/norm_e(Z) < ε_max
#equ
'Convergence reached at iteration -'n'. Relative error -'ε
'<h4>Results</h4>
'<p><strong>Joint displacements</strong> -'Z(j) = round(z_J(j)/δz)*δz'</p>
#novar
#for j = 1 : n_J
	#if norm_1(Z(j)) > δz*mm
		#val
		'<p>Joint <b>J'j'</b> - <var>Z</var><sub>'j
		#equ
		'</sub> ='Z(j)'</p>
	#end if
#loop
#varsub
'<p><strong>Reactions in supports</strong></p>
r(i) = row(c_J; i)','j(i) = take(1; r(i))', ' _
R(i) = -z_J(j(i))*last(r(i); 2)
#novar
#for i = 1 : n_c
	#val
	'<p>Joint <b>J'j(i)'</b> - <var>R</var><sub>'i
	#equ
	'</sub> ='R(i)'</p>
#loop
#varsub
'Element yield capacity  -'N_y = f_y*A
'Element ultimate capacity -'N_u = f_u*A
'<p><strong>Element axial forces</strong> -'N_e(e) = N(e; z_e(e))'</p> 
#novar
#for e = 1 : n_E
	#val
	'<p>Element <b>Е'e'</b> - <var>N</var><sub>'e'
	#equ
	'</sub> ='N_e(e)'</p>
#loop
#varsub
#val
#hide
k = W/w
δ = 0.5m*k
#show
svg$
'<polygon points="'0','0' '0','-3*δ' '-0.2*δ','-2.4*δ' '0.2*δ','-2.4*δ' '0','-3*δ'" load_style$/>
'<polygon points="'W','0' 'W','-3*δ' 'W - 0.2*δ','-2.4*δ' 'W + 0.2*δ','-2.4*δ' 'W','-3*δ'" load_style$/>
'<polygon points="'r','0' '-3*δ','0' '-2.4*δ','-0.2*δ' '-2.4*δ','0.2*δ' '-3*δ','0'" load_style$/>
'<polygon points="'W + r','0' 'W + 3*δ','0' 'W + 2.4*δ','-0.2*δ' 'W + 2.4*δ','0.2*δ' 'W + 3*δ','0'" load_style$/>
texth$(-2*δ;-0.5*δ;Rx='take(1; R(1))'kN)
texth$(W+2*δ;-0.5*δ;Rx='take(1; R(2))'kN)
texth$(2*δ;-3.2*δ;Ry='take(2; R(1))'kN)
texth$(W-2*δ;-3.2*δ;Ry='take(2; R(2))'kN)
frame$
'</svg>
#equ