Browaeys and Chevrot [work in progress]

examples/decomp.ipynb

@@ -0,0 +1,228 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Decomposing elastic constants\n",
+    "\n",
+    "This example provides an illustration of the implementation of the decomposition of elastic constants following the approach of Browaeys and Chevrot (2004). Stishovite elastic constants are used as for the example. The implementation closely follows that in MSAT (Walker and Wookey, 2012) and decomposition follows three stages. First, the input elastic constants are rotated such that the axes lie along high symmetry directions. Second, the matrix is decomposed into contributions for each symmetry. Third, the magnitude of each symmetry component is calculated. \n",
+    "\n",
+    "### References\n",
+    "\n",
+    "J. T. Browaeys and S. Chevrot. Decomposition of the elastic tensor and geophysical applications. Geophysical Journal International, 159:667 – 678, 2004.\n",
+    "\n",
+    "A. M. Walker and J. Wookey. MSAT – a new toolkit for the analysis of elastic and seismic anisotropy. Computers and Geosciences, 49:81 – 90, 2012."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# Module imports etc\n",
+    "import sys\n",
+    "sys.path.append('..')\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# We need to import decompose and rotate\n",
+    "import pytasa.decompose\n",
+    "import pytasa.rotate\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "## Setup the input\n",
+    "\n",
+    "Here we use the elastic constants of stishovite, but we rotate them in an arbitary way so that the symmetry is not obvious."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Stishovite elastic constants:\n",
+      "[[ 453.  211.  203.    0.    0.    0.]\n",
+      " [ 211.  453.  203.    0.    0.    0.]\n",
+      " [ 203.  203.  776.    0.    0.    0.]\n",
+      " [   0.    0.    0.  252.    0.    0.]\n",
+      " [   0.    0.    0.    0.  252.    0.]\n",
+      " [   0.    0.    0.    0.    0.  302.]]\n",
+      "Rotated stishovite elastic constants:\n",
+      "[[ 635.    87.3  150.9   77.6    1.6  -20.6]\n",
+      " [  87.3  679.3  200.6   31.8  -13.4  -14.9]\n",
+      " [ 150.9  200.6  724.2   35.2  -24.2   10.6]\n",
+      " [  77.6   31.8   35.2  250.4   13.9   -8.6]\n",
+      " [   1.6  -13.4  -24.2   13.9  213.2   58.3]\n",
+      " [ -20.6  -14.9   10.6   -8.6   58.3  164.1]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "stishovite_cij = np.array([[453.0, 211.0, 203.0,   0.0,   0.0,   0.0],\n",
+    "                           [211.0, 453.0, 203.0,   0.0,   0.0,   0.0],\n",
+    "                           [203.0, 203.0, 776.0,   0.0,   0.0,   0.0],\n",
+    "                           [  0.0,   0.0,   0.0, 252.0,   0.0,   0.0],\n",
+    "                           [  0.0,   0.0,   0.0,   0.0, 252.0,   0.0],\n",
+    "                           [  0.0,   0.0,   0.0,   0.0,   0.0, 302.0]])\n",
+    "\n",
+    "print(\"Stishovite elastic constants:\")\n",
+    "print(np.array_str(stishovite_cij, precision=1, suppress_small=True))\n",
+    "\n",
+    "# Rotate around three axes\n",
+    "stishovite_cij = pytasa.rotate.rot_3(stishovite_cij, 20, 30, 40)\n",
+    "stishovite_rho = 4290.0\n",
+    "\n",
+    "print(\"Rotated stishovite elastic constants:\")\n",
+    "print(np.array_str(stishovite_cij, precision=1, suppress_small=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Find the high symmetry oritentation\n",
+    "This works by computing the Voigt stiffness tensor and dilational stiffness tensor from the input elastic constants tensor (represented as a matrix using Voigt notation). For high symmetry materials the eigenvectors allign. The number of distinct eigenvalues provides information about the symmetry. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "FIXME: need to issue a warning\n",
+      "Elastic constants matrix alligned with symmetry directions:\n",
+      "[[ 592.2   71.8  203.    -0.    -0.   -76.3]\n",
+      " [  71.8  592.2  203.     0.    -0.    76.3]\n",
+      " [ 203.   203.   776.     0.    -0.    -0. ]\n",
+      " [  -0.     0.     0.   252.     0.    -0. ]\n",
+      " [  -0.    -0.     0.    -0.   252.    -0. ]\n",
+      " [ -76.3   76.3    0.    -0.    -0.   162.8]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "cij_on_axes, rotation_matrix = pytasa.decompose.axes(stishovite_cij)\n",
+    "print(\"Elastic constants matrix alligned with symmetry directions:\")\n",
+    "print(np.array_str(cij_on_axes, precision=1, suppress_small=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# FIXME! - this is clearly wrong!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Decompose and calculate norms\n",
+    "Project the rotated matrix into isotropic, hexagonal, tetragonal, orthorhombic, monoclinic and triclinic components, then calculate the norms."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "decomposed = pytasa.decompose.decompose(cij_on_axes)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# FIXME: this is a bit rubbish\n",
+    "norms = pytasa.decompose.norms(cij_on_axes, decomposed[0], decomposed[1], decomposed[2], \n",
+    "                               decomposed[3], decomposed[4], decomposed[5])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "ElasticNorms(isotropic=0.76694217270623366, hexagonal=0.058050715699155098, tetragonal=0.027521562819149636, orthorhombic=0.0, monoclinic=0.14748554877546138, triclinic=1.1102230246251565e-16)\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(norms)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.4.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}