#. ## Orthotropic free expansion of a cube
#.
#. To illustrate the point of the previous discussion, let us solve the thermal
#. expansion of an unrestrained unitary cube $[0,1~\text{mm}]\times[0,1~\text{mm}]\times[0,1~\text{mm}]$
#. subject to a linear radially-symmetric temperature field
#. $$
#. T(x,y,z) = 30 \text{ºC} + 150 \frac{\text{ºC}}{\text{mm}} \sqrt{x^2+y^2+z^2}
#. $$
#.
#. with a mean thermal expansion coefficient for each of the three directions $x$, $y$ and $z$
#. computed from each of the three columns of the ASME table TE-2, respectively.
#. If the data was consistent, the displacement at any point with the same
#. coordinates $x=y=z$ would be exactly equal.
#$ $ gmsh -3 cube-hex.geo
#$ [...]
#$ $ gmsh -3 cube-tet.geo
#$ [...]
#$ $ feenox cube-orthotropic-expansion.fee 
#$ displacement in x at (1,1,1) =  4.451e-03
#$ displacement in y at (1,1,1) =  4.449e-03
#$ displacement in z at (1,1,1) =  4.437e-03
#$ $ feenox cube-orthotropic-expansion.fee linear tet
#$ displacement in x at (1,1,1) =  4.451e-03
#$ displacement in y at (1,1,1) =  4.447e-03
#$ displacement in z at (1,1,1) =  4.438e-03
#$ $ feenox cube-orthotropic-expansion.fee akima hex
#$ displacement in x at (1,1,1) =  4.451e-03
#$ displacement in y at (1,1,1) =  4.451e-03
#$ displacement in z at (1,1,1) =  4.437e-03
#$ $ feenox cube-orthotropic-expansion.fee splines tet
#$ displacement in x at (1,1,1) =  4.451e-03
#$ displacement in y at (1,1,1) =  4.450e-03
#$ displacement in z at (1,1,1) =  4.438e-03
#$ $ 
DEFAULT_ARGUMENT_VALUE 1 steffen
DEFAULT_ARGUMENT_VALUE 2 hex

PROBLEM mechanical
READ_MESH cube-hex.msh

# aluminum-like linear isotropic material properties
E = 69e3
nu = 0.28

# free expansion
BC left   u=0
BC front  v=0
BC bottom w=0

# reference temperature is 20ºC
T0 = 20
# spatial temperature distribution symmetric wrt x,y & z
T(x,y,z) = 30+150*sqrt(x^2+y^2+z^2)

# read ASME data
FUNCTION A(T') FILE asme-expansion-table.dat COLUMNS 1 2 INTERPOLATION $1
FUNCTION B(T') FILE asme-expansion-table.dat COLUMNS 1 3 INTERPOLATION $1
FUNCTION C(T') FILE asme-expansion-table.dat COLUMNS 1 4 INTERPOLATION $1

# remember that the thermal expansion coefficients have to be
#  1. the mean value between T0 and T
#  2. functions of space, so temperature has to be written as T(x,y,z)

# in the x direction, we use column B directly
alpha_x(x,y,z) = 1e-6*B(T(x,y,z))

# in the y direction, we convert column A to mean
alpha_y(x,y,z) = 1e-6*integral(A(T'), T', T0, T(x,y,z))/(T(x,y,z)-T0)

# in the z direction, we convert column C to mean
alpha_z(x,y,z) = 1e-3*C(T(x,y,z))/(T(x,y,z)-T0)

SOLVE_PROBLEM

WRITE_MESH cube-orthotropic-expansion-$1-$2.vtk T VECTOR u v w
PRINT %.3e "displacement in x at (1,1,1) = " u(1,1,1)
PRINT %.3e "displacement in y at (1,1,1) = " v(1,1,1)
PRINT %.3e "displacement in z at (1,1,1) = " w(1,1,1)
#> ![Warped displacement ($\times 500$) of the cube using ASME's three columns.](cube-orthotropic-expansion.png){width=100%}
#>
#> Differences cannot be seen graphically, but they are there as the terminal mimic illustrates.
#> Yet, they are not as large nor as sensible to meshing and interpolation settings as one would
#> have expected after seeing the plots from the previous section. 
