{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Minimum tensor distance from a transverse isotropic and an isotropic one\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook we explore the minimum distance between stiffness/compliance of different symmetry clases; namely, transverse isotropic and isotropic. The former is defined by 5 elstic moduli, while the latter just 2. Stating this in an intuitive way: for a given transverse isotropic material we want to find an _isotropic_ equivalent one. We want to **isotropize** the material XD.\n", "\n", "Tha main idea is to find the minimum of\n", "\n", "$$\\min d(C_\\text{iso}, C_\\text{trans}) \\enspace ,$$\n", "\n", "being $d(A, B)$ a distance function between $A$ and $B$. The distance functions used are:\n", "1. Frobenius norm: $d_F(A,B) = \\Vert A - B\\Vert$\n", "2. Log-Euclidean norm: $d_L(A,B) = \\Vert \\log{A} - \\log{B}\\Vert$\n", "3. Riemman norm: $d_R(A,B) = \\Vert \\log{A^{-1/2}BA^{1/2}}\\Vert$\n", "\n", "where $\\Vert M\\Vert = [\\text{tr}{M^T M}]^{1/2}$. The last two are of interest since they are invariant under inversion, i.e., $d_{F/L}(A^{-1}, B^{-1}) = d_{F/L}(A, B)$, a property that is of interest for the physical nature of the tensors." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy\n", "from scipy.linalg import sqrtm, fractional_matrix_power\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from numba import jit" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "Distance functions for 4th order tensors, equivalently it can be\n", "used for (6,6) positive definite symmetric matrices.\n", "\"\"\"\n", "\n", "\n", "def voigt_to_mandel(A):\n", " r\"\"\"Convert from Voigt to Mandel notation\n", " \n", " In Voigt notation the tensor is represented as\n", " \n", " .. math::\n", " C = \\begin{pmatrix}\n", " C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16}\\\\\n", " C_{12} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26}\\\\\n", " C_{13} & C_{23} & C_{33} & C_{34} & C_{35} & C_{36}\\\\\n", " C_{14} & C_{24} & C_{34} & C_{44} & C_{45} & C_{46}\\\\\n", " C_{15} & C_{25} & C_{35} & C_{45} & C_{55} & C_{56}\\\\\n", " C_{16} & C_{26} & C_{36} & C_{46} & C_{56} & C_{66}\n", " \\end{pmatrix}\n", " \n", " and, ind Mandel notation it is represented as\n", " \n", " .. math::\n", " \\hat{C} = \\begin{pmatrix}\n", " C_{11} & C_{12} & C_{13} & \\sqrt{2}C_{14} & \\sqrt{2}C_{15} & \\sqrt{2}C_{16}\\\\\n", " C_{12} & C_{22} & C_{23} & \\sqrt{2}C_{24} & \\sqrt{2}C_{25} & \\sqrt{2}C_{26}\\\\\n", " C_{13} & C_{23} & C_{33} & C_{34} & C_{35} & C_{36}\\\\\n", " \\sqrt{2}C_{14} & \\sqrt{2}C_{24} & \\sqrt{2}C_{34} & 2C_{44} & 2C_{45} & 2C_{46}\\\\\n", " \\sqrt{2}C_{15} & \\sqrt{2}C_{25} & \\sqrt{2}C_{35} & 2C_{45} & 2C_{55} & 2C_{56}\\\\\n", " \\sqrt{2}C_{16} & \\sqrt{2}C_{26} & \\sqrt{2}C_{36} & 2C_{46} & 2C_{56} & 2C_{66}\n", " \\end{pmatrix}\n", " \n", " Parameters\n", " ----------\n", " A : (6,6) array_like\n", " Positive-semidefinite matrix representing stiffness/compliance\n", " tensors.\n", " \n", " References\n", " ----------\n", " .. [1] Voigt notation. (2015, January 9). In Wikipedia, The Free\n", " Encyclopedia. Retrieved May 12, 2015, from \n", " http://en.wikipedia.org/wiki/Voigt_notation\n", " \n", " .. [2] M. Mehrabadi, S.C Cowin. \"Eigentensor of linear anisotropic\n", " elastic materials\". Q.J. Mech. Appl. math, 43 (1990), 15-41.\n", " \n", " .. [3] P. Helnwein (2001). Some Remarks on the Compressed Matrix\n", " Representation of Symmetric Second-Order and Fourth-Order\n", " Tensors. Computer Methods in Applied Mechanics and Engineering,\n", " 190(22-23):2753-2770\n", " \n", " Examples\n", " --------\n", " Let us test with a matrix of ones.\n", " \n", " >>> A = np.ones((6,6))\n", " >>> print np.round(voigt_to_mandel(A), 4)\n", " [[ 1. 1. 1. 1.4142 1.4142 1.4142]\n", " [ 1. 1. 1. 1.4142 1.4142 1.4142]\n", " [ 1. 1. 1. 1.4142 1.4142 1.4142]\n", " [ 1.4142 1.4142 1.4142 2. 2. 2. ]\n", " [ 1.4142 1.4142 1.4142 2. 2. 2. ]\n", " [ 1.4142 1.4142 1.4142 2. 2. 2. ]]\n", " \n", " \n", " \"\"\"\n", " B = A.copy()\n", " B[0:3,3:6] = np.sqrt(2)*B[0:3,3:6]\n", " B[3:6,0:3] = np.sqrt(2)*B[3:6,0:3]\n", " B[3:6,3:6] = 2*B[3:6,3:6]\n", " \n", " return B\n", "\n", "def mat_norm(A):\n", " r\"\"\"Compute the norm of a matrix\n", " \n", " The norm is given by\n", " .. math::\n", " \\Vert A\\Vert = [\\text{tr}{A^T A}]^{1/2}\n", " \n", " Parameters\n", " ----------\n", " A : (n,n) ndarray\n", " Real matrix.\n", "\n", " Returns\n", " -------\n", " norm : float, nonnegative\n", " Norm of the matrix.\n", " \n", " Examples\n", " --------\n", " The following matrix has as Frobenius norm $\\sqrt{10}$.\n", " >>> A = np.array([[2,-1],\n", " ... [-1,2]])\n", " >>> print np.round(mat_norm(A), 6)\n", " 3.162278\n", " \n", " >>> A = np.array([\n", " ... [-3, 5, 7],\n", " ... [2, 6, 4],\n", " ... [0, 2, 8]])\n", " >>> print np.round(mat_norm(A), 6)\n", " 14.387495\n", " \n", " \"\"\"\n", " norm = np.sqrt(np.trace(np.dot(A.T, A)))\n", " return norm\n", "\n", "\n", "def elastic_tensor_dist(A, B, dist=\"riemman\"):\n", " r\"\"\"Compute the distant function for the tensor `A` and `B`\n", "\n", " The distance functions are\n", " \n", " .. math::\n", " \\begin{align}\n", " &d_F(A,B) = \\Vert A - B\\Vert\\\\\n", " &d_L(A,B) = \\Vert \\log{A} - \\log{B}\\Vert\\\\\n", " &d_R(A,B) = \\Vert \\log{A^{-1/2}BA^{1/2}}\\Vert\n", " \\end{align}\n", "\n", " where :math:`\\Vert M\\Vert = [\\text{tr}{M^T M}]^{1/2}`. \n", "\n", " References\n", " ----------\n", " .. [1] Norris, Andrew. \"The isotropic material closest to a given\n", " anisotropic material.\" Journal of Mechanics of Materials and\n", " Structures 1.2 (2006): 223-238.\n", "\n", " .. [2] Moakher, Maher, and Andrew N. Norris. \"The closest elastic\n", " tensor of arbitrary symmetry to an elasticity tensor of lower\n", " symmetry.\" Journal of Elasticity 85.3 (2006): 215-263.\n", " \n", " \n", " Examples\n", " --------\n", " >>> import homogenization as homog\n", " >>> C_Ag = homog.cubic_stiff(43.75, 46.1, 0.43)\n", " >>> C_W = homog.cubic_stiff(407.43, 160.8, 0.28)\n", " >>> print np.round(elastic_tensor_dist(C_Ag, C_W), 6)\n", " 14.597203\n", " \n", " >>> print np.round(elastic_tensor_dist(C_Ag, C_W, 'log'), 6)\n", " 4.102959\n", " \n", " >>> print np.round(elastic_tensor_dist(C_Ag, C_W, 'frob'), 6)\n", " 836.529651\n", "\n", " \"\"\"\n", " A = voigt_to_mandel(A)\n", " B = voigt_to_mandel(B)\n", " \n", " if dist==\"frob\":\n", " C = A - B\n", "\n", " if dist==\"riemman\":\n", " C = np.dot(B, sqrtm(A))\n", " C = np.dot(fractional_matrix_power(A, -0.5), C)\n", " C = scipy.linalg.logm(C)\n", "\n", " if dist==\"log\":\n", " C = scipy.linalg.logm(A) - scipy.linalg.logm(B) \n", " \n", " return mat_norm(C)\n", "\n", "\n", "def trans2iso(props, tensor_type=\"compliance\"):\n", " \"\"\"Return the closest isotropic tensor to an transverse isotropic one\n", "\n", " Parameters\n", " ----------\n", " props : list, [Ep, Et, nu_p, nu_tp, Gt]\n", " Material properties:\n", " \n", " - Ep : float, >0\n", " Young modulus in perpendicular to the axis of symmetry, i.e.,\n", " in the plane of isotropy.\n", " - Et : float, >0\n", " Young modulus in the direction of the axis of symmetry.\n", " - nu_p : float, (-1, 1)\n", " Poisson ratio in the plane of isotropy.\n", " - nu_tp : float, (-1, 1)\n", " Poisson ratio characterizing transverse cotnraction in the plane\n", " of isotropy when tension is applied normal to the plane.\n", " - Gt : float, >0\n", " Shear modulus for the plane perpendicular to the symmetry plane.\n", " tensor_type : string (optional)\n", " String showing the type of tensor used for the minimization:\n", " `stiffness` or `compliance`.\n", " \n", " Returns\n", " -------\n", " E : float, >0\n", " Young modulus for the istropic solid.\n", " nu : float, (-0.5, 1)\n", " Poisson ratio for the istropic solid.\n", " \n", " Examples\n", " --------\n", " Let us consider a transverse isotropic solid given by\n", " `Ep = 9.7e9`, ``Et = 153.75e9``, ``nu_p = 0.672, ``nu_tp = 0.344``,\n", " ``Gt = 5.79e9``.\n", " \n", " >>> Ep = 9.7e9\n", " >>> Et = 153.75e9\n", " >>> nu_p = 0.672\n", " >>> nu_tp = 0.344\n", " >>> Gt = 5.79e9\n", " >>> props = [Ep, Et, nu_p, nu_tp, Gt]\n", "\n", " Minimizing the Frobenius distance for the compliance tensor we should\n", " ``E = 1.225269e10`` and ``nu = 0.3682``.\n", "\n", " >>> E, nu = trans2iso(props)\n", " >>> np.abs(E - 1.225269e10)/1.225269e10 < 0.01\n", " True\n", " >>> np.abs(nu - 0.3682)/0.3682 < 0.01\n", " True\n", "\n", " In the case of stiffnes tensor, the values are ``E = 5.136497e10`` and\n", " ``nu = 0.2112``.\n", " \n", " >>> E, nu = trans2iso(props, tensor_type=\"stiffness\")\n", " >>> np.abs(E - 5.136497e10)/5.136497e10 < 0.01\n", " True\n", " >>> np.abs(nu - 0.2112)/0.2112 < 0.01\n", " True\n", " \n", " \"\"\"\n", " Ep, Et, nu_p, nu_tp, Gt = props\n", " if tensor_type==\"compliance\":\n", " E = -21*Et*Ep*Gt/(((8*nu_tp - 3)*Ep - 10*Et)*Gt - 4*Et*Ep)\n", " nu = -(((10*nu_tp - 2)*Ep + (7*nu_p - 2)*Et)*Gt + 2*Et*Ep)/(\n", " ((8*nu_tp - 3)*Ep - 10*Et)*Gt - 4*Et*Ep)\n", " if tensor_type==\"stiffness\":\n", " E = -((((96*nu_tp**3 + 48*nu_tp**2)*nu_p + 96*nu_tp**3 + 48*nu_tp**2)\n", " *Et*Ep**2 + ((-24*nu_tp**2 + 48*nu_tp + 24)*nu_p**2 + 24*nu_tp**2-\n", " 48*nu_tp-24)*Et**2*Ep + (-12*nu_p**3 + 12*nu_p**2 + 12*nu_p-12)*Et**3)\n", " *Gt + (120*nu_tp**3 + 60*nu_tp**2)*Et*Ep**3 + ((34*nu_tp**2 + 76*nu_tp\n", " + 22)*nu_p + 94*nu_tp**2-44*nu_tp-38)*Et**2*Ep**2 +\n", " ((16*nu_tp + 5)*nu_p**2 + 30*nu_p-16*nu_tp-35)*Et**3*Ep +\n", " (-8*nu_p**3 + 8*nu_p**2 + 8*nu_p-8)*Et**4)/(((48*nu_tp**4*nu_p\n", " + 48*nu_tp**4)*Ep**2 + (48*nu_tp**2*nu_p**2-48*nu_tp**2)*Et*Ep\n", " + (12*nu_p**3-12*nu_p**2-12*nu_p + 12)*Et**2)*Gt + 60*nu_tp**4*Ep**3 +\n", " ((-144*nu_tp**3-36*nu_tp**2)*nu_p-144*nu_tp**3-156*nu_tp**2)*Et*Ep**2\n", " + ((60*nu_tp**2-72*nu_tp-33)*nu_p**2-30*nu_p-60*nu_tp**2 + 72*nu_tp\n", " + 63)*Et**2*Ep + (30*nu_p**3-30*nu_p**2-30*nu_p + 30)*Et**3)\n", "\n", " nu = -(((8*nu_tp**2*nu_p + 8*nu_tp**2)*Ep + (4*nu_p**2-4)*Et)*Gt\n", " + 10*nu_tp**2*Ep**2 + ((20*nu_tp + 11)*nu_p + 20*nu_tp + 1)*Et*Ep\n", " + (1-nu_p**2)*Et**2)/(((8*nu_tp**2*nu_p + 8*nu_tp**2)*Ep\n", " + (4*nu_p**2-4)*Et)*Gt + 10*nu_tp**2*Ep**2\n", " + ((-24*nu_tp-11)*nu_p-24*nu_tp-21)*Et*Ep + (10*nu_p**2-10)*Et**2)\n", "\n", " return E, nu" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def iso_stiff(E, nu):\n", " r\"\"\"Form the stiffness tensor in Voigt notation.\n", "\n", " Parameters\n", " ----------\n", " E : float\n", " Young modulus.\n", " nu : float\n", " Poisson's ratio.\n", "\n", " Returns\n", " -------\n", " stiff : (6,6) ndarray\n", " Stiffness tensor in Voigt notation.\n", "\n", " Examples\n", " --------\n", " Let us consider `E=200` GPa and `nu=0.285`, like in steel, the\n", " units are GPa to have small numbers.\n", "\n", " >>> from homogenization import iso_stiff\n", " >>> E = 200\n", " >>> nu = 0.285\n", " >>> C = iso_stiff(E, nu)\n", " >>> print np.round(C,4)\n", " [[ 258.8001 103.1581 103.1581 0. 0. 0. ]\n", " [ 103.1581 258.8001 103.1581 0. 0. 0. ]\n", " [ 103.1581 103.1581 258.8001 0. 0. 0. ]\n", " [ 0. 0. 0. 77.821 0. 0. ]\n", " [ 0. 0. 0. 0. 77.821 0. ]\n", " [ 0. 0. 0. 0. 0. 77.821 ]]\n", "\n", " See Also\n", " --------\n", " iso_compl\n", "\n", " \"\"\"\n", "\n", " mu = E / (2.0 * (1.0 + nu))\n", " lam = E * nu / (1.0 + nu) / (1.0 - 2.0 * nu)\n", "\n", " stiff = np.zeros((6, 6))\n", "\n", " stiff[0, 0] = 2.0 * mu + lam\n", " stiff[1, 1] = 2.0 * mu + lam\n", " stiff[2, 2] = 2.0 * mu + lam\n", " stiff[0, 1] = lam\n", " stiff[0, 2] = lam\n", " stiff[1, 0] = lam\n", " stiff[1, 2] = lam\n", " stiff[2, 0] = lam\n", " stiff[2, 1] = lam\n", " stiff[3, 3] = mu\n", " stiff[4, 4] = mu\n", " stiff[5, 5] = mu\n", "\n", " return stiff" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def transiso_stiff(Ep, Et, nu_p, nu_tp, Gt):\n", " \"\"\"Form the stiffness tensor for transverse isotropic cases.\n", "\n", " Parameters\n", " ----------\n", " Ep : float\n", " Young modulus in the plane of isotropy.\n", " Et : float\n", " Young modulus perpendicular to the plane of isotropy.\n", " nu_p : float\n", " Poisson's ratio in the plane of isotropy.\n", " nu_tp : float\n", " Poisson's ration perpendicular to the plane of isotropy.\n", " G_t : float\n", " Shear modulus for the plane perpendicular to the symmetry\n", " plane.\n", "\n", " Returns\n", " -------\n", " stiff : (6,6) ndarray\n", " Compliance tensor in Voigt notation.\n", "\n", "\n", " The axis of symmetry in this case is x3 (z).\n", "\n", " Examples\n", " --------\n", " Let us consider a homogeneous material first `E=70` GPa and\n", " `nu=0.33`, like in aluminum, the units are GPa to have small\n", " numbers.\n", "\n", " >>> from homogenization import transiso_stiff\n", " >>> Ep = 70\n", " >>> Et = 70\n", " >>> nu_p = 0.33\n", " >>> nu_tp = 0.33\n", " >>> Gt = Ep/(2*(1 + nu_p))\n", " >>> C = transiso_stiff(Ep, Et, nu_p, nu_tp, Gt)\n", " >>> print np.round(C,4)\n", " [[ 103.7152 51.0836 51.0836 0. 0. 0. ]\n", " [ 51.0836 103.7152 51.0836 0. 0. 0. ]\n", " [ 51.0836 51.0836 103.7152 0. 0. 0. ]\n", " [ 0. 0. 0. 26.3158 0. 0. ]\n", " [ 0. 0. 0. 0. 26.3158 0. ]\n", " [ 0. 0. 0. 0. 0. 26.3158]]\n", "\n", " Now, let us consider apatite that has as parameters\n", " `Ep = 133.50` GPa, `Et = 91.63` GPa, `nu_p = -0.133`,\n", " `nu_tp = 0.366` and `Gt = 66.3` GPa.\n", "\n", " >>> from homogenization import transiso_stiff\n", " >>> Ep = 133.498\n", " >>> Et = 91.6269\n", " >>> nu_p = -0.1326\n", " >>> nu_tp = 0.3665\n", " >>> Gt = 66.3\n", " >>> C = transiso_stiff(Ep, Et, nu_p, nu_tp, Gt)\n", " >>> print np.round(C,4)\n", " [[ 167.0093 13.1033 66.0113 0. 0. 0. ]\n", " [ 13.1033 167.0093 66.0113 0. 0. 0. ]\n", " [ 66.0113 66.0113 140.0132 0. 0. 0. ]\n", " [ 0. 0. 0. 66.3 0. 0. ]\n", " [ 0. 0. 0. 0. 66.3 0. ]\n", " [ 0. 0. 0. 0. 0. 76.953 ]]\n", "\n", "\n", "\n", " \"\"\"\n", " C11 = (Ep**2 * nu_tp**2 - Ep * Et) / (Ep * (2 * nu_p + 2) * nu_tp**2 +\n", " Et * (nu_p**2 - 1.0))\n", " C12 = (-Ep**2 * nu_tp**2 - Ep * Et * nu_p) / \\\n", " (Ep * (2 * nu_p + 2) * nu_tp**2 + Et * (nu_p**2 - 1.0))\n", " C13 = -(Ep * Et * nu_tp) / (2 * Ep * nu_tp**2 + Et * (nu_p - 1.0))\n", " C33 = (Et**2 * (nu_p - 1.0)) / (2 * Ep * nu_tp**2 + Et * (nu_p - 1.0))\n", "\n", " stiff = np.array([\n", " [C11, C12, C13, 0.0, 0.0, 0.0],\n", " [C12, C11, C13, 0.0, 0.0, 0.0],\n", " [C13, C13, C33, 0.0, 0.0, 0.0],\n", " [0.0, 0.0, 0.0, Gt, 0.0, 0.0],\n", " [0.0, 0.0, 0.0, 0.0, Gt, 0.0],\n", " [0.0, 0.0, 0.0, 0.0, 0.0, Ep / (2 * nu_p + 2)]\n", " ])\n", "\n", " return stiff" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def compute_dist(props, E_vec, nu_vec):\n", " Ep, Et, nu_p, nu_tp, Gt = props\n", "\n", " Ctr = transiso_stiff(Ep,Et,nu_p,nu_tp,Gt)\n", "\n", " N = len(E_vec)\n", " dL = np.zeros((N, N))\n", " dR = np.zeros((N, N))\n", " dF = np.zeros((N, N))\n", " for row, E in enumerate(E_vec):\n", " for col, nu in enumerate(nu_vec):\n", " Ciso = iso_stiff(E,nu) \n", " dL[row,col] = elastic_tensor_dist(Ciso, Ctr, dist=\"log\")\n", " dR[row,col] = elastic_tensor_dist(Ciso, Ctr, dist=\"riemman\")\n", " dF[row,col] = elastic_tensor_dist(Ciso, Ctr, dist=\"frob\")\n", " \n", " return dL, dR, dF" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def plot_setup(nu_vec, E_vec, props, dFun, levels):\n", " \n", " N = len(E_vec)\n", " # Contour plot\n", " ax = fig.add_subplot(1, 2, 1)\n", " plt.contourf(nu_vec, E_vec, dFun, levels, cmap='YlGn', extend=\"max\")\n", " plt.colorbar()\n", " plt.contour(nu_vec, E_vec, dFun, levels, colors=\"k\")\n", " Em, num = trans2iso(props)\n", " plt.plot(num, Em, 'ob', label=\"Compliance based\")\n", " Em2, num2 = trans2iso(props, tensor_type=\"stiffness\")\n", " plt.plot(num2, Em2, 'sr', label=\"Stiffness based\")\n", " argmin = np.argmin(dFun)\n", " plt.plot(nu_vec[argmin%N], E_vec[argmin//N], '*k', label=\"Minimum\")\n", " plt.legend(loc=3)\n", " plt.xlabel(r\"$\\nu$\")\n", " plt.ylabel(r\"$E$\")\n", " \n", " \n", " ax = fig.add_subplot(1, 2, 2, projection='3d')\n", " surf = ax.plot_surface(nu_vec, E_vec, dFun, rstride=1, cstride=1, cmap='YlGn',\n", " linewidth=0, antialiased=False)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/nguarinz/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:14: ComplexWarning: Casting complex values to real discards the imaginary part\n", " \n" ] } ], "source": [ "# Carbon-epoxy\n", "Ep = 9.7\n", "Et = 153.75\n", "nu_p = 0.672\n", "nu_tp = 0.344\n", "Gt = 5.79\n", "\n", "props = [Ep, Et, nu_p, nu_tp, Gt]\n", "\n", "tol = 1e-3\n", "npts = 50\n", "E_vec = np.linspace(tol, npts, npts)\n", "nu_vec = np.linspace(-1 + tol, 0.5 - tol, npts)\n", "\n", "dL, dR, dF = compute_dist(props, E_vec, nu_vec)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('