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

# # Smooth manifolds, vector fields and tensor fields
# 
# This notebook accompanies the lecture
# [Symbolic tensor calculus on manifolds](http://sagemanifolds.obspm.fr/jncf2018/) at JNCF 2018.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/JNCF2018/jncf18_vector.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter Notebook server, via the command `sage -n jupyter`

# In[1]:


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


# Let us restore the 2-sphere manifold constructed in the [preceeding worksheet](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/JNCF2018/jncf18_scalar.ipynb):

# In[2]:


M = Manifold(2, 'M')
U = M.open_subset('U')
XU.<x,y> = U.chart()
V = M.open_subset('V')
XV.<xp,yp> = V.chart("xp:x' yp:y'")
M.declare_union(U,V)
XU_to_XV = XU.transition_map(XV, 
                             (x/(x^2+y^2), y/(x^2+y^2)), 
                             intersection_name='W',
                             restrictions1= x^2+y^2!=0, 
                             restrictions2= xp^2+yp^2!=0)
XV_to_XU = XU_to_XV.inverse()
M.atlas()


# We also reconstruct the point $p$:

# In[3]:


p = U((1,2), chart=XU, name='p')
print(p)


# and the embedding $\mathbb{S}^2 \to \mathbb{R}^3$:

# In[4]:


R3 = Manifold(3, 'R^3', r'\mathbb{R}^3')
XR3.<X,Y,Z> = R3.chart()
Phi = M.diff_map(R3, {(XU, XR3): 
                       [2*x/(1+x^2+y^2), 2*y/(1+x^2+y^2),
                        (x^2+y^2-1)/(1+x^2+y^2)],
                      (XV, XR3): 
                       [2*xp/(1+xp^2+yp^2), 2*yp/(1+xp^2+yp^2),
                        (1-xp^2-yp^2)/(1+xp^2+yp^2)]},
                  name='Phi', latex_name=r'\Phi')
Phi.display()


# In[5]:


graph = XU.plot(chart=XR3, mapping=Phi, number_values=25, 
                label_axes=False) + \
        XV.plot(chart=XR3, mapping=Phi, number_values=25, 
                color='green', label_axes=False) + \
        p.plot(chart=XR3, mapping=Phi, label_offset=0.05)
show(graph, viewer='threejs', online=True)


# Finally we reconstruct the scalar field $f$:

# In[6]:


f = M.scalar_field({XU: 1/(1+x^2+y^2), XV: (xp^2+yp^2)/(1+xp^2+yp^2)},
                   name='f')
f.display()


# and assign the Python variable `CM` to the algebra of scalar fields:

# In[7]:


CM = M.scalar_field_algebra()
CM


# ## Tangent vectors
# 
# The tangent space at the point $p$ introduced above is generated by

# In[8]:


Tp = M.tangent_space(p)
Tp


# It is a vector space over $\mathbb{R}$, which is represented by Sage's Symbolic Ring:

# In[9]:


print(Tp.category())


# The dimension of $T_p M$ is the same as that of $M$:

# In[10]:


dim(Tp)


# Tangent spaces are implemented as a class inherited from `TangentSpace` via the category framework:

# In[11]:


type(Tp)


# The class `TangentSpace` actually inherits from the generic class
# `FiniteRankFreeModule`, which, in SageMath, is devoted to free modules of finite rank
# without any distinguished basis:

# In[12]:


isinstance(Tp, FiniteRankFreeModule)


# Two bases of $T_p M$ are already available: those generated by the derivations
# at $p$ along the coordinates of charts `XU` and `XV` respectively:

# In[13]:


Tp.bases()


# None of these bases is distinguished, but one if the default one, which
# simply means that it is the basis to be considered if the basis argument
# is skipped in some methods:

# In[14]:


Tp.default_basis()


# A tangent vector is created as an element of the tangent space by the
# standard SageMath procedure
# `new\_element = parent(...)`, where `...`
# stands for some material sufficient to construct the element:

# In[15]:


vp = Tp((-3, 2), name='v')
print(vp)


# Since the basis is not specified, the pair $(-3,2)$ refers to components
# with respect to the default basis:

# In[16]:


vp.display()


# We have of course

# In[17]:


vp.parent()


# In[18]:


vp in Tp


# As other manifold objects, tangent vectors have some plotting capabilities:

# In[19]:


graph += vp.plot(chart=XR3, mapping=Phi, scale=0.5, color='gold')
show(graph, viewer='threejs', online=True)


# The main attribute of the object `vp` representing the vector $v$ is the dictionary `_components`, which stores the components of $v$ in various bases of $T_p M$:

# In[20]:


vp._components


# As we can see above, the keys of the dictionary `_components` are the bases of $T_p M$, while the values belongs to the class [Components](http://doc.sagemath.org/html/en/reference/tensor_free_modules/sage/tensor/modules/comp.html) devoted to store ring elements indexed by integers or tuples of integers:

# In[21]:


vpc = vp._components[Tp.default_basis()]
vpc


# In[22]:


type(vpc)


# The components themselves are stored in the dictionary `_comp` of the `Components` object, with the indices as keys:

# In[23]:


vpc._comp


# ## Module of vector fields
# 
# The $C^\infty(M)$-module of vector fields on $M$, $\mathfrak{X}(M)$, is obtained as

# In[24]:


YM = M.vector_field_module()
YM


# In[25]:


YM.category()


# In[26]:


YM.base_ring() is CM


# $\mathfrak{X}(M)$ is not a free module (at least its SageMath implementation does not
# belong to the class `FiniteRankFreeModule`):

# In[27]:


isinstance(YM, FiniteRankFreeModule)


# This is because $M=\mathbb{S}^2$ is not a parallelizable manifold:

# In[28]:


M.is_manifestly_parallelizable()


# Via the category mechanism,
# the module $\mathfrak{X}(M)$ is implemented by a dynamically-generated subclass
# of the class `VectorFieldModule`, which is devoted to modules of vector fields
# on non-parallelizable manifolds:

# In[29]:


type(YM)


# On the contrary, the set $\mathfrak{X}(U)$ of vector fields on $U$ is a free module of finite rank over the algebra $C^\infty(U)$:

# In[30]:


YU = U.vector_field_module()
isinstance(YU, FiniteRankFreeModule)


# In[31]:


YU.base_ring()


# This is because the open subset $U$ is a parallelizable manifold:

# In[32]:


U.is_manifestly_parallelizable()


# being a coordinate chart domain:

# In[33]:


U.is_manifestly_coordinate_domain()


# We can check that in the atlas of $U$, at least one chart has $U$ for domain:

# In[34]:


U.atlas()


# The rank of $\mathfrak{X}(U)$ is the manifold's dimension:

# In[35]:


rank(YU)


# Via the category mechanism,
# the free module $\mathfrak{X}(U)$ is implemented by a dynamically-generated subclass
# of the class `VectorFieldFreeModule`, which is devoted to modules of vector fields
# on parallelizable manifolds:

# In[36]:


type(YU)


# The class `VectorFieldFreeModule` is itself a subclass
# of the generic class `FiniteRankFreeModule`:

# In[37]:


class_graph(
sage.manifolds.differentiable.vectorfield_module.VectorFieldFreeModule
).plot()


# Since $U$ is a chart domain, the free module $\mathfrak{X}(U)$ is automatically endowed with a basis: the coordinate frame associated to the chart:

# In[38]:


YU.bases()


# Let us denote by `eU` this frame. We can set `eU = YU.bases()[0]` or
# alternatively

# In[39]:


eU = YU.default_basis()
eU


# Another equivalent instruction would have been `eU = U.default_frame()`.
# 
# Similarly, $\mathfrak{X}(V)$ is a free module, endowed with the coordinate frame
# associated to stereographic coordinates from the South pole, which we
# denote by `eV`:

# In[40]:


YV = V.vector_field_module()
YV.bases()


# In[41]:


eV = YV.default_basis()
eV


# If we consider the intersection $W=U\cap V$, we notice its module
# of vector fields is endowed with two bases, reflecting the fact that
# $W$ is covered by two charts: $(W,(x,y))$ and $(W,(x',y'))$:

# In[42]:


W = U.intersection(V)
YW = W.vector_field_module()
YW.bases()


# Let us denote by `eUW` and `eUV` these two bases, which are
# actually the restrictions of the vector frames `eU` and `eV` to $W$:

# In[43]:


eUW = eU.restrict(W)
eVW = eV.restrict(W)
YW.bases() == [eUW, eVW]


# The free module $\mathfrak{X}(W)$ is also automatically endowed with automorphisms
# connecting the two bases, i.e. change-of-frame operators:

# In[44]:


W.changes_of_frame()


# The first of them is

# In[45]:


P = W.change_of_frame(eUW, eVW)
P


# It belongs to the general linear group of the free module $\mathfrak{X}(W)$:

# In[46]:


P.parent()


# and its matrix is deduced from the Jacobian matrix of the transition map `XV` $\to$ `XU`:

# In[47]:


P[:]


# ## An example of vector field
# 
# We introduce a vector field $v$ on $M$ by

# In[48]:


v = M.vector_field(name='v')
v[eU, 0] = f.restrict(U)
v[eU, 1] = -2
v.display(eU)


# Notice that at this stage, we have defined $v$ only on $U$, by setting
# its components in the vector frame `eU`, either explicitely as scalar
# fields, like the component $v^0$ set to the restriction of $f$ to $U$ or
# to some symbolic expression, like for the component $v^1$: the $-2$
# will be coerced to the constant scalar field of value $-2$.
# We can ask for the scalar-field value of a components via the double-bracket
# operator, since `eU` is the default frame on $M$, we don't have to specify
# it:

# In[49]:


v[[0]]


# In[50]:


v[[0]].display()


# Note that the single bracket operator returns a chart function
# of the component:

# In[51]:


v[0]


# The restriction of $v$ to $W$ is of course

# In[52]:


v.restrict(W).display(eUW)


# Since we have a second vector frame on $W$, namely `eVW`, and the
# change-of-frame automorphisms are known, we can ask for the components
# of $v$ with respect to that frame:

# In[53]:


v.restrict(W).display(eVW)


# Notice that the components are expressed in terms of the coordinates $(x,y)$
# since they form the default chart on $W$. To have them expressed in
# terms of the coordinates $(x',y')$, we have to add the restriction of
# the chart
# $(V,(x',y'))$ to $W$ as the second argument of the method `display()`:

# In[54]:


v.restrict(W).display(eVW, XV.restrict(W))


# We extend the expression of $v$ to the full vector frame `XV`
# by continuation of this expression:

# In[55]:


v.add_comp_by_continuation(eV, W, chart=XV)


# We have then

# In[56]:


v.display(eV)


# At this stage, the vector field $v$ is defined in all $M$.
# According to the hairy ball theorem\index{hairy ball theorem}, it has to vanish somewhere.
# Let us show that this occurs at the North pole, by first introducing the
# latter, as the point of stereographic coordinates $(x',y')=(0,0)$:

# In[57]:


N = M((0,0), chart=XV, name='N')
print(N)


# As a check, we verify that the image of $N$ by the canonical embedding
# $\Phi: \mathbb{S}^2 \to \mathbb{R}^3$ is the point of Cartesian coordinates $(0,0,1)$:

# In[58]:


XR3(Phi(N))


# The vanishing of $\left. v\right| _N$:

# In[59]:


v.at(N).display()


# On the other hand, $v$ does not vanish at the point $p$ introduced above:

# In[60]:


v.at(p).display()


# We may plot the vector field $v$ in terms of the stereographic coordinates
# from the North pole:

# In[61]:


v.plot(chart=XU, chart_domain=XU, max_range=2, 
       number_values=5, scale=0.4, aspect_ratio=1)


# or in term of those from the South pole:

# In[62]:


v.plot(chart=XV, chart_domain=XV, max_range=2, 
       number_values=9, scale=0.05, aspect_ratio=1)


# Thanks to the embedding $\Phi$, we may also have a 3D plot of $v$
# atop of the 3D plot already obtained:

# In[63]:


graph_v = v.plot(chart=XR3, mapping=Phi, chart_domain=XU, 
                 number_values=7, scale=0.2) + \
          v.plot(chart=XR3, mapping=Phi, chart_domain=XV, 
                 number_values=7, scale=0.2)
show(graph + graph_v, viewer='threejs', online=True)


# Note that the sampling, performed on the two charts `XU` and `XV`
# is not uniform on the sphere. A better sampling would be acheived by introducing
# spherical coordinates.

# ### Some details about the implementation of vector fields

# Vector fields on $M$ are implemented via the class `VectorField`:

# In[64]:


isinstance(v, sage.manifolds.differentiable.vectorfield.VectorField)


# Since $M$ is not parallelizable, the representation of the vector field
# $v$ via its restrictions to parallelizable open subsets; they are stored in the dictionary `_restrictions`, whose keys are the open subsets:

# In[65]:


v._restrictions


# Let us consider one of these restrictions, for instance the restriction to $U$:

# In[66]:


vU = v._restrictions[U]
vU is v.restrict(U)


# Since $U$ is a parallelizable open subset, `vU` belongs to the class
# `VectorFieldParal`:

# In[67]:


isinstance(vU, sage.manifolds.differentiable.vectorfield.VectorFieldParal)


# The main attribute of `vU` is the dictionary `_components`, which stores the components of `vU` with respect to (possibly various) vector frames on $U$, the keys of that dictionary being the vector frames:

# In[68]:


vU._components


# In[69]:


v._restrictions[W]._components


# The values of the dictionary `_components` belong to the same class [Components](http://doc.sagemath.org/html/en/reference/tensor_free_modules/sage/tensor/modules/comp.html) as that used for the storage of components of tangent vectors (cf. the example of `vp` above):

# In[70]:


vUc = vU._components[eU]
vUc


# In[71]:


type(vUc)


# As already mentioned above, the components themselves are stored in a dictionary, `_comp`, whose keys are the indices:

# In[72]:


vUc._comp


# The difference with the tangent vector case is that the values are now scalar fields, i.e. elements of $C^\infty(U)$ in the present case. This is of course in agreement with the treatment of $\mathfrak{X}(U)$ as a free module over $C^\infty(U)$.

# Let us perform some algebraic operation on $v$:

# In[73]:


w = v + f*v
w


# The code for the addition is accessible via

# In[74]:


get_ipython().run_line_magic('pinfo2', 'v.__add__')


# As for the addition of scalar field (cf. [this notebook](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/JNCF2018/jncf18_scalar.ipynb)), we see that `__add__()` is implemented at the level
# of the generic class `Element` from which `VectorField` inherits.
# When both operands of the addition have the same parent, as here, `__add__()` invokes the
# method `_add_()`
# (note the single underscore on each side of `add`). This operator is
# implemented at the level of `TensorField`, as it can be checked from the source code:

# In[75]:


get_ipython().run_line_magic('pinfo2', 'v._add_')


# The first step in the addition of two vector fields is to search in the
# restrictions of both vector fields for common domains: this is performed via the method `_common_subdomains`. Then the addition is performed at the level of the restrictions. The rest of the code is simply the set up of the vector field object containing the result.
# Recursively, the addition performed will reach a level at which
# the domains are parallelizable. Then a different method `_add_()`, will
# be involved, as we can check on `vU`:

# In[76]:


get_ipython().run_line_magic('pinfo2', 'vU._add_')


# We see that this method `_add_()` is implemented
# at the level of tensors on free modules, i.e. in the class
# [FreeModuleTensor](http://doc.sagemath.org/html/en/reference/tensor_free_modules/sage/tensor/modules/free_module_tensor.html), from which `VectorFieldParal` inherits. Here the free module is clearly $\mathfrak{X}(U)$. The addition amounts to adding the components in a
# a basis of the free module in which both operands have known components. Such
# a basis is returned by the method `common_basis` invoked in the first line.
# If necessary, this method can use change-of-basis formulas to compute the
# components of `self` or `other` in a common basis.
# The addition of the components in the returned basis involves the method `__add__()` of 
# class `Components`; we can examine the corresponding code via the `Components` object `vUc`:

# In[77]:


get_ipython().run_line_magic('pinfo2', 'vUc.__add__')


# We note that the computation can be parallelized on the
# components if the user has turned on parallelization (this is done with the command
# `Parallelism().set(nproc=8)` (for 8 threads)). Focusing on the sequential code, we see that the addition is
# performed component by component. Each component being an element of
# $C^\infty(U)$ (the base ring of $\mathfrak{X}(U)$), this addition is that
# of scalar fields, as presented in [this notebook](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/JNCF2018/jncf18_scalar.ipynb).
# 

# ### Action of a vector field on a scalar field
# 
# The action of $v$ on $f$ is defined pointwise by
# considering $v$ at each point $p\in M$ as a derivation (the very definition of a tangent vector); the result is then a
# scalar field $v(f)$ on $M$:

# In[78]:


vf = v(f)
vf


# In[79]:


vf.display()


# # Tensor fields
# 
# ### 1-forms
# 
# Let us start with a 1-form, namely the differential of $f$:

# In[80]:


df = f.differential()
df


# In[81]:


print(df)


# A 1-form is a tensor field of type $(0,1)$:

# In[82]:


df.tensor_type()


# while a vector field is a tensor field of type $(1,0)$:

# In[83]:


v.tensor_type()


# Specific 1-forms are those forming the dual basis (coframe) of a given vector frame: for instance for the vector frame `eU` = $(U, (\partial_{x}, \partial_{y}))$, we have

# In[84]:


eU.dual_basis()


# In[85]:


print(eU.dual_basis()[0])


# Since `eU` is the default frame on $M$, the default display of $\mathrm{d}f$ is performed in terms of `eU`'s coframe:

# In[86]:


df.display()


# In[87]:


df[0] == diff(f.expr(), x)


# In[88]:


df[1] == diff(f.expr(), y)


# In[89]:


df.display(eV)


# In[90]:


df[eV,0,XV]


# In[91]:


df[eV,0,XV] == diff(f.expr(XV), xp)


# In[92]:


df[eV,1,XV] == diff(f.expr(XV), yp)


# The parent of $\mathrm{d}f$ is the set $\Omega^1(M)$ of all 1-forms on $M$, considered as a $C^\infty(M)$-module:

# In[93]:


print(df.parent())
df.parent()


# In[94]:


df.parent().base_ring()


# This module is actually the dual of `YM` = $\mathfrak{X}(M)$:

# In[95]:


YM.dual()


# A 1-form acts on vector fields, yielding a scalar field:

# In[96]:


print(df(v))


# This scalar field is nothing but the result of the action of $v$ on $f$,
# considering $v$ at each point $p\in M$ as a derivation (the very definition of a tangent vector):

# In[97]:


df(v) == v(f)


# ### More general tensor fields

# We construct a tensor of type $(1,1)$ by taking the tensor product $v\otimes \mathrm{d}f$:

# In[98]:


t = v * df
t


# In[99]:


t.display()


# In[100]:


t.display(eV)


# In[101]:


t.display_comp()


# The parent of `t` is the set $\mathcal{T}^{(1,1)}(M)$ of all type-$(1,1)$ tensor fields on $M$,
# considered as a $C^\infty(M)$-module:

# In[102]:


print(t.parent())
t.parent()


# In[103]:


t.parent().base_ring()


# In[104]:


t._restrictions


# As for vector fields, since $M$ is not parallelizable, the $C^\infty(M)$-module
# $\mathcal{T}^{(1,1)}(M)$ is not free and the tensor fields are described by their restrictions to parallelizable subdomains:

# In[105]:


print(t._restrictions[U].parent())


# In[106]:


t._restrictions[U].parent().base_ring()


# ### Riemannian metric on $M$
# 
# The standard metric on $M=\mathbb{S}^2$ is that induced by the Euclidean metric of $\mathbb{R}^3$. Let us start by defining the latter:

# In[107]:


h = R3.metric('h')
h[0,0], h[1,1], h[2, 2] = 1, 1, 1
h.display()


# The metric $g$ on $M$ is the pullback of $h$ associated with the embedding $\Phi$:

# In[108]:


g = M.metric('g')
g.set( Phi.pullback(h) )
print(g)


# Note that we could have defined $g$ intrinsically, i.e. by providing its components in the two frames `eU` and `eV`, as we did for the metric $h$ on $\mathbb{R}^3$. Instead, we have chosen to get it as the pullback of $h$, as an example of pullback associated with some differential map. 
# 
# The metric is a symmetric tensor field of type (0,2):

# In[109]:


g.tensor_type()


# The expression of the metric in terms of the default frame on $M$ (`eU`):

# In[110]:


g.display()


# We may factorize the metric components:

# In[111]:


g[0,0].factor() ; g[1,1].factor()


# In[112]:


g.display()


# A matrix view of the components of $g$ in the manifold's default frame:

# In[113]:


g[:]


# Display in terms of the vector frame $(V, (\partial_{x'}, \partial_{y'}))$:

# In[114]:


g.display(eV)


# The metric acts on vector field pairs, resulting in a scalar field:

# In[115]:


print(g(v,v))


# In[116]:


g(v,v).parent()


# In[117]:


g(v,v).display()


# The **Levi-Civita connection** associated with the metric $g$:

# In[118]:


nab = g.connection()
print(nab)
nab


# The nonzero Christoffel symbols of $g$ (skipping those that can be deduced by symmetry on the last two indices) w.r.t. the chart `XU`:

# In[119]:


g.christoffel_symbols_display(chart=XU)


# $\nabla_g$ acting on the vector field $v$:

# In[120]:


Dv = nab(v)
print(Dv)


# In[121]:


Dv.display()


# ### Curvature
# 
# The Riemann tensor associated with the metric $g$:

# In[122]:


Riem = g.riemann()
print(Riem)
Riem.display()


# The components of the Riemann tensor in the default frame on $M$:

# In[123]:


Riem.display_comp()


# In[124]:


print(Riem.parent())


# In[125]:


Riem.symmetries()


# The Riemann tensor associated with the Euclidean metric $h$ on $\mathbb{R}^3$ is identically zero:

# In[126]:


h.riemann().display()


# The Ricci tensor and the Ricci scalar:

# In[127]:


Ric = g.ricci()
Ric.display()


# In[128]:


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


# Hence we recover the fact that $(\mathbb{S}^2,g)$ is a Riemannian manifold of constant positive curvature.
# 
# In dimension 2, the Riemann curvature tensor is entirely determined by the Ricci scalar $R$ according to
# $$ R^i_{\ \, jlk} = \frac{R}{2} \left( \delta^i_{\ \, k} g_{jl} - \delta^i_{\ \, l} g_{jk} \right)$$
# Let us check this formula here, under the form $R^i_{\ \, jlk} = -R g_{j[k} \delta^i_{\ \, l]}$:

# In[129]:


delta = M.tangent_identity_field()
Riem == - R*(g*delta).antisymmetrize(2,3)


# Similarly the relation $\mathrm{Ric} = (R/2)\; g$ must hold:

# In[130]:


Ric == (R/2)*g


# The **Levi-Civita tensor** associated with $g$:

# In[131]:


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


# The exterior derivative of the 2-form $\epsilon_g$:

# In[132]:


print(eps.exterior_derivative())


# Of course, since $M$ has dimension 2, all 3-forms vanish identically:

# In[133]:


eps.exterior_derivative().display()


# In[ ]: