# from https://github.com/seamplex/pipe-linearize

# convergence study of linearized stresses in an infinite pipe
# with respect to the number of elements in the pipe thickness
PROBLEM mechanical PC mumps
# PETSC_OPTIONS -mg_levels_pc_type sor

# t0 = wall_time()     # initial wall time

# read mesh according to shape $1, order $2 and number of elements through thickness $2
DEFAULT_ARGUMENT_VALUE 2 2
DEFAULT_ARGUMENT_VALUE 3 2
READ_MESH pipe-$1-$2-$3.msh 

# problem parameters for
# 12"-inch schedule 100
b = 323.8/2    # external radius [ mm ]
a = b-21.5     # internal radius [ mm ]
l = 2*(b-a)    # axial length [ mm ]

E = 200e3      # Young modulus [ MPa ]
nu = 0.3       # Poisson's ratio [ non-dimensional ]
p = 10         # internal pressure [ MPa ]


# ------------------------------------------------------------------------
# definition of analytical solutions for comparison from
# <http://eprints.whiterose.ac.uk/110536/1/art%253A10.1007%252Fs00707-016-1762-7.pdf>
ur(x,y,z) = (p*a^2*sqrt(y^2+z^2))/(E*(b^2-a^2)) * ((1-2*nu)*(1+nu) + (1+nu)*b^2/(y^2+z^2))
sigmal(x,y,z)     = 2*nu*p*a^2/(b^2-a^2)
sigmar(x,y,z)     = p*a^2/(b^2-a^2) * (1 - b^2/(y^2+z^2))
sigmatheta(x,y,z) = p*a^2/(b^2-a^2) * (1 + b^2/(y^2+z^2))


# principal stresses along the radial coordinate (may be y or z)
sigma1_anal(r) = sigmatheta(0,0,r)
sigma2_anal(r) = sigmal(0,0,r)
sigma3_anal(r) = sigmar(0,0,r)

# computation of main membrane stresses
M1 = 1/(b-a)*integral(sigma1_anal(r), r, a, b)
M2 = 1/(b-a)*integral(sigma2_anal(r), r, a, b)
M3 = 1/(b-a)*integral(sigma3_anal(r), r, a, b)

# von mises and tresca membrane stresses
Mt_anal = max(abs(M1-M2), abs(M2-M3), abs(M3-M1))
Mv_anal = sqrt(((M1-M2)^2 + (M2-M3)^2 + (M3-M1)^2)/2)

# computation of membrane plus bending stresses
MB1 = M1 + 6/(b-a)^2*integral(sigma1_anal(r)*((a+b)/2-r), r, a, b)
MB2 = M2 + 6/(b-a)^2*integral(sigma2_anal(r)*((a+b)/2-r), r, a, b)
MB3 = M3 + 6/(b-a)^2*integral(sigma3_anal(r)*((a+b)/2-r), r, a, b)

# von mises and tresca
MBv_anal = sqrt(((MB1-MB2)^2 + (MB2-MB3)^2 + (MB3-MB1)^2)/2)
MBt_anal = max(abs(MB1-MB2), abs(MB2-MB3), abs(MB3-MB1))



# set boundary conditions
BC pressure  pressure=p
BC front     tangential radial
BC back      tangential radial

SOLVE_PROBLEM

# write distribution of results in gmsh format (optional)
# WRITE_MESH pipe-out-$1-$2-$3.msh VECTOR u v w sigma1 sigma2 sigma3 sigmax sigmay sigmaz tauxy tauyz tauzx

# write same thing as ASCII data
# PRINT_FUNCTION v sigma1 sigma2 sigma3  MIN 0 a 0 MAX 0 b 0 NSTEPS 1 200 1  FILE pipe-dist-$1-$2-$3.dat

# compute linearized stresses for both von Mises and Tresca
LINEARIZE_STRESS FROM 0 a 0 TO 0 b 0     M Mv   MB MBv        Mt Mt   MBt MBt

# compute L2 errors
volume = pi*(b^2-a^2)*l
h = (volume/cells)^(1/3)

ur_fea(x,y,z) = sqrt(v(x,y,z)^2+w(x,y,z)^2)
INTEGRATE (ur_fea(x,y,z)-ur(x,y,z))^2         RESULT e2ur
INTEGRATE (sigma1(x,y,z)-sigmatheta(x,y,z))^2 RESULT e2sigma1
INTEGRATE (sigma3(x,y,z)-sigmar(x,y,z))^2     RESULT e2sigma3

error_ur = sqrt(e2ur)/volume
error_sigma1 = sqrt(e2sigma1)/volume
error_sigma3 = sqrt(e2sigma3)/volume
error_Mv = abs(Mv-Mv_anal)/Mv_anal
error_MBv = abs(MBv-MBv_anal)/Mv_anal
error_Mt = abs(Mt-Mt_anal)/Mv_anal
error_MBt = abs(MBt-MBt_anal)/Mv_anal

PRINT error_ur+error_sigma1+error_sigma3+error_Mv+error_MBv+error_Mt+error_MBt