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

# # Manifold tutorial
# 
# This notebook provides a short introduction to differentiable manifolds in SageMath. The tools described below have been implemented through the [SageManifolds](https://sagemanifolds.obspm.fr) project.

# This notebook is valid for version 9.2 or higher of SageMath:

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


# ## Defining a manifold
# 
# As an example let us define a differentiable manifold of dimension 3 over $\mathbb{R}$:

# In[3]:


M = Manifold(3, 'M', latex_name=r'\mathcal{M}', start_index=1)


# - The first argument, `3`, is the manifold dimension; it can be any
#   positive integer.
# - The second argument, `'M'`, is a string defining the manifold's name; it may be 
#   different from the symbol set on the left-hand side of the = sign (here `M`): the latter
#   stands for the Python variable that refers to the manifold object in the computer 
#   memory, while the string `'M'` is the mathematical symbol chosen for the manifold.
# - The optional argument `latex_name=r'\mathcal{M}'` sets the LaTeX
#   symbol to display the manifold. Note the letter 'r' in front on the first quote: 
#   it indicates that the  string is a *raw* one, so that the backslash character 
#   in `\mathcal` is considered as an ordinary character (otherwise, the backslash is 
#   used to escape some special characters). If the argument `latex_name` is not 
#   provided by the user, it is set to the string used as the second argument (here `'M'`)
# - The optional argument `start_index=1` defines the range of indices to be used for 
#   tensor components on the manifold: setting it to 1 means that indices will range 
#   in $\{1,2,3\}$. The default value is `start_index=0`.
# 
# Note that the default base field is $\mathbb{R}$. If we would have used the optional
# argument `field='complex'`, we would have defined a manifold over $\mathbb{C}$. See the
# [list of all options](http://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/manifold.html#sage.manifolds.manifold.Manifold) for more details. 
# 
# If we ask for M, it is displayed via its LaTeX symbol:

# In[4]:


M


# If we use the function `print()` instead, we get a short description of the object:

# In[5]:


print(M)


# Via the function `type()`, we get the type of the Python object corresponding to M (here the Python class `DifferentiableManifold_with_category`):

# In[6]:


print(type(M))


# We can also ask for the category of M and see that it is the category of smooth manifolds over $\mathbb{R}$:

# In[7]:


category(M)


# The indices on the manifold are generated by the method `irange()`, to be used in loops:

# In[8]:


[i for i in M.irange()]


# If the parameter `start_index` had not been specified, the default range of the indices would have been $\{0,1,2\}$ instead:

# In[9]:


M0 = Manifold(3, 'M', latex_name=r'\mathcal{M}')
[i for i in M0.irange()]


# ## Defining a chart on the manifold
# 
# Let us assume that the manifold $\mathcal{M}$ can be covered by a single chart (other cases are discussed below); the chart is declared as follows:

# In[10]:


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


# The writing `.<x,y,z>` in the left-hand side means that the Python variables `x`, `y` and `z` are set to the three coordinates of the chart. This allows one to refer subsequently to the coordinates by their names.
# 
# In this example, the function `chart()` has no arguments, which implies that the coordinate symbols will be `x`, `y` and `z` (i.e. exactly the characters set in the `<...>` operator) and that each coordinate range is $(-\infty,+\infty)$. For other cases, an argument must be passed to `chart()`  to specify the coordinate symbols and range, as well as the LaTeX symbol of a coordinate if the latter is different from the coordinate name (an example will be provided below).

# <p>The chart is displayed as a pair formed by the open set covered by it (here the whole manifold) and the coordinates:</p>

# In[11]:


print(X)


# In[12]:


X


# The coordinates can be accessed individually, by means of their indices, following the convention defined by `start_index=1` in the manifold's definition:

# In[13]:


X[1]


# In[14]:


X[2]


# In[15]:


X[3]


# The full set of coordinates is obtained by means of the operator `[:]`:

# In[16]:


X[:]


# Thanks to the operator `<x,y,z>` used in the chart declaration, each coordinate can be accessed directly via its name:

# In[17]:


z is X[3]


# Coordinates are SageMath symbolic expressions:

# In[18]:


print(type(z))


# ### Functions of the chart coordinates
# 
# Real-valued functions of the chart coordinates (mathematically speaking, *functions defined on the chart codomain*) are generated via the method `function()` acting on the chart:

# In[19]:


f = X.function(x+y^2+z^3)
f


# In[20]:


f.display()


# In[21]:


f(1,2,3)


# They belong to the class `ChartFunction` (actually the subclass `𝙲𝚑𝚊𝚛𝚝𝙵𝚞𝚗𝚌𝚝𝚒𝚘𝚗𝚁𝚒𝚗𝚐_𝚠𝚒𝚝h_𝚌𝚊𝚝𝚎𝚐𝚘𝚛𝚢.𝚎𝚕𝚎𝚖𝚎𝚗t_𝚌𝚕𝚊𝚜𝚜`):

# In[22]:


print(type(f))


# and differ from SageMath standard symbolic functions by automatic simplifications in all operations. For instance, adding the two symbolic functions

# In[23]:


f0(x,y,z) = cos(x)^2; g0(x,y,z) = sin(x)^2


# results in

# In[24]:


f0 + g0


# while the sum of the corresponding functions in the class `ChartFunction` is automatically simplified:

# In[25]:


f1 = X.function(cos(x)^2); g1 = X.function(sin(x)^2)
f1 + g1


# To get the same output with symbolic functions, one has to invoke the method `simplify_trig()`:

# In[26]:


(f0 + g0).simplify_trig()


# Another difference regards the display; if we ask for the symbolic function `f0`, we get

# In[27]:


f0


# while if we ask for the chart function `f1`, we get only the coordinate expression:

# In[28]:


f1


# To get an output similar to that of `f0`, one should call the method `display()`:

# In[29]:


f1.display()


# Note that the method `expr()` returns the underlying symbolic expression:

# In[30]:


f1.expr()


# In[31]:


print(type(f1.expr()))


# ### Introducing a second chart on the manifold
# 
# Let us first consider an open subset of $\mathcal{M}$, for instance the complement $U$ of the region defined by $\{y=0, x\geq 0\}$ (note that `(y!=0, x<0)` stands for $y\not=0$ OR $x<0$; the condition $y\not=0$ AND $x<0$ would have been written `[y!=0, x<0]` instead):

# In[32]:


U = M.open_subset('U', coord_def={X: (y!=0, x<0)})


# Let us call `X_U` the restriction of the chart `X` to the open subset $U$:

# In[33]:


X_U = X.restrict(U)
X_U


# <p>We introduce another chart on $U$, with spherical-type coordinates $(r,\theta,\phi)$:</p>

# In[34]:


Y.<r,th,ph> = U.chart(r'r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
Y


# The method `chart()` is now used with an argument; it is a string, which contains specific LaTeX symbols, hence the prefix 'r' to it (for *raw* string). It also contains the coordinate ranges, since they are different from the default value, which is $(-\infty, +\infty)$. For a given coordinate, the various fields are separated by the character ':' and a space character separates the coordinates. Note that for the coordinate $r$, there are only two fields, since the LaTeX symbol has not to be specified. The LaTeX symbols are used for the outputs:

# In[35]:


th, ph


# In[36]:


Y[2], Y[3]


# The declared coordinate ranges are now known to Sage, as we may check by means of the command `assumptions()`:

# In[37]:


assumptions()


# <p>They are used in simplifications:</p>

# In[38]:


simplify(abs(r))


# In[39]:


simplify(abs(x)) # no simplification occurs since x can take any value in R


# <p>After having been declared, the chart Y can be fully specified by its relation to the chart X_U, via a transition map:</p>

# In[40]:


transit_Y_to_X = Y.transition_map(X_U, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
transit_Y_to_X


# In[41]:


transit_Y_to_X.display()


# The inverse of the transition map can be specified by means of the method `set_inverse()`:

# In[42]:


transit_Y_to_X.set_inverse(sqrt(x^2+y^2+z^2), atan2(sqrt(x^2+y^2),z), atan2(y, x))


# A check of the provided inverse is performed by composing it with the original transition map, on the left and on the right respectively. As indicated, the reported failure for `th` and `ph` is actually due to a lack of simplification of expressions involving `arctan2`.  

# We have then

# In[43]:


transit_Y_to_X.inverse().display()


# <p>At this stage, the manifold's <strong>atlas</strong> (the "user atlas", not the maximal atlas!) contains three charts:</p>

# In[44]:


M.atlas()


# The first chart defined on the manifold is considered as the manifold's default chart (this can be changed by the method `set_default_chart()`):

# In[45]:


M.default_chart()


# <p>Each open subset has its own atlas (since an open subset of a manifold is a manifold by itself):</p>

# In[46]:


U.atlas()


# In[47]:


U.default_chart()


# We can draw the chart $Y$ in terms of the chart $X$ via the command `Y.plot(X)`, which shows the lines of constant coordinates from the $Y$ chart in a "Cartesian frame" based on the $X$ coordinates:

# In[48]:


Y.plot(X)


# The method `plot()` allows for many options, to control the number of coordinate lines to be drawn, their style and color, as well as the coordinate ranges (cf. the [list of all options](http://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/chart.html#sage.manifolds.chart.RealChart.plot)):

# In[49]:


Y.plot(X, ranges={r:(1,2), th:(0,pi/2)}, number_values=4, 
       color={r:'blue', th:'green', ph:'red'}, aspect_ratio=1)


# <p>Conversly, the chart $X|_{U}$ can be plotted in terms of the chart $Y$ (this is not possible for the whole chart $X$ since its domain is larger than that of chart $Y$):</p>

# In[50]:


graph = X_U.plot(Y)
show(graph, axes_labels=['r','theta','phi'])


# ## Points on the manifold
# 
# A point on $\mathcal{M}$ is defined by its coordinates in a given chart:

# In[51]:


p = M.point((1,2,-1), chart=X, name='p')
print(p)
p


# <p>Since $X=(\mathcal{M}, (x,y,z))$ is the manifold's default chart, its name can be omitted:</p>

# In[52]:


p = M.point((1,2,-1), name='p')
print(p)
p


# <p>Of course, $p$ belongs to $\mathcal{M}$:</p>

# In[53]:


p in M


# <p>It is also in $U$:</p>

# In[54]:


p in U


# <p>Indeed the coordinates of $p$ have $y\not=0$:</p>

# In[55]:


p.coord(X)


# <p>Note in passing that since $X$ is the default chart on $\mathcal{M}$, its name can be omitted in the arguments of coord():</p>

# In[56]:


p.coord()


# <p>The coordinates of $p$ can also be obtained by letting the chart acting of the point (from the very definition of a chart!):</p>

# In[57]:


X(p)


# <p>Let $q$ be a point with $y = 0$ and $x \geq 0$:</p>

# In[58]:


q = M.point((1,0,2), name='q')


# <p>This time, the point does not belong to $U$:</p>

# In[59]:


q in U


# <p>Accordingly, we cannot ask for the coordinates of $q$ in the chart $Y=(U, (r,\theta,\phi))$:</p>

# In[60]:


try:
    q.coord(Y)
except ValueError as exc:
    print("Error: " + str(exc))


# <p>but we can for point $p$:</p>

# In[61]:


p.coord(Y)


# In[62]:


Y(p)


# <p>Points can be compared:</p>

# In[63]:


q == p


# In[64]:


p1 = U.point((sqrt(3)*sqrt(2), pi-atan(sqrt(5)), atan(2)), chart=Y)
p1 == p


# <p>In SageMath's terminology, points are <strong>elements</strong>, whose <strong>parents</strong> are the manifold on which they have been defined:</p>

# In[65]:


p.parent()


# In[66]:


q.parent()


# In[67]:


p1.parent()


# ## Scalar fields
# 
# A **scalar field** is a differentiable map $U \to \mathbb{R}$, where $U$ is an open subset of $\mathcal{M}$.
# 
# A scalar field is defined by its expressions in terms of charts covering its domain (in general more than one chart is necessary to cover all the domain):

# In[68]:


f = U.scalar_field({X_U: x+y^2+z^3}, name='f')
print(f)


# The coordinate expressions of the scalar field are passed as a Python dictionary, with the charts as keys, hence the writing `{X_U: x+y^2+z^3}`.
# Since in the present case, there is only one chart in the dictionary, an alternative writing is

# In[69]:


f = U.scalar_field(x+y^2+z^3, chart=X_U, name='f')
print(f)


# Since `X_U` is the domain's default chart, it can be omitted in the above declaration:

# In[70]:


f = U.scalar_field(x+y^2+z^3, name='f')
print(f)


# As a mapping $U\subset\mathcal{M}\longrightarrow\mathbb{R}$, a scalar field acts on points, not on coordinates:

# In[71]:


f(p)


# The method `display()` provides the expression of the scalar field in terms of a given chart:

# In[72]:


f.display(X_U)


# If no argument is provided, the method `display()` shows the coordinate expression of the scalar field in all the charts defined on the domain (except for *subcharts*, i.e. the restrictions of some chart to a subdomain):

# In[73]:


f.display()


# Note that the expression of $f$ in terms of the coordinates $(r,\theta,\phi)$ has not been provided by the user but has been automatically computed by means of the change-of-coordinate formula declared above in the transition map.

# In[74]:


f.display(Y)


# In each chart, the scalar field is represented by a function of the chart coordinates (an object of the type `ChartFunction` described above), which is accessible via the method `coord_function()`:

# In[75]:


f.coord_function(X_U)


# In[76]:


f.coord_function(X_U).display()


# In[77]:


f.coord_function(Y)


# In[78]:


f.coord_function(Y).display()


# <p>The "raw" symbolic expression is returned by the method <span style="font-family: courier new,courier;">expr()</span>:</p>

# In[79]:


f.expr(X_U)


# In[80]:


f.expr(Y)


# In[81]:


f.expr(Y) is f.coord_function(Y).expr()


# <p>A scalar field can also be defined by some unspecified function of the coordinates:</p>

# In[82]:


h = U.scalar_field(function('H')(x, y, z), name='h')
print(h)


# In[83]:


h.display()


# In[84]:


h.display(Y)


# In[85]:


h(p) # remember that p is the point of coordinates (1,2,-1) in the chart X_U


# <p>The parent of $f$ is the set $C^\infty(U)$ of all smooth scalar fields on $U$, which is a commutative algebra over $\mathbb{R}$:</p>

# In[86]:


CU = f.parent()
CU


# In[87]:


print(CU)


# In[88]:


CU.category()


# The base ring of the algebra is the field $\mathbb{R}$, which is represented here by SageMath's Symbolic Ring (`SR`):

# In[89]:


CU.base_ring()


# Arithmetic operations on scalar fields are defined through the algebra structure:

# In[90]:


s = f + 2*h
print(s)


# In[91]:


s.display()


# ## Tangent spaces
# 
# The tangent vector space to the manifold at point $p$ is obtained as follows:

# In[92]:


Tp = M.tangent_space(p)
Tp


# In[93]:


print(Tp)


# $T_p\, \mathcal{M}$ is a 2-dimensional vector space over $\mathbb{R}$ (represented here by SageMath's Symbolic Ring (`SR`)):

# In[94]:


print(Tp.category())


# In[95]:


Tp.dim()


# <p>$T_p\, \mathcal{M}$ is automatically endowed with vector bases deduced from the vector frames defined around the point:</p>

# In[96]:


Tp.bases()


# <p>For the tangent space at the point $q$, on the contrary, there is only one pre-defined basis, since $q$ is not in the domain $U$ of the frame associated with coordinates $(r,\theta,\phi)$:</p>

# In[97]:


Tq = M.tangent_space(q)
Tq.bases()


# <p>A random element:</p>

# In[98]:


v = Tp.an_element()
print(v)


# In[99]:


v.display()


# In[100]:


u = Tq.an_element()
print(u)


# In[101]:


u.display()


# <p>Note that, despite what the above simplified writing may suggest (the mention of the point $p$ or $q$ is omitted in the basis vectors), $u$ and $v$ are different vectors, for they belong to different vector spaces:</p>

# In[102]:


v.parent()


# In[103]:


u.parent()


# <p>In particular, it is not possible to add $u$ and $v$:</p>

# In[104]:


try:
    s = u + v
except TypeError as exc:
    print("Error: " + str(exc))


# ## Vector Fields
# 
# Each chart defines a vector frame on the chart domain: the so-called **coordinate basis**:

# In[105]:


X.frame()


# In[106]:


X.frame().domain()  # this frame is defined on the whole manifold


# In[107]:


Y.frame()


# In[108]:


Y.frame().domain() # this frame is defined only on U


# The list of frames defined on a given open subset is returned by the method `frames()`:

# In[109]:


M.frames()


# In[110]:


U.frames()


# In[111]:


M.default_frame()


# Unless otherwise specified (via the command `set_default_frame()`), the default frame is that associated with the default chart:

# In[112]:


M.default_frame() is M.default_chart().frame()


# In[113]:


U.default_frame() is U.default_chart().frame()


# <p>Individual elements of a frame can be accessed by means of their indices:</p>

# In[114]:


e = U.default_frame()
e2 = e[2]
e2


# In[115]:


print(e2)


# <p><span id="cell_outer_109">We may define a new vector field as follows:</span></p>

# In[116]:


v = e[2] + 2*x*e[3]
print(v)


# In[117]:


v.display()


# <p>A vector field can be defined by its components with respect to a given vector frame. When the latter is not specified, the open set's default frame is of course assumed:</p>

# In[118]:


v = U.vector_field(name='v') # vector field defined on the open set U
v[1] = 1+y
v[2] = -x
v[3] = x*y*z
v.display()


# Since version 8.8 of SageMath, it is possible to initialize the components of the vector field while declaring it, so that the above is equivalent to

# In[119]:


v = U.vector_field(1+y, -x, x*y*z, name='v')  # valid only in SageMath 8.8 and higher
v.display() 


# <p>Vector fields on $U$ are Sage <em>element</em> objects, whose <em>parent</em> is the set $\mathfrak{X}(U)$ of vector fields defined on $U$:</p>

# In[120]:


v.parent()


# <p>The set $\mathfrak{X}(U)$ is a module over the commutative algebra $C^\infty(U)$ of scalar fields on $U$:</p>

# In[121]:


print(v.parent())


# In[122]:


print(v.parent().category())


# In[123]:


v.parent().base_ring()


# <p>A vector field acts on scalar fields:</p>

# In[124]:


f.display()


# In[125]:


s = v(f)
print(s)


# In[126]:


s.display()


# In[127]:


e[3].display()


# In[128]:


e[3](f).display()


# <p>Unset components are assumed to be zero:</p>

# In[129]:


w = U.vector_field(name='w')
w[2] = 3
w.display()


# <p>A vector field on $U$ can be expanded in the vector frame associated with the chart $(r,\theta,\phi)$:</p>

# In[130]:


v.display(Y.frame())


# By default, the components are expressed in terms of the default coordinates $(x,y,z)$. To express them in terms of the coordinates $(r,\theta,\phi)$, one should add the corresponding chart as the second argument of the method `display()`:

# In[131]:


v.display(Y.frame(), Y)


# In[132]:


for i in M.irange(): 
    show(e[i].display(Y.frame(), Y))


# The components of a tensor field w.r.t. the default frame can also be obtained as a list, thanks to the operator `[:]`: 

# In[133]:


v[:]


# An alternative is to use the method `display_comp()`:

# In[134]:


v.display_comp()


# To obtain the components w.r.t. another frame, one may go through the method `comp()` and specify the frame:

# In[135]:


v.comp(Y.frame())[:]


# <p>However a shortcut is to provide the frame as the first argument of the square brackets:</p>

# In[136]:


v[Y.frame(), :]


# In[137]:


v.display_comp(Y.frame())


# <p>Components are shown expressed in terms of the default's coordinates; to get them in terms of the coordinates $(r,\theta,\phi)$ instead, add the chart name as the last argument in the square brackets:</p>

# In[138]:


v[Y.frame(), :, Y]


# or specify the chart in `display_comp()`:

# In[139]:


v.display_comp(Y.frame(), chart=Y)


# To get some vector component as a scalar field instead of a coordinate expression, use double square brackets:

# In[140]:


print(v[[1]])


# In[141]:


v[[1]].display()


# In[142]:


v[[1]].expr(X_U)


# <p>A vector field can be defined with components being unspecified functions of the coordinates:</p>

# In[143]:


u = U.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()


# In[144]:


s = v + u
s.set_name('s')
s.display()


# ### Values of vector field at a given point
# 
# The value of a vector field at some point of the manifold is obtained via the method `at()`:

# In[145]:


vp = v.at(p)
print(vp)


# In[146]:


vp.display()


# <p>Indeed, recall that, w.r.t. chart X_U=$(x,y,z)$,  the coordinates of the point $p$ and the components of the vector field $v$ are</p>

# In[147]:


p.coord(X_U)


# In[148]:


v.display(X_U.frame(), X_U)


# Note that to simplify the writing, the symbol used to denote the value of the vector field at point $p$ is the same as that of the vector field itself (namely $v$); this can be changed by the method `set_name()`:

# In[149]:


vp.set_name(latex_name='v|_p')
vp.display()


# <p>Of course, $v|_p$ belongs to the tangent space at $p$:</p>

# In[150]:


vp.parent()


# In[151]:


vp in M.tangent_space(p)


# In[152]:


up = u.at(p)
print(up)


# In[153]:


up.display()


# ## 1-forms
# 
# A **1-form** on $\mathcal{M}$ is a field of linear forms. For instance, it can be the *differential of a scalar field*:

# In[154]:


df = f.differential()
print(df)


# An equivalent writing is

# In[155]:


df = diff(f)


# The method `display()` shows the expansion on the default coframe:

# In[156]:


df.display()


# In the above writing, the 1-form is expanded over the basis $(\mathrm{d}x, \mathrm{d}y, \mathrm{d}z)$ associated with the chart $(x,y,z)$. This basis can be accessed via the method `coframe()`:

# In[157]:


dX = X.coframe()
dX


# The list of all coframes defined on a given manifold open subset is returned by the method `coframes()`:

# In[158]:


M.coframes()


# As for a vector field, the value of the differential form at some point on the manifold is obtained by the method `at()`:

# In[159]:


dfp = df.at(p)
print(dfp)


# In[160]:


dfp.display()


# <p>Recall that</p>

# In[161]:


p.coord()


# <p>The linear form $\mathrm{d}f|_p$ belongs to the dual of the tangent vector space at $p$:</p>

# In[162]:


dfp.parent()


# In[163]:


dfp.parent() is M.tangent_space(p).dual()


# <p>As such, it is acting on vectors at $p$, yielding a real number:</p>

# In[164]:


print(vp)
vp.display()


# In[165]:


dfp(vp)


# In[166]:


print(up)
up.display()


# In[167]:


dfp(up)


# <p>The differential 1-form of the unspecified scalar field $h$:</p>

# In[168]:


dh = h.differential()
dh.display()


# <p>A 1-form <span id="cell_outer_132">can also be defined from scratch:</span></p>

# In[169]:


om = U.one_form(name='omega', latex_name=r'\omega')
print(om)


# <p>It can be specified by providing its components in a given coframe:</p>

# In[170]:


om[:] = [x^2+y^2, z, x-z]    # components in the default coframe (dx,dy,dz)
om.display()


# Since version 8.8 of SageMath, it is possible to initialize the components of the 1-form while declaring it, so that the above is equivalent to

# In[171]:


om = U.one_form(x^2+y^2, z, x-z, name='omega',  # valid only in 
                latex_name=r'\omega')           # SageMath 8.8 and higher
om.display()


# <p>Of course, one may set the components in a frame different from the default one:</p>

# In[172]:


om[Y.frame(), :, Y] = [r*sin(th)*cos(ph), 0, r*sin(th)*sin(ph)]
om.display(Y.frame(), Y)


# <p>The components in the coframe $(\mathrm{d}x,\mathrm{d}y,\mathrm{d}z)$ are updated automatically:</p>

# In[173]:


om.display()


# <p>Let us revert to the values set previously:</p>

# In[174]:


om[:] = [x^2+y^2, z, x-z]
om.display()


# <p>This time, the components in the coframe $(\mathrm{d}r, \mathrm{d}\theta,\mathrm{d}\phi)$ are those that are updated:</p>

# In[175]:


om.display(Y.frame(), Y)


# <p>A 1-form acts on vector fields, resulting in a scalar field:</p>

# In[176]:


print(om(v))
om(v).display()


# In[177]:


print(df(v))
df(v).display()


# In[178]:


om(u).display()


# <p>In the case of a differential 1-form, the following identity holds:</p>

# In[179]:


df(v) == v(f)


# 1-forms are Sage *element* objects, whose *parent* is the $C^\infty(U)$-module $\Omega^{1}(U)$ of all 1-forms defined on $U$:

# In[180]:


df.parent()


# In[181]:


print(df.parent())


# In[182]:


print(om.parent())


# $\Omega^{1}(U)$ is actually the dual of the free module $\mathfrak{X}(U)$:

# In[183]:


df.parent() is v.parent().dual()


# ## Differential forms and exterior calculus
# 
# The **exterior product** of two 1-forms is taken via the method `wedge()` and results in a 2-form:

# In[184]:


a = om.wedge(df)
print(a)
a.display()


# <p>A matrix view of the components:</p>

# In[185]:


a[:]


# <p>Displaying only the non-vanishing components, skipping the redundant ones (i.e. those that can be deduced by antisymmetry):</p>

# In[186]:


a.display_comp(only_nonredundant=True)


# <p>The 2-form $\omega\wedge\mathrm{d}f$ can be expanded on the $(\mathrm{d}r,\mathrm{d}\theta,\mathrm{d}\phi)$ coframe:</p>

# In[187]:


a.display(Y.frame(), Y)


# <p>As a 2-form, $A:=\omega\wedge\mathrm{d}f$ can be applied to a pair of vectors and is antisymmetric:</p>

# In[188]:


a.set_name('A')
print(a(u,v))
a(u,v).display()


# In[189]:


a(u,v) == - a(v,u)


# In[190]:


a.symmetries()


# <p>The <strong>exterior derivative</strong>  of a differential form:</p>

# In[191]:


dom = om.exterior_derivative()
print(dom)
dom.display()


# Instead of invoking the method `exterior_derivative()`, one can use the function `diff()` (available in SageMath 9.2 or higher):

# In[192]:


dom = diff(om)  


# In[193]:


da = diff(a)
print(da)
da.display()


# <p>The exterior derivative is nilpotent:</p>

# In[194]:


ddf = diff(df)
ddf.display()


# In[195]:


ddom = diff(dom)
ddom.display()


# ## Lie derivative
# 
# The Lie derivative of any tensor field with respect to a vector field is computed by the method `lie_derivative()`, with the vector field as the argument:

# In[196]:


lv_om = om.lie_derivative(v)
print(lv_om)
lv_om.display()


# In[197]:


lu_dh = dh.lie_derivative(u)
print(lu_dh)
lu_dh.display()


# <p>Let us check <strong>Cartan identity</strong> on the 1-form $\omega$:</p>
# <p>$\mathcal{L}_v \omega = v\cdot \mathrm{d}\omega + \mathrm{d}\langle \omega, v\rangle$</p>
# <p>and on the 2-form $A$:</p>
# <p>$\mathcal{L}_v A = v\cdot \mathrm{d}A + \mathrm{d}(v\cdot A)$</p>

# In[198]:


om.lie_derivative(v) == v.contract(diff(om)) + diff(om(v))


# In[199]:


a.lie_derivative(v) == v.contract(diff(a)) + diff(v.contract(a))


# <p>The Lie derivative of a vector field along another one is the <strong>commutator</strong> of the two vectors fields:</p>

# In[200]:


v.lie_derivative(u)(f) == u(v(f)) - v(u(f))


# ## Tensor fields of arbitrary rank
# 
# Up to now, we have encountered tensor fields
# 
# - of type (0,0) (i.e. scalar fields),
# - of type (1,0) (i.e. vector fields),
# - of type (0,1) (i.e. 1-forms),
# - of type (0,2) and antisymmetric (i.e. 2-forms).
# 
# More generally, tensor fields of any type $(p,q)$ can be introduced in SageMath. For instance a tensor field of type (1,2) on the open subset $U$ is declared as follows:

# In[201]:


t = U.tensor_field(1, 2, name='T')
print(t)


# <p>As for vectors or 1-forms, the tensor's components with respect to the domain's default frame are set by means of square brackets:</p>

# In[202]:


t[1,2,1] = 1 + x^2
t[3,2,1] = x*y*z


# <p>Unset components are zero:</p>

# In[203]:


t.display()


# In[204]:


t[:]


# <p>Display of the nonzero components:</p>

# In[205]:


t.display_comp()


# <p>Double square brackets return the component (still w.r.t. the default frame) as a scalar field, while single square brackets return the expression of this scalar field in terms of the domain's default coordinates:</p>

# In[206]:


print(t[[1,2,1]])
t[[1,2,1]].display()


# In[207]:


print(t[1,2,1])
t[1,2,1]


# <p>A tensor field of type (1,2) maps a 3-tuple (1-form, vector field, vector field) to a scalar field:</p>

# In[208]:


print(t(om, u, v))
t(om, u, v).display()


# <p>As for vectors and differential forms, the tensor components can be taken in any frame defined on the manifold:</p>

# In[209]:


t[Y.frame(), 1,1,1, Y]


# ## Tensor calculus
# 
# The **tensor product** $\otimes$ is denoted by `*`:

# In[210]:


print(v.tensor_type())
print(a.tensor_type())


# In[211]:


b = v*a
print(b)
b


# <p>The tensor product preserves the (anti)symmetries: since $A$ is a 2-form, it is antisymmetric with respect to its two arguments (positions 0 and 1); as a result, b is antisymmetric with respect to its last two arguments (positions 1 and 2):</p>

# In[212]:


a.symmetries()


# In[213]:


b.symmetries()


# <p>Standard<strong> tensor arithmetics</strong> is implemented:</p>

# In[214]:


s = - t + 2*f* b
print(s)


# **Tensor contractions** are dealt with by the methods`trace()` and `contract()`: for instance, let us contract the tensor $T$ w.r.t. its first two arguments (positions 0 and 1), i.e. let us form the tensor $c$ of components $c_i = T^k_{\ \, k i}$:

# In[215]:


c = t.trace(0,1)
print(c)


# An alternative to the writing `trace(0,1)` is to use the **index notation** to denote the contraction: the indices are given in a string inside the `[]` operator, with `'^'` in front of the contravariant indices and `'_'` in front of the covariant ones:

# In[216]:


c1 = t['^k_ki']
print(c1)
c1 == c


# The contraction is performed on the repeated index (here `k`); the letter denoting the remaining index (here `i`) is arbitrary:

# In[217]:


t['^k_kj'] == c


# In[218]:


t['^b_ba'] == c


# <p>It can even be replaced by a dot:</p>

# In[219]:


t['^k_k.'] == c


# <p>LaTeX notations are allowed:</p>

# In[220]:


t['^{k}_{ki}'] == c


# as well as Greek letters (only for SageMath 9.2 or higher):

# In[221]:


t['^μ_μα'] == c


# <p>The contraction $T^i_{\ j k} v^k$ of the tensor fields $T$ and $v$ is taken as follows (2 refers to the last index position of $T$ and 0 to the only index position of v):</p>

# In[222]:


tv = t.contract(2, v, 0)
print(tv)


# <p>Since 2 corresponds to the last index position of $T$ and 0 to the first index position of $v$, a shortcut for the above is</p>

# In[223]:


tv1 = t.contract(v)
print(tv1)


# In[224]:


tv1 == tv


# <p>Instead of <span style="font-family: courier new,courier;">contract()</span>, the <strong>index notation</strong>, combined with the <strong>*</strong> operator, can be used to denote the contraction:</p>

# In[225]:


t['^i_jk']*v['^k'] == tv


# <p>The non-repeated indices can be replaced by dots:</p>

# In[226]:


t['^._.k']*v['^k'] == tv


# ## Metric structures
# 
# A **Riemannian metric** on the manifold $\mathcal{M}$ is declared as follows:

# In[227]:


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


# <p>It is a symmetric tensor field of type (0,2):</p>

# In[228]:


g.parent()


# In[229]:


print(g.parent())


# In[230]:


g.symmetries()


# <p>The metric is initialized by its components with respect to some vector frame. For instance, using the default frame of $\mathcal{M}$:</p>

# In[231]:


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


# <p>The components w.r.t. another vector frame are obtained as for any tensor field:</p>

# In[232]:


g.display(Y.frame(), Y)


# <p>Of course, the metric acts on vector pairs:</p>

# In[233]:


print(g(u,v))
g(u,v).display()


# <p>The <strong>Levi-Civita connection</strong> associated to the metric $g$:</p>

# In[234]:


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


# <p>The Christoffel symbols with respect to the manifold's default coordinates:</p>

# In[235]:


nabla.coef()[:]


# <p>The Christoffel symbols with respect to the coordinates $(r,\theta,\phi)$:</p>

# In[236]:


nabla.coef(Y.frame())[:, Y]


# A nice view is obtained via the method `display()` (by default, only the nonzero connection coefficients are shown):

# In[237]:


nabla.display(frame=Y.frame(), chart=Y)


# One may also use the method `christoffel_symbols_display()` of the metric, which (by default) displays only the non-redundant Christoffel symbols:

# In[238]:


g.christoffel_symbols_display(Y)


# <p>The connection acting as a covariant derivative:</p>

# In[239]:


nab_v = nabla(v)
print(nab_v)
nab_v.display()


# <p>Being a Levi-Civita connection, $\nabla_g$ is torsion.free:</p>

# In[240]:


print(nabla.torsion())
nabla.torsion().display()


# <p>In the present case, it is also flat:</p>

# In[241]:


print(nabla.riemann())
nabla.riemann().display()


# <p>Let us consider a non-flat metric, by changing $g_{rr}$ to $1/(1+r^2)$:</p>

# In[242]:


g[Y.frame(), 1,1, Y] = 1/(1+r^2)
g.display(Y.frame(), Y)


# <p>For convenience, we change the default chart on the domain $U$ to Y=$(U,(r,\theta,\phi))$:</p>

# In[243]:


U.set_default_chart(Y)


# <p>In this way, we do not have to specify Y when asking for coordinate expressions in terms of $(r,\theta,\phi)$:</p>

# In[244]:


g.display(Y.frame())


# <p>We recognize the metric of the hyperbolic space $\mathbb{H}^3$. Its expression in terms of the chart $(U,(x,y,z))$ is</p>

# In[245]:


g.display(X_U.frame(), X_U)


# <p>A matrix view of the components may be more appropriate:</p>

# In[246]:


g[X_U.frame(), :, X_U]


# <p>We extend these components, a priori defined only on $U$, to the whole manifold $\mathcal{M}$, by demanding the same coordinate expressions in the frame associated to the chart X=$(\mathcal{M},(x,y,z))$:</p>

# In[247]:


g.add_comp_by_continuation(X.frame(), U, X)
g.display()


# <p>The Levi-Civita connection is automatically recomputed, after the change in $g$:</p>

# In[248]:


nabla = g.connection()


# <p>In particular, the Christoffel symbols are different:</p>

# In[249]:


nabla.display(only_nonredundant=True)


# In[250]:


nabla.display(frame=Y.frame(), chart=Y, only_nonredundant=True)


# <p>The <strong>Riemann tensor</strong> is now</p>

# In[251]:


Riem = nabla.riemann()
print(Riem)
Riem.display(Y.frame())


# <p>Note that it can be accessed directely via the metric, without any explicit mention of the connection:</p>

# In[252]:


g.riemann() is nabla.riemann()


# <p>The <strong>Ricci tensor</strong> is</p>

# In[253]:


Ric = g.ricci()
print(Ric)
Ric.display(Y.frame())


# <p>The<strong> Weyl tensor</strong> is:</p>

# In[254]:


C = g.weyl()
print(C)
C.display()


# <p>The Weyl tensor vanishes identically because the dimension of $\mathcal{M}$ is 3.</p>
# <p>Finally, the <strong>Ricci scalar</strong> is </p>

# In[255]:


R = g.ricci_scalar()
print(R)
R.display()


# We recover the fact that $\mathbb{H}^3$ is a Riemannian manifold of constant negative curvature.

# ## Tensor transformations induced by a metric
# 
# The most important tensor transformation induced by the metric $g$ is the so-called **musical isomorphism**, or **index raising** and **index lowering**:

# In[256]:


print(t)


# In[257]:


t.display()


# In[258]:


t.display(X_U.frame(), X_U)


# Raising the last index (position 2) of $T$ with $g$:

# In[259]:


s = t.up(g, 2)
print(s)


# Note that the raised index becomes the *last* one among the contravariant indices,  i.e. the tensor $s$ returned by the method `up` is 
# $$ s^{ab}_{\ \ \ \, c} = g^{bi} T^a_{\ \ ic}$$

# See the [up() documentation](https://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/differentiable/tensorfield.html#sage.manifolds.differentiable.tensorfield.TensorField.up) for more details.

# <p>Raising all the covariant indices of $T$ (i.e. those at the positions 1 and 2):</p>

# In[260]:


s = t.up(g)
print(s)


# Lowering all contravariant indices of $T$ (i.e. the index at position 0):

# In[261]:


s = t.down(g)
print(s)


# Note that the lowered index becomes the *first* one among the covariant indices,  i.e. the tensor $s$ returned by the method `down` is 
# $$ s_{abc} = g_{ai} T^i_{\ \ bc}$$

# See the [down() documentation](https://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/differentiable/tensorfield.html#sage.manifolds.differentiable.tensorfield.TensorField.down) for more details.

# ## Hodge duality
# 
# The volume 3-form (Levi-Civita tensor) associated with the metric $g$ is

# In[262]:


epsilon = g.volume_form()
print(epsilon)
epsilon.display()


# In[263]:


epsilon.display(Y.frame())


# In[264]:


print(f)
f.display()


# In[265]:


sf = f.hodge_dual(g.restrict(U))
print(sf)
sf.display()


# <p>We check the classical formula $\star f = f\, \epsilon_g$, or, more precisely, $\star f = f\, \epsilon_g|_U$ (for $f$ is defined on $U$ only):</p>

# In[266]:


sf == f * epsilon.restrict(U)


# <p>The Hodge dual of a 1-form is a 2-form:</p>

# In[267]:


print(om)
om.display()


# In[268]:


som = om.hodge_dual(g)
print(som)
som.display()


# <p>The Hodge dual of a 2-form is a 1-form:</p>

# In[269]:


print(a)


# In[270]:


sa = a.hodge_dual(g)
print(sa)
sa.display()


# <p>Finally, the Hodge dual of a 3-form is a 0-form:</p>

# In[271]:


print(da)
da.display()


# In[272]:


sda = da.hodge_dual(g)
print(sda)
sda.display()


# <p>In dimension 3 and for a Riemannian metric, the Hodge star is idempotent:</p>

# In[273]:


sf.hodge_dual(g) == f


# In[274]:


som.hodge_dual(g) == om


# In[275]:


sa.hodge_dual(g) == a


# In[276]:


sda.hodge_dual(g.restrict(U)) == da


# ## Getting help
# 
# To get the list of functions (methods) that can be called on a object, type the name of the object, followed by a dot and the TAB key, e.g. `sa.<TAB>`.
# 
# To get information on an object or a method, use the question mark:

# In[277]:


get_ipython().run_line_magic('pinfo', 'nabla')


# In[278]:


get_ipython().run_line_magic('pinfo', 'g.ricci_scalar')


# Using a double question mark leads directly to the **Python source code** (SageMath is **open source**, isn't it?):

# In[279]:


get_ipython().run_line_magic('pinfo2', 'g.ricci_scalar')


# ## Going further
# 
# Have a look at the [examples on SageManifolds page](https://sagemanifolds.obspm.fr/examples.html), especially the 
# [2-dimensional sphere](https://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Notebooks/SM_sphere_S2.ipynb) for usage on a non-parallelizable manifold (each scalar field has to be defined in at least two coordinate charts, the $C^\infty(\mathcal{M})$-module $\mathfrak{X}(\mathcal{M})$ is no longer free and each tensor field has to be defined in at least two vector frames).
# 
# You may also take a look at the [tutorial videos](https://www.youtube.com/playlist?list=PLnrOCYZpQUuJlnQbQ48zgGk-Ks1t145Yw) by Christian Bär.