add generation of vel gradient scalars

meta/Velocity gradient.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n",
+      "-Axx*(Axx**2/3 + Axy*Ayx + Axz*Azx) - Axy*(Ayx*Ayy + Ayz*Azx) - Axz*(Ayx*Azy + Azx*Azz) - Ayy*(Ayy**2/3 + Ayz*Azy) - Ayz*Azy*Azz - Azz**3/3\n",
+      "Axx**2 + Axy*(Axy/2 + Ayx) + Axz*(Axz/2 + Azx) + Ayx**2/2 + Ayy**2 + Ayz*(Ayz/2 + Azy) + Azx**2/2 + Azy**2/2 + Azz**2\n"
+     ]
+    }
+   ],
+   "source": [
+    "import sympy as sp\n",
+    "\n",
+    "A = []\n",
+    "for deriv in ['x', 'y', 'z']:\n",
+    "    A.append([])\n",
+    "    for field in ['x', 'y', 'z']:\n",
+    "        A[-1].append(sp.Symbol('A' + deriv + field))\n",
+    "\n",
+    "A = sp.Matrix(A)\n",
+    "\n",
+    "A2 = A**2\n",
+    "A3 = A**3\n",
+    "Q = -sp.horner(sp.simplify(sum(A2[i, i] for i in range(3))/2))\n",
+    "R = -sp.horner(sp.simplify(sum(A3[i, i] for i in range(3))/3))\n",
+    "print(Q)\n",
+    "print(R)\n",
+    "\n",
+    "S = (A + A.T)/2\n",
+    "S2 = S**2\n",
+    "trS2 = sp.horner(sp.simplify(sum(S2[i, i] for i in range(3))))\n",
+    "print(trS2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n",
+      "3*DIV + 3*MUL + NEG + 3*POW + 5*SUB\n",
+      "-Axx**3/3 - Axx*Axy*Ayx - Axx*Axz*Azx - Axy*Ayx*Ayy - Axy*Ayz*Azx - Axz*Ayx*Azy - Axz*Azx*Azz - Ayy**3/3 - Ayy*Ayz*Azy - Ayz*Azy*Azz - Azz**3/3\n",
+      "3*DIV + 16*MUL + NEG + 3*POW + 10*SUB\n"
+     ]
+    }
+   ],
+   "source": [
+    "def alt_measure(expr):\n",
+    "    POW = sp.Symbol('POW')\n",
+    "    count = sp.count_ops(expr, visual=True).subs(POW, 10)\n",
+    "    count = count.replace(sp.Symbol, type(sp.S.One))\n",
+    "    return count\n",
+    "\n",
+    "Qalt = sp.simplify(Q, measure = alt_measure)\n",
+    "print(Qalt)\n",
+    "print(sp.count_ops(Qalt, visual=True))\n",
+    "Ralt = sp.simplify(R, measure = alt_measure)\n",
+    "print(Ralt)\n",
+    "print(sp.count_ops(Ralt, visual=True))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n",
+      "0\n"
+     ]
+    }
+   ],
+   "source": [
+    "Qalt = - sum(A[i, (i+1)%3] * A[(i+1)%3, i] for i in range(3)) - sum(A[i, i]**2 for i in range(3))/2\n",
+    "print Qalt\n",
+    "print(sp.simplify(Qalt - Q))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-Axx*(Axx**2/3 + Axy*Ayx + Axz*Azx) - Axy*Ayz*Azx - Axz*Ayx*Azy - Ayy*(Axy*Ayx + Ayy**2/3 + Ayz*Azy) - Azz*(Axz*Azx + Ayz*Azy + Azz**2/3)\n",
+      "6*ADD + 3*DIV + 13*MUL + NEG + 3*POW + 4*SUB\n",
+      "0\n"
+     ]
+    }
+   ],
+   "source": [
+    "Ralt = - (sum(A[i, i]*(A[i, i]**2/3 + sum(A[i, (i+j)%3]*A[(i+j)%3, i]\n",
+    "                                          for j in range(1, 3)))\n",
+    "              for i in range(3)) +\n",
+    "          A[0, 1]*A[1, 2]*A[2, 0] + A[0, 2]*A[1, 0]*A[2, 1])\n",
+    "print Ralt\n",
+    "print(sp.count_ops(Ralt, visual=True))\n",
+    "print(sp.simplify(Ralt - R))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0\n",
+      "0\n",
+      "0\n"
+     ]
+    }
+   ],
+   "source": [
+    "AxxAxx = A[0, 0]**2\n",
+    "AyyAyy = A[1, 1]**2\n",
+    "AzzAzz = A[2, 2]**2\n",
+    "AxyAyx = A[0, 1]*A[1, 0]\n",
+    "AyzAzy = A[1, 2]*A[2, 1]\n",
+    "AzxAxz = A[2, 0]*A[0, 2]\n",
+    "\n",
+    "Qalt = - (AxxAxx + AyyAyy + AzzAzz)/2 - AxyAyx -AyzAzy - AzxAxz\n",
+    "print(sp.simplify(Qalt - Q))\n",
+    "Ralt = - (A[0, 0]*(AxxAxx/3 + AxyAyx + AzxAxz) +\n",
+    "          A[1, 1]*(AyyAyy/3 + AxyAyx + AyzAzy) +\n",
+    "          A[2, 2]*(AzzAzz/3 + AzxAxz + AyzAzy) +\n",
+    "          A[0, 1]*A[1, 2]*A[2, 0] +\n",
+    "          A[0, 2]*A[1, 0]*A[2, 1])\n",
+    "print sp.simplify(Ralt - R)\n",
+    "trS2alt = (AxxAxx + AyyAyy + AzzAzz +\n",
+    "           ((A[0, 1] + A[1, 0])**2 +\n",
+    "            (A[1, 2] + A[2, 1])**2 +\n",
+    "            (A[2, 0] + A[0, 2])**2)/2)\n",
+    "print sp.simplify(trS2alt - trS2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "language": "python",
+   "name": "python2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}