{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "import numpy as np\n", "import itertools" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "#ndim = 2#3 # number of dimensions\n", "N = 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# tensor operations/products: np.einsum enables index notation, avoiding loops\n", "# e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid\n", "trans2 = lambda A2 : np.einsum('ij... ->ji... ',A2 )\n", "ddot42 = lambda A4,B2: np.einsum('ijkl...,lk... ->ij... ',A4,B2)\n", "ddot44 = lambda A4,B4: np.einsum('ijkl...,lkmn...->ijmn...',A4,B4)\n", "dot11 = lambda A1,B1: np.einsum('ixyz ,ixyz ->xyz ',A1,B1)\n", "dot22 = lambda A2,B2: np.einsum('ij... ,jk... ->ik... ',A2,B2)\n", "dot24 = lambda A2,B4: np.einsum('ij... ,jkmn...->ikmn...',A2,B4)\n", "dot42 = lambda A4,B2: np.einsum('ijkl...,lm... ->ijkm...',A4,B2)\n", "dyad22 = lambda A2,B2: np.einsum('ij... ,kl... ->ijkl...',A2,B2)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# material parameters + function to convert to grid of scalars\n", "phase = .1\n", "param = lambda M0,M1: M0*np.ones([N,N,N])*(1.-phase)+M1*np.ones([N,N,N])*phase\n", "Kv = param(0.833,8.33) # bulk modulus [grid of scalars]\n", "mu = param(0.386,3.86) # shear modulus [grid of scalars]\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# bulk modulus K = lambda - 2*µ/3\n", "# shear modulus µ = lamé's µ\n", "Young = 9*Kv*mu/(3*Kv + mu)\n", "Poisson = (3*Kv-2*mu)/(2*(3*Kv+ mu))\n", "print(\"Young = {}\".format(Young.flatten()[0]))\n", "print(\"Poisson = {}\".format(Poisson.flatten()[0]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def mat(dim):\n", " if dim == 2:\n", " return np.array([[0, 0],\n", " [1, 1],\n", " [0, 1],\n", " [1, 0]])\n", " elif dim == 3:\n", " return np.array([[0, 0],\n", " [1, 1],\n", " [2, 2],\n", " [1, 2],\n", " [0, 2],\n", " [0, 1],\n", " [2, 1],\n", " [2, 0],\n", " [1, 0]])\n", " else:\n", " raise Exception(\"unsupported dimension\")\n", "\n", "def vsize(dim):\n", " return dim * (dim - 1) // 2 + dim\n", "def V4(arr, sym=True):\n", "\n", " dim = arr.shape[0]\n", " if sym:\n", " siz = vsize(dim)\n", " else:\n", " siz = arr.shape[0]**2\n", "\n", " vout = np.zeros([siz, siz])\n", " ids = mat(dim)\n", " for I in range(siz):\n", " i, j = ids[I]\n", " for J in range(siz):\n", " k, l = ids[J]\n", " vout[I, J] = arr [i, j, k, l]\n", " pass\n", " pass\n", " return vout\n", "\n", "def print_eigen_input(arr):\n", " dim = len(arr.shape)\n", " assert dim in {1, 2}, \"only 1d and 2d arrays\"\n", "\n", " if dim == 1:\n", " return \"<<\\n\" + \", \".join((str(i) for i in arr.tolist()))\n", " else:\n", " return \"<<\\n\" + \",\\n\".join(\n", " \", \".join((str(i) for i in row))\n", " for row in arr.tolist())\n", " pass\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def constitutive(F):\n", " ndim = F.shape[0]\n", " # identity tensor [single tensor]\n", " i = np.eye(ndim)\n", " # identity tensors [grid of tensors]\n", " I = i\n", "\n", " I4 = np.einsum('il,jk',i,i)\n", "\n", " I4rt = np.einsum('ik,jl',i,i)\n", "\n", " I4s = (I4+I4rt)/2.\n", "\n", " II = dyad22(I,I)\n", "\n", " C4 = Kv*II+2.*mu*(I4s-1./3.*II)\n", " S = ddot42(C4,.5*(dot22(trans2(F),F)-I))\n", " P = dot22(F,S)\n", " K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt)\n", " return P,K4,C4,S\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "\n", "Fs = [np.array([[1, .1],\n", " [.2, 1]]),\n", " np.array([[1, .1, 0],\n", " [.2, 1, 0],\n", " [0, 0, 1]]),\n", " np.array([[1, .1, .3],\n", " [.2, 1, .4],\n", " [.5, .6, 1]])]\n", "\n", "for F in Fs:\n", " P, K, C, S = constitutive(F)\n", " print(\"F =\\n{}\".format(print_eigen_input(F )))\n", " print(\"P =\\n{}\".format(print_eigen_input(P )))\n", " print(\"K =\\n{}\".format(print_eigen_input(V4(K, sym=False))))\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def illustrate(F):\n", " ndim = F.shape[0]\n", " # identity tensor [single tensor]\n", " i = np.eye(ndim)\n", " # identity tensors [grid of tensors]\n", " I = i\n", "\n", " I4 = np.einsum('il,jk',i,i)\n", "\n", " I4rt = np.einsum('ik,jl',i,i)\n", "\n", " I4s = (I4+I4rt)/2.\n", "\n", " II = dyad22(I,I)\n", "\n", " C4 = Kv*II+2.*mu*(I4s-1./3.*II)\n", " S = ddot42(C4,.5*(dot22(trans2(F),F)-I))\n", " P = dot22(F,S)\n", " def printa(name, var):\n", " print(\"{} =\\n{}\".format(name, V4(var, sym=False)))\n", " C_FT = dot42(C4, trans2(F))\n", " printa('C_FT', C_FT)\n", " F_C_FT = dot24(F, C_FT)\n", " printa('F_C_FT', F_C_FT)\n", " IRT_F_C_FT = ddot44(I4rt, F_C_FT)\n", " printa('IRT_F_C_FT', IRT_F_C_FT)\n", " IRT_F_C_FT_IRT = ddot44(IRT_F_C_FT, I4rt)\n", " printa('IRT_F_C_FT_IRT', IRT_F_C_FT_IRT)\n", " S_I = dot24(S,I4)\n", " printa('S_I', S_I)\n", " my_K4 = S_I + IRT_F_C_FT_IRT\n", " K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt)\n", " print(\"max Error = {}\".format(abs(my_K4 - K4).max()))\n", "\n", "illustrate(Fs[0])" ] }, { "cell_type": "markdown", "metadata": { "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "source": [ "Computation of Projection operator" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def Ghat(ndim, N):\n", " delta = lambda i,j: np.float(i==j) # Dirac delta function\n", " freq = np.arange(-(N-1)/2.,+(N+1)/2.) # coordinate axis -> freq. axis\n", " Ns = [N for _ in range(ndim)]\n", " Ghat4 = np.zeros([ndim,ndim,ndim,ndim]+Ns) # zero initialize\n", " # - compute\n", " for i,j,l,m in itertools.product(range(ndim),repeat=4):\n", " for xyz in itertools.product(range(N), repeat=ndim):\n", " q = np.array([freq[coord] for coord in xyz]) # frequency vector\n", " # zero freq. -> mean\n", " if not q.dot(q) == 0:\n", " if len(xyz) == 2:\n", " Ghat4[i,j,l,m, xyz[0], xyz[1]] = delta(i,m)*q[j]*q[l]/(q.dot(q))\n", " elif len(xyz) == 3:\n", " Ghat4[i,j,l,m, xyz[0], xyz[1], xyz[2]] = delta(i,m)*q[j]*q[l]/(q.dot(q))\n", " else:\n", " print(xyz)\n", " raise Exception(\"'{}' is an illegal dim\".format(len(xyz)))\n", " return Ghat4\n", "\n", "\n", "print(Ghat(2,3).shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def Ghat2(ndim, Nx, Ny, Nz):\n", " shape = Nx, Ny, Nz\n", " Ghat4 = np.zeros([ndim,ndim,ndim,ndim,Nx,Ny,Nz]) # projection operator\n", " x = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # position vectors\n", " q = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # frequency vectors\n", " q2 = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # frequency vectors\n", " delta = lambda i,j: np.float(i==j) # Dirac delta function\n", " # - set \"x\" as position vector of all grid-points [grid of vector-components]\n", " x[0],x[1],x[2] = np.mgrid[:Nx,:Ny,:Nz]\n", " # - convert positions \"x\" to frequencies \"q\" [grid of vector-components]\n", " for i in range(3):\n", " #freq = np.arange(-(shape[i]-1)/2,+(shape[i]+1)/2,dtype='int64')\n", " freq = np.fft.fftfreq(shape[i], 1/shape[i])\n", " freq2 = np.fft.fftfreq(shape[i])\n", " q[i] = freq[x[i]]\n", " q2[i] = freq[x[i]]\n", " # - compute \"Q = ||q||\", and \"norm = 1/Q\" being zero for the mean (Q==0)\n", " # NB: avoid zero division\n", " q = q.astype(np.float)\n", " Q = dot11(q,q)\n", " Z = Q==0\n", " Q[Z] = 1.\n", " norm = 1./Q\n", " norm[Z] = 0.\n", " # - set projection operator [grid of tensors]\n", " for i, j, l, m in itertools.product(range(ndim), repeat=4):\n", " Ghat4[i,j,l,m] = norm*delta(i,m)*q[j]*q[l]\n", " pass\n", " return Ghat4\n", "\n", "def GhatBlackwiseRowMajor(G):\n", " mainshape = G.shape[:4]\n", " flatshape = G.shape[4:]\n", " ni, nj, nk = flatshape\n", " for i in range(ni):\n", " for j in range(nj):\n", " for k in range(nk):\n", " print(\"//({},{},{})\".format(i,j,k))\n", " print(\"retval.emplace_back(std::move(std::make_unique>()));\")\n", " print(\"*retval.back() {};\".format(print_eigen_input(V4(G[:,:,:,:,i,j,k], sym=False))))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "G = Ghat2(2, 5,3,1)\n", "print(\"empty-ratio: {}\".format((G==0).sum()/G.size))\n", "print(\"size: {}:\".format(G.size))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "G = Ghat2(2, 5,3,1)\n", "GhatBlackwiseRowMajor(G)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "G = Ghat2(3, 5,3,1)\n", "GhatBlackwiseRowMajor(G)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "G = Ghat2(3, 7, 5,3)\n", "GhatBlackwiseRowMajor(G)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "sc.sparse.linalg\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "for n in (8, 9):\n", " print(print_eigen_input(np.fft.fftfreq(n, 15.2/n)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "np.fft.fftfreq??" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "import scipy as sc\n", "import scipy.sparse.linalg as sp" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "sp.cg??" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "autoscroll": "json-false", "ein.tags": [ "worksheet-0" ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "name": "python3" }, "name": "material_hyper_elastic.ipynb" }, "nbformat": 4, "nbformat_minor": 2 }