Demonstration of covariant derivatives in sympy branch

Covariant_derivative_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc3ddb64-0824-4100-8f61-3c1f2afe1299",
   "metadata": {},
   "source": [
    "This is a motivating example for https://github.com/Jason-S-Ross/sympy/tree/add-covariant-derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de4dca3-27c1-431e-ae42-9f210da13f8a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Problem Statement\n",
    "\n",
    "In the theory of elasticity, when considering infinitesimal deformations of isotropic bodies, the stress tensor is given in terms of the strain tensor by the following equation:\n",
    "\n",
    "$$\\frac{\\tau^{ij}}{2 \\mu} = \\gamma^{ij} + \\frac{\\eta}{1 - 2 \\eta} g^{ij}\\gamma^r_r$$\n",
    "\n",
    "where $\\mu$ is the shear modulus and $\\eta$ is Poisson's ratio.\n",
    "\n",
    "The strain tensor can be expressed in terms of the displacements with the following equation:\n",
    "\n",
    "$$\\gamma_{ij} = \\frac{1}{2}\\left(\\left.v_i\\right|_j + \\left.v_j\\right|_i\\right)$$\n",
    "\n",
    "The compatibility of displacements requires that\n",
    "\n",
    "$$\\left(1 - 2 \\eta \\right) \\left.v^i\\right|^j_j + \\left. v^j\\right|^i_j = 0$$\n",
    "\n",
    "* Show that the equation $2 \\mu v_i = \\left. F \\right|_i$ satisfies the compatibility equation when $F$ is a harmonic function.\n",
    "* Find the components of the stress tensor in cartesian, cylindrical, and spherical coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "694e8d6d-2fed-44ba-8c34-6aec8cf53868",
   "metadata": {},
   "source": [
    "# Procedure\n",
    "\n",
    "First, create Sympy expressions for the above equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "02e2d3dc-0514-4f8f-9651-2dfc42f7847a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(1 - 2 \\eta\\right)\\left.v{}^{i}\\right|{{}^{E_{0}}{}_{E_{0}}} + \\left.v{}^{E_{0}}\\right|{{}^{i}{}_{E_{0}}} = 0$"
      ],
      "text/plain": [
       "Eq(TensDiff(v(E_0), [i, -E_0]) + (1 - 2*eta)*TensDiff(v(i), [E_0, -E_0]), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\gamma{}_{ij} = \\left(\\frac{1}{2}\\right)\\left(\\left.v{}_{i}\\right|{{}_{j}} + \\left.v{}_{j}\\right|{{}_{i}}\\right)$"
      ],
      "text/plain": [
       "Eq(gamma(-i, -j), (1/2)*(TensDiff(v(-i), [-j]) + TensDiff(v(-j), [-i])))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(\\frac{1}{2}\\right)\\frac{1}{\\mu}\\tau{}^{ij} = \\eta\\frac{1}{1 - 2 \\eta}g{}^{ij}\\gamma{}^{E_{0}}{}_{E_{0}} + \\gamma{}^{ij}$"
      ],
      "text/plain": [
       "Eq((1/2)*1/mu*tau(i, j), eta*1/(1 - 2*eta)*g(i, j)*gamma(E_0, -E_0) + gamma(i, j))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle v{}_{i} = \\left(\\frac{1}{2}\\right)\\frac{1}{\\mu}\\left.F\\right|{{}_{i}}$"
      ],
      "text/plain": [
       "Eq(v(-i), (1/2)*1/mu*TensDiff(F, [-i]))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import eta, mu, phi, rho, theta, x, y, z\n",
    "from sympy.tensor.tensor import (\n",
    "    TensorCoordinateSystem,\n",
    "    TensorIndexType,\n",
    "    tensor_heads,\n",
    "    tensor_indices,\n",
    ")\n",
    "\n",
    "F_func = symbols(\"F\", cls=Function)\n",
    "Euclid = TensorIndexType(\"Euclid\", dummy_name=\"E\")\n",
    "i, j, r = tensor_indices(\"i j r\", Euclid)\n",
    "\n",
    "g = tensor_heads(\"g\", [Euclid, Euclid])\n",
    "F_tens = tensor_heads(\"F\", [])\n",
    "v = tensor_heads(\"v\", [Euclid])\n",
    "g, gamma, tau = tensor_heads(\"g gamma tau\", [Euclid, Euclid])\n",
    "\n",
    "compatibility = Eq((1 - 2 * eta) * v(i).covar_diff(j, -j) + v(j).covar_diff(i, -j), 0)\n",
    "display(compatibility)\n",
    "strain = Eq(gamma(-i, -j), (v(-i).covar_diff(-j) + v(-j).covar_diff(-i)) / 2)\n",
    "display(strain)\n",
    "stress = Eq(\n",
    "    tau(i, j) / (2 * mu), gamma(i, j) + eta / (1 - 2 * eta) * g(i, j) * gamma(r, -r)\n",
    ")\n",
    "display(stress)\n",
    "displacement = Eq(v(-i), F_tens().covar_diff(-i) / (2 * mu))\n",
    "display(displacement)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf815b8-7b9c-46c6-85ea-578ba0470339",
   "metadata": {},
   "source": [
    "Next, show that the conditions of compatibility are satisfied by a harmonic function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a55eb7d5-0bd4-4a33-ac21-cc3316d2b0f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(1 - 2 \\eta\\right)\\left.F\\right|{{}^{iE_{0}}{}_{E_{0}}} + \\left(\\frac{1}{2}\\right)\\frac{1}{\\mu}\\left.F\\right|{{}^{E_{0}i}{}_{E_{0}}} = 0$"
      ],
      "text/plain": [
       "Eq((1 - 2*eta)*TensDiff(F, [i, E_0, -E_0]) + (1/2)*1/mu*TensDiff(F, [E_0, i, -E_0]), 0)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compatibility.subs(displacement.lhs, displacement.rhs).doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36e84240-1c47-47aa-8564-eddffc94a052",
   "metadata": {},
   "source": [
    "Currently, I don't have a way to express whether or not the order of covariant derivatives matters, so it's a manual process to verify that the above equations are true when $F$ is harmonic.  This proceeds as follows:\n",
    "\n",
    "The Laplacian operation $\\nabla^2 F$ can be written in tensor notation as $\\left.F\\right|^i_i$. \n",
    "Since we are in a euclidean space, the order of partial differentiation does not matter, so we can see that we have a Laplacian operation inside of each term. Therefore, the condition is satisfied if $F$ is harmonic.\n",
    "\n",
    "Next, substitute our stress function into the stress equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32fcf611-734e-4c27-993d-b06fbd9fe13f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(\\frac{1}{2}\\right)\\frac{1}{\\mu}\\tau{}^{ij} = \\left(\\frac{1}{4}\\right)\\eta\\frac{1}{\\mu}\\frac{1}{1 - 2 \\eta}g{}^{ij}\\left.F\\right|{{}^{E_{0}}{}_{E_{0}}} + \\left(\\frac{1}{4}\\right)\\eta\\frac{1}{\\mu}\\frac{1}{1 - 2 \\eta}g{}^{ij}\\left.F\\right|{{}_{E_{0}}{}^{E_{0}}} + \\left(\\frac{1}{4}\\right)\\frac{1}{\\mu}\\left.F\\right|{{}^{ij}} + \\left(\\frac{1}{4}\\right)\\frac{1}{\\mu}\\left.F\\right|{{}^{ji}}$"
      ],
      "text/plain": [
       "Eq((1/2)*1/mu*tau(i, j), (1/4)*1/mu*TensDiff(F, [i, j]) + (1/4)*1/mu*TensDiff(F, [j, i]) + (1/4)*eta*1/mu*1/(1 - 2*eta)*g(i, j)*TensDiff(F, [-E_0, E_0]) + (1/4)*eta*1/mu*1/(1 - 2*eta)*g(i, j)*TensDiff(F, [E_0, -E_0]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stress_new = Eq(\n",
    "    stress.lhs,\n",
    "    stress.rhs.subs(strain.lhs, strain.rhs)\n",
    "    .subs(displacement.lhs, displacement.rhs)\n",
    "    .doit()\n",
    "    .expand(),\n",
    ")\n",
    "stress_new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f97d1b08-a9d0-42f3-9c69-ca450583f000",
   "metadata": {},
   "source": [
    "We can simplify further by noting that $F$ is harmonic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5d191407-d08e-438f-a592-dad2d9706d38",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\tau{}^{ij} = \\left(\\frac{1}{2}\\right)\\left.F\\right|{{}^{ij}} + \\left(\\frac{1}{2}\\right)\\left.F\\right|{{}^{ji}}$"
      ],
      "text/plain": [
       "Eq(tau(i, j), (1/2)*TensDiff(F, [i, j]) + (1/2)*TensDiff(F, [j, i]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stress_new = Eq(stress_new.lhs, stress_new.rhs.subs(F_tens().covar_diff(i, -i), S.Zero)).doit()\n",
    "stress_new = Eq(stress_new.lhs * 2 * mu, stress_new.rhs * 2 * mu).expand()\n",
    "stress_new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58203100-2001-4d25-a4ec-6d32ea7ad2f0",
   "metadata": {},
   "source": [
    "Right now, I can't simplify the above equation any more because I don't have a way to express that the order of covariant derivatives doesn't matter in this space. This could probably be a new property of `TensorIndexType`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdc32e96-f553-4ce4-8b7d-15b0c7dd19f2",
   "metadata": {},
   "source": [
    "# Stress Tensor Components"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "238e4594-1246-4fdf-8e52-b08fd13f759a",
   "metadata": {},
   "source": [
    "In order to evaluate covariant derivatives, a new class is used called `TensorCoordinateSystem`. This contains the metric, the coordinate variables, and the Christoffel symbols in a single object. To use it, substitute for the assocated `TensorIndexType` in the replacement dict, instead of the metric array as before (using an array still works, but won't let you take covariant derivatives)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7570b707-2762-4b4e-9274-c399f5b6ccc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\partial^{2}}{\\partial x^{2}} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial y\\partial x} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial z\\partial x} F{\\left(x,y,z \\right)}\\\\\\frac{\\partial^{2}}{\\partial y\\partial x} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial y^{2}} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial z\\partial y} F{\\left(x,y,z \\right)}\\\\\\frac{\\partial^{2}}{\\partial z\\partial x} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial z\\partial y} F{\\left(x,y,z \\right)} & \\frac{\\partial^{2}}{\\partial z^{2}} F{\\left(x,y,z \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[Derivative(F(x, y, z), (x, 2)), Derivative(F(x, y, z), x, y), Derivative(F(x, y, z), x, z)], [Derivative(F(x, y, z), x, y), Derivative(F(x, y, z), (y, 2)), Derivative(F(x, y, z), y, z)], [Derivative(F(x, y, z), x, z), Derivative(F(x, y, z), y, z), Derivative(F(x, y, z), (z, 2))]]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cartesian = TensorCoordinateSystem(\"R\", symbols=[x, y, z])  # Metric is identity by default\n",
    "cartesian_replacement = {\n",
    "    Euclid: cartesian,\n",
    "    F_tens(): F_func(x, y, z),\n",
    "    g(i, -j): eye(3),  # Hack to get the metric tensor to appear in the expression\n",
    "}\n",
    "simplify(stress_new.rhs.replace_with_arrays(cartesian_replacement))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cce6561d-9f55-4da6-801c-4a52c3bea2c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\partial^{2}}{\\partial \\rho^{2}} F{\\left(\\rho,\\theta,z \\right)} & \\frac{\\rho \\frac{\\partial^{2}}{\\partial \\theta\\partial \\rho} F{\\left(\\rho,\\theta,z \\right)} - \\frac{\\partial}{\\partial \\theta} F{\\left(\\rho,\\theta,z \\right)}}{\\rho^{3}} & \\frac{\\partial^{2}}{\\partial z\\partial \\rho} F{\\left(\\rho,\\theta,z \\right)}\\\\\\frac{\\rho \\frac{\\partial^{2}}{\\partial \\theta\\partial \\rho} F{\\left(\\rho,\\theta,z \\right)} - \\frac{\\partial}{\\partial \\theta} F{\\left(\\rho,\\theta,z \\right)}}{\\rho^{3}} & \\frac{\\rho \\frac{\\partial}{\\partial \\rho} F{\\left(\\rho,\\theta,z \\right)} + \\frac{\\partial^{2}}{\\partial \\theta^{2}} F{\\left(\\rho,\\theta,z \\right)}}{\\rho^{4}} & \\frac{\\frac{\\partial^{2}}{\\partial z\\partial \\theta} F{\\left(\\rho,\\theta,z \\right)}}{\\rho^{2}}\\\\\\frac{\\partial^{2}}{\\partial z\\partial \\rho} F{\\left(\\rho,\\theta,z \\right)} & \\frac{\\frac{\\partial^{2}}{\\partial z\\partial \\theta} F{\\left(\\rho,\\theta,z \\right)}}{\\rho^{2}} & \\frac{\\partial^{2}}{\\partial z^{2}} F{\\left(\\rho,\\theta,z \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[Derivative(F(rho, theta, z), (rho, 2)), (rho*Derivative(F(rho, theta, z), rho, theta) - Derivative(F(rho, theta, z), theta))/rho**3, Derivative(F(rho, theta, z), rho, z)], [(rho*Derivative(F(rho, theta, z), rho, theta) - Derivative(F(rho, theta, z), theta))/rho**3, (rho*Derivative(F(rho, theta, z), rho) + Derivative(F(rho, theta, z), (theta, 2)))/rho**4, Derivative(F(rho, theta, z), theta, z)/rho**2], [Derivative(F(rho, theta, z), rho, z), Derivative(F(rho, theta, z), theta, z)/rho**2, Derivative(F(rho, theta, z), (z, 2))]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cylindrical = TensorCoordinateSystem(\n",
    "    \"C\",\n",
    "    symbols=[rho, theta, z],\n",
    "    metric=Array([[1, 0, 0], [0, rho**2, 0], [0, 0, 1]])\n",
    ")\n",
    "cylindrical_replacement = {\n",
    "    Euclid: cylindrical,\n",
    "    F_tens(): F_func(rho, theta, z),\n",
    "    g(i, -j): eye(3),  # Hack to get the metric tensor to appear in the expression\n",
    "}\n",
    "simplify(stress_new.rhs.replace_with_arrays(cylindrical_replacement))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1f218207-dd0c-472c-be93-138e0e167b9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\partial^{2}}{\\partial \\rho^{2}} F{\\left(\\rho,\\theta,\\phi \\right)} & \\frac{\\rho \\frac{\\partial^{2}}{\\partial \\theta\\partial \\rho} F{\\left(\\rho,\\theta,\\phi \\right)} - \\frac{\\partial}{\\partial \\theta} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{3}} & \\frac{\\rho \\frac{\\partial^{2}}{\\partial \\rho\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)} - \\frac{\\partial}{\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{3} \\sin^{2}{\\left(\\theta \\right)}}\\\\\\frac{\\rho \\frac{\\partial^{2}}{\\partial \\theta\\partial \\rho} F{\\left(\\rho,\\theta,\\phi \\right)} - \\frac{\\partial}{\\partial \\theta} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{3}} & \\frac{\\rho \\frac{\\partial}{\\partial \\rho} F{\\left(\\rho,\\theta,\\phi \\right)} + \\frac{\\partial^{2}}{\\partial \\theta^{2}} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{4}} & \\frac{\\sin{\\left(\\theta \\right)} \\frac{\\partial^{2}}{\\partial \\theta\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)} - \\cos{\\left(\\theta \\right)} \\frac{\\partial}{\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{4} \\sin^{3}{\\left(\\theta \\right)}}\\\\\\frac{\\rho \\frac{\\partial^{2}}{\\partial \\rho\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)} - \\frac{\\partial}{\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{3} \\sin^{2}{\\left(\\theta \\right)}} & \\frac{\\sin{\\left(\\theta \\right)} \\frac{\\partial^{2}}{\\partial \\theta\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)} - \\cos{\\left(\\theta \\right)} \\frac{\\partial}{\\partial \\phi} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{4} \\sin^{3}{\\left(\\theta \\right)}} & \\frac{\\rho \\sin^{2}{\\left(\\theta \\right)} \\frac{\\partial}{\\partial \\rho} F{\\left(\\rho,\\theta,\\phi \\right)} + \\frac{\\sin{\\left(2 \\theta \\right)} \\frac{\\partial}{\\partial \\theta} F{\\left(\\rho,\\theta,\\phi \\right)}}{2} + \\frac{\\partial^{2}}{\\partial \\phi^{2}} F{\\left(\\rho,\\theta,\\phi \\right)}}{\\rho^{4} \\sin^{4}{\\left(\\theta \\right)}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[[Derivative(F(rho, theta, phi), (rho, 2)), (rho*Derivative(F(rho, theta, phi), rho, theta) - Derivative(F(rho, theta, phi), theta))/rho**3, (rho*Derivative(F(rho, theta, phi), phi, rho) - Derivative(F(rho, theta, phi), phi))/(rho**3*sin(theta)**2)], [(rho*Derivative(F(rho, theta, phi), rho, theta) - Derivative(F(rho, theta, phi), theta))/rho**3, (rho*Derivative(F(rho, theta, phi), rho) + Derivative(F(rho, theta, phi), (theta, 2)))/rho**4, (sin(theta)*Derivative(F(rho, theta, phi), phi, theta) - cos(theta)*Derivative(F(rho, theta, phi), phi))/(rho**4*sin(theta)**3)], [(rho*Derivative(F(rho, theta, phi), phi, rho) - Derivative(F(rho, theta, phi), phi))/(rho**3*sin(theta)**2), (sin(theta)*Derivative(F(rho, theta, phi), phi, theta) - cos(theta)*Derivative(F(rho, theta, phi), phi))/(rho**4*sin(theta)**3), (rho*sin(theta)**2*Derivative(F(rho, theta, phi), rho) + sin(2*theta)*Derivative(F(rho, theta, phi), theta)/2 + Derivative(F(rho, theta, phi), (phi, 2)))/(rho**4*sin(theta)**4)]]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spherical = TensorCoordinateSystem(\n",
    "    \"S\",\n",
    "    [rho, theta, phi],\n",
    "    Array([[1, 0, 0], [0, rho ** 2, 0], [0, 0, rho ** 2 * sin(theta) ** 2]]),\n",
    ")\n",
    "spherical_replacement = {\n",
    "    Euclid: spherical,\n",
    "    F_tens(): F_func(rho, theta, phi),\n",
    "    g(i, -j): eye(3),  # Hack to get the metric tensor to appear in the expression\n",
    "}\n",
    "\n",
    "simplify(stress_new.rhs.replace_with_arrays(spherical_replacement))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Interesting. Do you plan on merging your changes back to SymPy?

I would need some help to get this merged. I don't have a strong enough grasp of tensor calculus to know if all the assumptions I'm making are correct.

sympy.tensor.tensor was originally meant to only support the Penrose abstract index notation (which is independent of coordinate systems). Unfortunately that was too restrictive as partial derivatives cannot be represented. Adding support for indices with coordinate systems is quite complicated OTOH.

Another great shortcoming of sympy.tensor.tensor is that it is complicated to use, and its applications are quite limited (physics and Riemannian geometry only).

I have since switched to the development of sympy.tensor.array.expressions, which is meant to be more general as possible (it's just N-dimensional arrays, not even tensors, but it's probably useful to many more users of SymPy).

Also have a look at sympy.diffgeom, that's a quite robust implementation of Riemannian geometry.

sympy.tensor.tensor was originally meant to only support the Penrose abstract index notation (which is independent of coordinate systems). Unfortunately that was too restrictive as partial derivatives cannot be represented. Adding support for indices with coordinate systems is quite complicated OTOH.

To my understanding this does add support for indices with coordinate systems, but probably not in a rigorous way.

Another great shortcoming of sympy.tensor.tensor is that it is complicated to use, and its applications are quite limited (physics and Riemannian geometry only).

I see. To me as an engineering student, physics and Riemannian geometry are basically the entire universe. I didn't find the tensor module to be overly complicated to use, for what it's worth.

Also have a look at sympy.diffgeom, that's a quite robust implementation of Riemannian geometry.

I tried that module first, but I don't have any background in differential geometry and I didn't understand how I could use it to solve problems like this. I asked the mailing list and didn't get much help, so I wrote this.

---

Unfortunately I am ignorant of the extent of my own ignorance of this subject, so I'm not qualified to speculate on whether or not this would be good for SymPy. I assume I have made some crucial mistakes in my assumptions that would not be appropriate for the level of generality desired for the tensor module. If not, I would be happy to submit a pull request.