This worksheet demonstrates a few capabilities of SageManifolds (version 1.0, as included in SageMath 7.5) in computations regarding elasticity theory in Cartesian coordinates.

Click here to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, via the command `sage -n jupyter`

*NB:* a version of SageMath at least equal to 7.5 is required to run this worksheet:

In [1]:

```
version()
```

Out[1]:

First we set up the notebook to display mathematical objects using LaTeX rendering:

In [2]:

```
%display latex
```

We introduce the Euclidean space as a 3-dimensional differentiable manifold:

In [3]:

```
M = Manifold(3, 'M', start_index=1)
print(M)
```

We shall make use of spherical coordinates $(r,\theta,\phi)$:

In [4]:

```
spher.<r,th,ph> = M.chart(r'r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
print(spher)
spher
```

Out[4]:

The natural vector frame of spherical coordinates is

In [5]:

```
spher.frame()
```

Out[5]:

In [6]:

```
to_orthonormal = M.automorphism_field()
to_orthonormal[1,1] = 1
to_orthonormal[2,2] = 1/r
to_orthonormal[3,3] = 1/(r*sin(th))
to_orthonormal.display()
```

Out[6]:

In other words, the change-of-basis matrix is

In [7]:

```
to_orthonormal[:]
```

Out[7]:

In [8]:

```
e = spher.frame().new_frame(to_orthonormal, 'e')
e
```

Out[8]:

In [9]:

```
e[1].display()
```

Out[9]:

In [10]:

```
e[2].display()
```

Out[10]:

In [11]:

```
e[3].display()
```

Out[11]:

In [12]:

```
M.default_frame()
```

Out[12]:

Since we prefer the orthonormal frame, we declare

In [13]:

```
M.set_default_frame(e)
```

Then, by default, all vector and tensor fields are displayed with respect to that frame:

In [14]:

```
e[3].display()
```

Out[14]:

`Out[10]`

, one should specify the frame for display, since this is no longer the default one:

In [15]:

```
e[3].display(spher.frame())
```

Out[15]:

Let us define the **displacement vector** $U$ in terms of its components w.r.t. the orthonormal spherical frame:

In [16]:

```
U = M.vector_field(name='U')
U[:] = [function('U_1')(r,th,ph), function('U_2')(r,th,ph),
function('U_3')(r,th,ph)]
U.display()
```

Out[16]:

In [17]:

```
g = M.riemannian_metric('g')
print(g)
```

In [18]:

```
g[1,1], g[2,2], g[3,3] = 1, 1, 1
g.display()
```

Out[18]:

The expression of $g$ with respect to the natural frame of spherical coordinates is then

In [19]:

```
g.display(spher.frame())
```

Out[19]:

In [20]:

```
nabla = g.connection()
print(nabla)
nabla
```

Out[20]:

The connection coefficients with respect to the spherical orthonormal frame $e$ are

In [21]:

```
nabla.display()
```

Out[21]:

while those with respect to the natural frame of spherical coordinates (Christoffel symbols) are:

In [22]:

```
nabla.display(spher.frame())
```

Out[22]:

The covariant derivative of the displacement vector $U$ is

In [23]:

```
nabU = nabla(U)
print(nabU)
```

In [24]:

```
nabU.display()
```

Out[24]:

In [25]:

```
nabU_form = nabU.down(g)
print(nabU_form)
```

In [26]:

```
nabU_form.display()
```

Out[26]:

The **strain tensor** $\varepsilon$ is defined as the symmetrized part of this tensor:

In [27]:

```
E = nabU_form.symmetrize()
print(E)
```

In [28]:

```
E.set_name('E', latex_name=r'\varepsilon')
E.display()
```

Out[28]:

Let us display the components of $\varepsilon$, skipping those that can be deduced by symmetry:

In [29]:

```
E.display_comp(only_nonredundant=True)
```

Out[29]:

To form the stress tensor according to Hooke's law, we introduce first the LamÃ© constants:

In [30]:

```
var('ll', latex_name=r'\lambda')
```

Out[30]:

In [31]:

```
var('mu', latex_name=r'\mu')
```

Out[31]:

`pos=0`

) by means of $g$ and (ii) by taking the trace of the resulting endomorphism:

In [32]:

```
trE = E.up(g, pos=0).trace()
print(trE)
```

In [33]:

```
trE.display()
```

Out[33]:

**stress tensor** $S$ is obtained via Hooke's law for isotropic material:
$$ S = \lambda \, \mathrm{tr}\varepsilon \; g + 2\mu \, \varepsilon$$

In [34]:

```
S = ll*trE*g + 2*mu*E
print(S)
```

In [35]:

```
S.set_name('S')
S.display()
```

Out[35]:

In [36]:

```
S.display_comp(only_nonredundant=True)
```

Out[36]:

Each component can be accessed individually:

In [37]:

```
S[1,2]
```

Out[37]:

The divergence of the stress tensor (with respect to $g$) is the 1-form:
$$ f_i = \nabla_j S^j_{\ \, i} $$
In a next version of SageManifolds, there will be a function `divergence()`

. For the moment, to evaluate $f$,
we first form the tensor $S^j_{\ \, i}$ by raising the first index (`pos=0`

) of $S$ with $g$:

In [38]:

```
SU = S.up(g, pos=0)
print(SU)
```

`0`

) and the third one (`2`

) of the tensor
$(\nabla S)^j_{\ \, ik} = \nabla_k S^j_{\ \, i}$:

In [39]:

```
divS = nabla(SU).trace(0,2)
print(divS)
```

In [40]:

```
divS.set_name('f')
divS.display()
```

Out[40]:

In [41]:

```
divS.display_comp()
```

Out[41]:

In [42]:

```
divS[1]
```

Out[42]:

In [43]:

```
divS[2]
```

Out[43]:

In [44]:

```
divS[3]
```

Out[44]: