Tensor Calculus

[Jason Ross]

Hello,

(I'm having problems getting this to format in google groups so here's a link to a gist: https://gist.github.com/Jason-S-Ross/11336ddb94f59cf8377e1ffe09df1db5)

I'm trying to use Sympy to help solve problems in solid mechanics, where I have
a lot of tensor differential equations, expressed in arbitrary curvilinear
coordinates in euclidean space.

Here's an example problem:

> The compatibility of displacements requires that
> $$\left(1 - 2\eta\right) g^{rj}v^i|_{jr} + g^{ri} v^j|_{jr} = 0$$
> Show that the following is displacement field satisfies the
> compatibility requirement and compute the components of stress
> in spherical coordinates:
> $$2 \mu v_i = \phi|_i$$
> where $\phi$ is a scalar field that is harmonic, i.e. $(g^{ir}\phi|_r)|_i=0$
> and the stress $\tau^{ij}$ is given by the following equation:
> $$\frac{\tau^{ij}}{\mu} = g^{js} v^i|_s + g^{ir} v^j|_r + \frac{2\eta}{1 - 2\eta}g^{ij} v^r|_r$$

This creates a few challenges:

* How do we manipulate differential equations involving covariant derivatives
 like $\left(1 - 2\eta\right) g^{rj}v^i|_{jr} + g^{ri} v^j|_{jr} = 0$
 symbolically, to show that a given equation like $v_i = \frac{1}{2\mu}\phi|_i$
 satisfies the differential equation without abandoning the index notation?
* How do we take a known solution like $\phi = \frac{1}{\sqrt{x^2 + y^2 + z^2}}$
and compute the corresponding stress in an arbitrary coordinate system like
spherical coordinates? This involves computing the metric tensor for the
coordinate system, computing Christoffel symbols, and evaluating covariant
derivatives, all of which are possible, but the Tensor api doesn't seem to
support.
* How do we operate on tensor expressions involving the metric tensor itself,
when `replace_with_arrays` relies on the metric being a property of the
`TensorIndexType`? For instance, consider the following:

```
    from sympy import *
    from sympy.tensor.tensor import TensorHead, TensorSymmetry, TensorIndexType, tensor_indices
    from sympy.abc import eta
    euclid = TensorIndexType('Euclid', dummy_name='R')
    i, j, r, s = tensor_indices('i j r s', euclid)
    v = TensorHead('v', [euclid])
    # Hack to handle the first covariant derivative of v; populate with values "by hand"
    v1 = TensorHead('v', [euclid, euclid])
    # Hack to handle the second covariant derivative of v; populate with values "by hand"
    v2 = TensorHead('v', [euclid, euclid, euclid])
    g = TensorHead('g', [euclid, euclid], TensorSymmetry.fully_symmetric(2))
    expr = g(j, s) * v1(i, -s) + g(i, r) * v1(j, -r) + 2 * eta / (1 - 2 * eta) * g(i, j) * v1(r, -r)
```
How do we evaluate this expression in spherical coordinates? `g` is the metric
tensor, but the substitution we give in the replacement dictionary for `euclid`
is what's actually used to raise and lower indices.
* How do we compute Christoffel symbols for tensor expressions? They're part of
the `diffgeom` module, but they're not represented there in a way that's
compatible with the `tensor` module.

In light of these I'm forced to work with piecemeal solutions that require
finicky expression manipulation using combinations of `tensorproduct` and
`tensorcontraction` and lose a lot of the expressiveness of tensor notation in
the process. Can anyone provide some pointers on working with tensor
differential equations in sympy?

[Nicolas Guarin]

I developed a package for continuum mechanics (https://continuum-mechanics.readthedocs.io/en/latest/) on top of SymPy. It uses a different approach, though, since I write the expressions in coordinates explicitly. I tried to write on top of the Mechanics module, but didn't understand the use of coordinate systems and dyads at the moment (I haven't look in a while). The same goes for the tensor module, I find it cumbersome for most of the calculations.

[Jason Ross]

This looks like an interesting application. Thanks for letting me know!