# 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 ]
n = 4          # number of elements through thickness 

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))


# analytical results in [ MPa ]	
M_tresca    = 	69.9467
M_vonmises  = 	61.7785
M_1         = 	65.3023
M_2         = 	18.1974
M_3         = 	-4.64439

MB_tresca   = 	79.9062
MB_vonmises = 	70.2562
MB_1        = 	70.2821
MB_2        = 	18.1974
MB_3        = 	-9.62415


ur1(x) = ur(0,x,0)
sigmal1(x) = sigmal(0,x,0)
sigmar1(x) = sigmar(0,x,0)
sigmatheta1(x) = sigmatheta(0,x,0)

delta = 0.05

set preamble "\usepackage{amsmath}"

set width 14*unit(cm)

set axis x arrow nomirrored
set axis y arrow nomirrored


set grid
set terminal pdf


set output "ur.pdf"
set xrange [a:b]
set yrange [*:*]
set key top right
set xlabel "$r$ [mm]"
set ylabel "$u_r(r)$ [mm]"
plot ur1(x)                        w l  lw 3  lt 1   color olivegreen   ti "Analytical solution $u_r(r)$",\
     "pipe-dist-1-3.dat" u 2:4 w l  lw 1  lt 2   color maroon       ti "",\
     "pipe-dist-1-3.dat" u 2:4 w p  pt 18 ps 1.5 color maroon       every (100)/(1*3) ti "$n=3$ first order",\
     "pipe-dist-1-6.dat" u 2:4 w l  lw 1  lt 4   color tan          ti "",\
     "pipe-dist-1-6.dat" u 2:4 w p  pt 16 ps 1.5 color tan          every (100)/(1*6) ti "$n=6$ first order",\
     "pipe-dist-2-2.dat" u 2:4 w l  lw 2  lt 3   color cerulean     ti "",\
     "pipe-dist-2-2.dat" u 2:4 w p  pt 17 ps 1.5 color cerulean     every (100)/(2*2) ti "$n=2$ second order"

set output "sigmatheta.pdf"
set xrange [a:b]
set yrange [*:*]
set key top right
set xlabel "$r$ [mm]"
set ylabel "$\sigma_\theta(r)$ [mm]"
plot sigmatheta1(x)                w l  lw 3  lt 1   color springgreen ti "Analytical solution $\sigma_\theta(r)$",\
     "pipe-dist-1-3.dat" u 2:5 w l  lw 1  lt 2   color melon       ti "",\
     "pipe-dist-1-3.dat" u 2:5 w p  pt 18 ps 1.5 color melon       every (100)/(1*3) ti "$n=3$ first order",\
     "pipe-dist-1-6.dat" u 2:5 w l  lw 1  lt 2   color redorange   ti "",\
     "pipe-dist-1-6.dat" u 2:5 w p  pt 16 ps 1.5 color redorange   every (100)/(1*6) ti "$n=6$ first order",\
     "pipe-dist-2-2.dat" u 2:5 w l  lw 2  lt 3   color turquoise   ti "",\
     "pipe-dist-2-2.dat" u 2:5 w p  pt 17 ps 1.5 color turquoise   every (100)/(2*2) ti "$n=2$ second order"

set output "sigmal.pdf"
set xrange [a:b]
set yrange [*:*]
set key top right
set xlabel "$r$ [mm]"
set ylabel "$\sigma_l(r)$ [mm]"
plot sigmal1(x)                    w l  lw 3  lt 1   color forestgreen    ti "Analytical solution $\sigma_l(r)$",\
     "pipe-dist-1-3.dat" u 2:6 w l  lw 1  lt 2   color peach          ti "",\
     "pipe-dist-1-3.dat" u 2:6 w p  pt 18 ps 1.5 color peach          every (100)/(1*3) ti "$n=3$ first order",\
     "pipe-dist-1-6.dat" u 2:6 w l  lw 1  lt 2   color dandelion      ti "",\
     "pipe-dist-1-6.dat" u 2:6 w p  pt 16 ps 1.5 color dandelion      every (100)/(1*6) ti "$n=6$ first order",\
     "pipe-dist-2-3.dat" u 2:6 w l  lw 2  lt 3   color periwinkle     ti "",\
     "pipe-dist-2-3.dat" u 2:6 w p  pt 17 ps 1.5 color periwinkle     every (100)/(2*2) ti "$n=2$ second order"

set output "sigmar.pdf"
set xrange [a:b]
set yrange [*:*]
set key top left
set xlabel "$r$ [mm]"
set ylabel "$\sigma_r(r)$ [mm]"
plot sigmar1(x)                    w l  lw 3  lt 1   color emerald      ti "Analytical solution $\sigma_r(r)$",\
     "pipe-dist-1-3.dat" u 2:7 w l  lw 1  lt 2   color bittersweet  ti "",\
     "pipe-dist-1-3.dat" u 2:7 w p  pt 18 ps 1.5 color bittersweet  every (100)/(1*3) ti "$n=3$ first order",\
     "pipe-dist-1-6.dat" u 2:7 w l  lw 1  lt 2   color orange       ti "",\
     "pipe-dist-1-6.dat" u 2:7 w p  pt 16 ps 1.5 color orange       every (100)/(1*6) ti "$n=6$ first order",\
     "pipe-dist-2-3.dat" u 2:7 w l  lw 2  lt 3   color royalblue    ti "",\
     "pipe-dist-2-3.dat" u 2:7 w p  pt 17 ps 1.5 color royalblue    every (100)/(2*2) ti "$n=2$ second order"