{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Concepts for Diagonalization of strain\n",
    "\n",
    "This jupyter notebook is a worksheet that demostrates how coordinate transformation of strain works in the shear elastic modulus part works this program. You could clone this repo and change the `key` below and rerun it to see how the shear modulus solving mechanism is working for each of the elastic modulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding the transformation matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a shear elastic modulus $c_{ij}$ below as `key`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from cij.util import c_\n",
    "\n",
    "key = c_(46)\n",
    "\n",
    "if not key.is_shear: raise RuntimeError(f\"The strain c_{key} is not a shear modulus.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the corresponding fictitous strain $e$ for $c_{ij}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 1. 0.]\n",
      " [1. 0. 1.]\n",
      " [0. 1. 0.]]\n"
     ]
    }
   ],
   "source": [
    "e = numpy.zeros((3, 3))\n",
    "\n",
    "e[key.i[0] - 1, key.i[1] - 1] = 1\n",
    "e[key.i[1] - 1, key.i[0] - 1] = 1\n",
    "\n",
    "e[key.j[0] - 1, key.j[1] - 1] = 1\n",
    "e[key.j[1] - 1, key.j[0] - 1] = 1\n",
    "\n",
    "print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transformation matrix $T$ is the column-wise eigenvectors of the strain matrix $e$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.00000000e-01  7.07106781e-01  5.00000000e-01]\n",
      " [-7.07106781e-01  9.02056208e-17  7.07106781e-01]\n",
      " [ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]]\n"
     ]
    }
   ],
   "source": [
    "evals, T = numpy.linalg.eig(e)\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $e'$ in the rotated coordinate system is the matrix with eigenvalues of $e$ in its diagonals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.41421356e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  9.77950360e-17  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  1.41421356e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(numpy.diag(evals))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can retrive the original strain with\n",
    "\n",
    "$$\n",
    "e' = T^{-1} e T\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e = \n",
      "[[-1.11022302e-16  1.00000000e+00 -1.11022302e-16]\n",
      " [ 1.00000000e+00 -1.11022302e-16  1.00000000e+00]\n",
      " [-1.11022302e-16  1.00000000e+00  1.11022302e-16]]\n",
      "e = T e' T^{-1}: True\n"
     ]
    }
   ],
   "source": [
    "print(\"e = \")\n",
    "print(T @ numpy.diag(evals) @ T.T)\n",
    "\n",
    "print(\"e = T e' T^{-1}:\", numpy.allclose(e, T @ numpy.diag(evals) @ T.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or backwards\n",
    "\n",
    "$$\n",
    "e = T e' T^{-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e' = \n",
      "[[-1.41421356e+00 -1.11022302e-16  0.00000000e+00]\n",
      " [-6.68036251e-17  4.93038066e-32  2.47214867e-16]\n",
      " [ 0.00000000e+00  3.33066907e-16  1.41421356e+00]]\n",
      "e' = T^{-1} e T: True\n"
     ]
    }
   ],
   "source": [
    "print(\"e' = \")\n",
    "print(T.T @ e @ T)\n",
    "\n",
    "print(\"e' = T^{-1} e T:\", numpy.allclose(numpy.diag(evals), T.T @ e @ T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing with symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "\n",
    "def round_expr(expr, num_digits=4):\n",
    "    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(sympy.Number)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The strain energy invariance equation\n",
    "\n",
    "The strain energy $E = \\sum \\frac{1}{2} c_{ijkl}e_{ij}e_{kl}$ should be same in the original and rotated coordinate system. Here we show the energy in both coordinate systems with the help of `sympy`.  The function used to calculate the strain energy is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cij.core.phonon_contribution import calculate_fictitious_strain_energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **left hand side** of the strain energy invariance equation, which is under the **original coordinate system**, the expression for the ficticious strain $e$ defined above is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.0 c_{44} + 4.0 c_{46} + 2.0 c_{66}$"
      ],
      "text/plain": [
       "2.0*c_44 + 4.0*c_46 + 2.0*c_66"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def resolve_elastic_modulus(key):\n",
    "    return sympy.symbols(\"c_%d%d\" % key.v)\n",
    "    \n",
    "calculate_fictitious_strain_energy(\n",
    "    e,\n",
    "    resolve_elastic_modulus\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and on the **right hand side**, which is under the **rotated coordinate system**, the expression for the rotated ficticious strain $e'$ defined above is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.0 c_{1'1'} - 2.0 c_{1'3'} + 1.0 c_{3'3'}$"
      ],
      "text/plain": [
       "1.0*c_1'1' - 2.0*c_1'3' + 1.0*c_3'3'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def resolve_elastic_modulus_rotated(key):\n",
    "    return sympy.symbols(\"c_%d'%d'\" % key.v)\n",
    "\n",
    "e_prime = T.T @ e @ T\n",
    "\n",
    "calculate_fictitious_strain_energy(\n",
    "    numpy.diag(numpy.diag(e_prime)),\n",
    "    resolve_elastic_modulus_rotated,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The strain used for calculating strain-Grüneisen parameter from volume-Grüneisen parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.25 e_{11} + 0.5 e_{22} + 0.25 e_{33} & 0.3536 e_{11} - 0.3536 e_{33} & 0.25 e_{11} - 0.5 e_{22} + 0.25 e_{33}\\\\0.3536 e_{11} - 0.3536 e_{33} & 0.5 e_{11} + 0.5 e_{33} & 0.3536 e_{11} - 0.3536 e_{33}\\\\0.25 e_{11} - 0.5 e_{22} + 0.25 e_{33} & 0.3536 e_{11} - 0.3536 e_{33} & 0.25 e_{11} + 0.5 e_{22} + 0.25 e_{33}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0.25*e_11 + 0.5*e_22 + 0.25*e_33, 0.3536*e_11 - 0.3536*e_33, 0.25*e_11 - 0.5*e_22 + 0.25*e_33],\n",
       "[       0.3536*e_11 - 0.3536*e_33,       0.5*e_11 + 0.5*e_33,        0.3536*e_11 - 0.3536*e_33],\n",
       "[0.25*e_11 - 0.5*e_22 + 0.25*e_33, 0.3536*e_11 - 0.3536*e_33, 0.25*e_11 + 0.5*e_22 + 0.25*e_33]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = sympy.symbols(\"e_11\")\n",
    "b = sympy.symbols(\"e_22\")\n",
    "c = sympy.symbols(\"e_33\")\n",
    "\n",
    "m = sympy.Matrix([\n",
    "    [a, 0, 0],\n",
    "    [0, b, 0],\n",
    "    [0, 0, c]\n",
    "])\n",
    "\n",
    "round_expr(T.T @ m @ T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.25 α + 0.5 β + 0.25 γ + 1.0 & 0.3536 α - 0.3536 γ & 0.25 α - 0.5 β + 0.25 γ\\\\0.3536 α - 0.3536 γ & 0.5 α + 0.5 γ + 1.0 & 0.3536 α - 0.3536 γ\\\\0.25 α - 0.5 β + 0.25 γ & 0.3536 α - 0.3536 γ & 0.25 α + 0.5 β + 0.25 γ + 1.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0.25*α + 0.5*β + 0.25*γ + 1.0, 0.3536*α - 0.3536*γ,       0.25*α - 0.5*β + 0.25*γ],\n",
       "[          0.3536*α - 0.3536*γ, 0.5*α + 0.5*γ + 1.0,           0.3536*α - 0.3536*γ],\n",
       "[      0.25*α - 0.5*β + 0.25*γ, 0.3536*α - 0.3536*γ, 0.25*α + 0.5*β + 0.25*γ + 1.0]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = sympy.symbols(\"α\")\n",
    "b = sympy.symbols(\"β\")\n",
    "c = sympy.symbols(\"γ\")\n",
    "\n",
    "m = sympy.Matrix([\n",
    "    [1 + a, 0.   , 0.   ],\n",
    "    [0.   , 1 + b, 0.   ],\n",
    "    [0.   , 0.   , 1 + c]\n",
    "])\n",
    "\n",
    "round_expr(T.T @ m @ T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Order for modulus calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[11(1111),\n",
       " 22(2222),\n",
       " 33(3333),\n",
       " 12(1122),\n",
       " 13(1133),\n",
       " 23(2233),\n",
       " 44(2323),\n",
       " 55(1313),\n",
       " 66(1212),\n",
       " 15(1113),\n",
       " 25(2213),\n",
       " 35(3313),\n",
       " 46(2312)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from cij.util.voigt import _cij_sort_key\n",
    "from IPython.display import display\n",
    "\n",
    "sorted([\n",
    "    c_(11), c_(12), c_(13), c_(22), c_(23), c_(33), c_(15), c_(25), c_(35), c_(46), c_(44), c_(55), c_(66)\n",
    "], key=_cij_sort_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}