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

# # Laplace-Beltrami operator
# 
# This worksheet explores the issue raised by this [ask.sagemath question](https://ask.sagemath.org/question/37513/better-implementation-for-laplace-beltrami-operator/)

# In[1]:


version()


# In[2]:


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


# In[3]:


Parallelism().set(nproc=1)


# ## Manifold and coordinate charts

# In[4]:


M = Manifold(2*3,'R^6',field='real',start_index=1)


# In[5]:


m1, m2 = var('m1 m2', domain='real')


# In[6]:


m_CM = m1+m2; mu1 = m1/m_CM; mu2 = m2/m_CM; mu = m1*m2/m_CM


# In[7]:


c_Cart.<x1,y1,z1,x2,y2,z2> = M.chart()
c_Cart


# In[8]:


g = M.riemannian_metric('g')
g[1,1],g[2,2],g[3,3], g[4,4],g[5,5],g[6,6] = m1,m1,m1, m2,m2,m2; 
g.display()


# In[9]:


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


# In[10]:


ch_Cart_CM = c_Cart.transition_map(c_CM, [mu1*x1+mu2*x2,mu1*y1+mu2*y2,mu1*z1+mu2*z2, x1-x2,y1-y2,z1-z2],
                                   restrictions2 = [x!=0, y!=0, z!=0])
ch_Cart_CM.display()


# In[11]:


ch_Cart_CM.inverse().display()


# In[12]:


M.set_default_chart(c_CM)
M.set_default_frame(c_CM.frame())


# ## The scalar field $\psi$

# In[13]:


psiX = M.scalar_field({c_CM: function('Psi_X')(X,Y,Z)}, name='psiX',
                      latex_name='\Psi_X')
psiX.display()


# In[14]:


psix = M.scalar_field({c_CM: function('psi_x')(x,y,z)}, name='psix',
                      latex_name='\psi_x')
psix.display()


# In[15]:


psi = psiX*psix
psi.set_name('psi', latex_name='\psi')
psi.display()


# ## Laplace-Beltrami operator as $*\mathrm{d}*\mathrm{d}\psi$

# In[16]:


get_ipython().run_line_magic('time', 'dpsi = psi.exterior_derivative()')
print(dpsi)


# In[17]:


dpsi.display()


# In[18]:


get_ipython().run_line_magic('time', 'sdpsi = dpsi.hodge_dual(g)')
print(sdpsi)


# In[19]:


get_ipython().run_line_magic('time', 'dsdpsi = sdpsi.exterior_derivative()')
print(dsdpsi)


# In[20]:


dsdpsi.display()


# In[21]:


get_ipython().run_line_magic('time', 'LBpsi = dsdpsi.hodge_dual(g)')
print(LBpsi)


# ## Laplace-Beltrami operator as $\nabla_\mu \nabla^\mu \psi$

# In[22]:


nabla = g.connection()
print(nabla)


# In[23]:


get_ipython().run_line_magic('time', 'LBpsi1 = nabla(nabla(psi).up(g)).trace()')
LBpsi1


# In[24]:


LBpsi1.coord_function()


# Test:

# In[25]:


(LBpsi1 - LBpsi).display()


# In[ ]: