#!/usr/bin/env python
# coding: utf-8

# # Strain and stress tensors in spherical coordinates
# 
# This worksheet demonstrates a few capabilities of [SageManifolds](http://sagemanifolds.obspm.fr/) (version 1.0, as included in SageMath 7.5) in computations regarding elasticity theory in Cartesian coordinates.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_elasticity_spher.ipynb) 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()


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

# In[2]:


get_ipython().run_line_magic('display', 'latex')


# ## Euclidean 3-space and spherical coordinates

# 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


# Spherical coordinates do not form a regular coordinate system of the Euclidean space. So declaring that they span $M$ means that, strictly speaking, the manifold $M$ is not the whole Euclidean space, but the Euclidean space minus some half plane (the azimuthal origin). However, in this worksheet, this difference will not matter. 

# The natural vector frame of spherical coordinates is

# In[5]:


spher.frame()


# We shall expand vector and tensor fields on the orthonormal frame $(e_1, e_2, e_3)$ associated with spherical coordinates, which is related to the natural frame $(\partial/\partial r, \partial/\partial\theta, \partial/\partial\phi)$ displayed above by means of the following field of automorphisms: 

# 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()


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

# In[7]:


to_orthonormal[:]


# We construct the orthonormal frame from the natural frame of spherical coordinates by this change of basis: 

# In[8]:


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


# In[9]:


e[1].display()


# In[10]:


e[2].display()


# In[11]:


e[3].display()


# At this stage, the default vector frame on $M$ is the first one introduced, namely the natural frame of spherical coordinates: 

# In[12]:


M.default_frame()


# 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()


# To get the same output as in `Out[10]`, one should specify the frame for display, since this is no longer the default one:

# In[15]:


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


# ## Displacement vector and strain tensor
# 
# 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()


# The following computations will involve the metric $g$ of the Euclidean space. At the current stage of SageManifolds, we need to introduce it explicitly, as a Riemannian metric on the manifold $M$ (in a future version of SageManifolds, one shall to declare $M$ as an Euclidean space, and not merely as a manifold, so that it will come equipped with $g$):

# In[17]:


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


# Since $e$ is supposed to be an orthonormal frame, we declare that the components of $g$ with respect to it are
# $\mathrm{diag}(1,1,1)$:

# In[18]:


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


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

# In[19]:


g.display(spher.frame())


# The covariant derivative operator $\nabla$ is introduced as the (Levi-Civita) connection associated with $g$: 

# In[20]:


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


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

# In[21]:


nabla.display()


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

# In[22]:


nabla.display(spher.frame())


# The covariant derivative of the displacement vector $U$ is

# In[23]:


nabU = nabla(U)
print(nabU)


# In[24]:


nabU.display()


# We convert it to a tensor field of type (0,2) (i.e. a bilinear form) by lowering the upper index with $g$:

# In[25]:


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


# In[26]:


nabU_form.display()


# 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()


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

# In[29]:


E.display_comp(only_nonredundant=True)


# ## Stress tensor and Hooke's law
# 
# To form the stress tensor according to Hooke's law, we introduce first the Lamé constants:

# In[30]:


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


# In[31]:


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


# The trace (with respect to $g$) of the bilinear form $\varepsilon$ is obtained by (i) raising the first index (`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()


# The **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()


# In[36]:


S.display_comp(only_nonredundant=True)


# Each component can be accessed individually:

# In[37]:


S[1,2]


# ## Divergence of the stress tensor
# 
# 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)


# The divergence is obtained by taking the trace on the first index (`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()


# In[41]:


divS.display_comp()


# Note that $f_1$ is quite badly displayed. We get a better view by displaying the components one by one:

# In[42]:


divS[1]


# In[43]:


divS[2]


# In[44]:


divS[3]