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

# # Strain and stress tensors in Cartesian 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_Cartesian.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 Cartesian coordinates

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

# In[3]:


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


# We then introduce the Cartesian coordinates $(x,y,z)$ as a chart $X$ on $M$:

# In[4]:


X.<x,y,z> = M.chart()
print(X)
X


# The associated vector frame is

# In[5]:


X.frame()


# We shall expand vector and tensor fields not on this frame, which is the default one on $M$:

# In[6]:


M.default_frame()


# ## Displacement vector and strain tensor
# 
# Let us define the **displacement vector** $U$ in terms of its components w.r.t. the orthonormal Cartesian frame:

# In[7]:


U = M.vector_field(name='U')
U[:] = [function('U_x')(x,y,z), function('U_y')(x,y,z), 
        function('U_z')(x,y,z)]
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[8]:


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


# We initialize $g$ by declaring that its components with respect to the frame of Cartesian coordinates are
# $\mathrm{diag}(1,1,1)$:

# In[9]:


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


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

# In[10]:


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


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

# In[11]:


nabU = nabla(U)
print(nabU)


# In[12]:


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[13]:


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


# In[14]:


nabU_form.display()


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

# In[15]:


E = nabU_form.symmetrize()
print(E)


# In[16]:


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[17]:


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[18]:


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


# In[19]:


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[20]:


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


# In[21]:


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[22]:


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


# In[23]:


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


# In[24]:


S.display_comp(only_nonredundant=True)


# Each component can be accessed individually:

# In[25]:


S[1,2]


# ## Divergence of the stress tensor
# 
# The divergence of the stress tensor 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[26]:


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[27]:


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


# In[28]:


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


# In[29]:


divS.display_comp()


# Displaying the components one by one:

# In[30]:


divS[1]


# In[31]:


divS[2]


# In[32]:


divS[3]


# In[ ]: