"""
Advanced Examples
=================

In this section, we demonstrate how to go from the inner product to the
discrete approximation for some special cases. We also show how all
necessary operators are constructed for each case.
    


"""

####################################################
#
# Import Packages
# ---------------
#
# Here we import the packages required for this tutorial
#

from discretize.utils import sdiag
from discretize import TensorMesh
import numpy as np
# Constitive Relations and Differential Operators
# -----------------------------------------------
#
# Where :math:`\psi` and :math:`\phi` are scalar quantities,
# :math:`\vec{u}` and :math:`\vec{v}` are vector quantities, and
# :math:`\sigma` defines a constitutive relationship, we may need to derive
# discrete approximations for the following inner products:
#
#     1. :math:`(\vec{u} , \sigma \nabla \phi)`
#     2. :math:`(\psi , \sigma \nabla \cdot \vec{v})`
#     3. :math:`(\vec{u} , \sigma \nabla \times \vec{v})`
#
# These cases effectively combine what was learned in the previous two
# tutorials. For each case, we must:
#
#     - Define discretized quantities at the appropriate mesh locations
#     - Define an inner product matrix that depends on a single constitutive parameter (:math:`\sigma`) or a tensor (:math:`\Sigma`)
#     - Construct differential operators that may require you to define boundary conditions
#
# Where :math:`\mathbf{M_e}(\sigma)` is the property dependent inner-product
# matrix for quantities on cell edges, :math:`\mathbf{M_f}(\sigma)` is the
# property dependent inner-product matrix for quantities on cell faces,
# :math:`\mathbf{G_{ne}}` is the nodes to edges gradient operator and
# :math:`\mathbf{G_{cf}}` is the centers to faces gradient operator:
#
# .. math::
#     (\vec{u} , \sigma \nabla \phi) &= \mathbf{u_f^T M_f}(\sigma) \mathbf{ G_{cf} \, \phi_c} \;\;\;\;\; (\vec{u} \;\textrm{on faces and} \; \phi \; \textrm{at centers}) \\
#     &= \mathbf{u_e^T M_e}(\sigma) \mathbf{ G_{ne} \, \phi_n} \;\;\;\; (\vec{u} \;\textrm{on edges and} \; \phi \; \textrm{on nodes})
#
# Where :math:`\mathbf{M_c}(\sigma)` is the property dependent inner-product
# matrix for quantities at cell centers and :math:`\mathbf{D}` is the faces
# to centers divergence operator:
#
# .. math::
#     (\psi , \sigma \nabla \cdot \vec{v}) = \mathbf{\psi_c^T M_c} (\sigma)\mathbf{ D v_f} \;\;\;\; (\psi \;\textrm{at centers and} \; \vec{v} \; \textrm{on faces} )
#
# Where :math:`\mathbf{C_{ef}}` is the edges to faces curl operator and
# :math:`\mathbf{C_{fe}}` is the faces to edges curl operator:
#
# .. math::
#     (\vec{u} , \sigma \nabla \times \vec{v}) &= \mathbf{u_f^T M_f} (\sigma) \mathbf{ C_{ef} \, v_e} \;\;\;\; (\vec{u} \;\textrm{on edges and} \; \vec{v} \; \textrm{on faces} )\\
#     &= \mathbf{u_e^T M_e} (\sigma) \mathbf{ C_{fe} \, v_f} \;\;\;\; (\vec{u} \;\textrm{on faces and} \; \vec{v} \; \textrm{on edges} )
#
# **With the operators constructed below, you can compute all of the
# aforementioned inner products.**
#


# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])
sig = np.random.rand(mesh.nC)  # isotropic
Sig = np.random.rand(mesh.nC, 6)  # anisotropic

# Inner product matricies
Mc = sdiag(mesh.vol * sig)  # Inner product matrix (centers)
# Mn = mesh.getNodalInnerProduct(sig)  # Inner product matrix (nodes)  (*functionality pending*)
Me = mesh.getEdgeInnerProduct(sig)  # Inner product matrix (edges)
Mf = mesh.getFaceInnerProduct(sig)  # Inner product matrix for tensor (faces)

# Differential operators
Gne = mesh.nodalGrad  # Nodes to edges gradient
mesh.setCellGradBC(["neumann", "dirichlet", "neumann"])  # Set boundary conditions
Gcf = mesh.cellGrad  # Cells to faces gradient
D = mesh.faceDiv  # Faces to centers divergence
Cef = mesh.edgeCurl  # Edges to faces curl
Cfe = mesh.edgeCurl.T  # Faces to edges curl

# EXAMPLE: (u, sig*Curl*v)
fig = plt.figure(figsize=(9, 5))

ax1 = fig.add_subplot(121)
ax1.spy(Mf * Cef, markersize=0.5)
ax1.set_title("Me(sig)*Cef (Isotropic)", pad=10)

Mf_tensor = mesh.getFaceInnerProduct(Sig)  # inner product matrix for tensor
ax2 = fig.add_subplot(122)
ax2.spy(Mf_tensor * Cef, markersize=0.5)
ax2.set_title("Me(sig)*Cef (Anisotropic)", pad=10)

#####################################################
# Divergence of a Scalar and a Vector Field
# -----------------------------------------
#
# Where :math:`\psi` and :math:`\phi` are scalar quantities, and
# :math:`\vec{u}` is a known vector field, we may need to derive
# a discrete approximation for the following inner product:
#
# .. math::
#     (\psi , \nabla \cdot \phi \vec{u})
#
# Scalar and vector quantities are generally discretized to lie on
# different locations on the mesh. As result, it is better to use the
# identity :math:`\nabla \cdot \phi \vec{u} = \phi \nabla \cdot \vec{u} + \vec{u} \cdot \nabla \phi`
# and separate the inner product into two parts:
#
# .. math::
#    (\psi , \phi \nabla \cdot \vec{u} ) + (\psi , \vec{u} \cdot \nabla \phi)
#
# **Term 1:**
#
# If the vector field :math:`\vec{u}` is divergence free, there is no need
# to evaluate the first inner product term. This is the case for advection when
# the fluid is incompressible.
#
# Where :math:`\mathbf{D_{fc}}` is the faces to centers divergence operator, and
# :math:`\mathbf{M_c}` is the basic inner product matrix for cell centered
# quantities, we can approximate this inner product as:
#
# .. math::
#     (\psi , \phi \nabla \cdot \vec{u} ) = \mathbf{\psi_c^T M_c} \textrm{diag} (\mathbf{D_{fc} u_f} ) \, \mathbf{\phi_c}
#
# **Term 2:**
#
# Let :math:`\mathbf{G_{cf}}` be the cell centers to faces gradient operator,
# :math:`\mathbf{M_c}` be the basic inner product matrix for cell centered
# quantities, and :math:`\mathbf{\tilde{A}_{fc}}` and averages *and* sums the
# cartesian contributions of :math:`\vec{u} \cdot \nabla \phi`, we can
# approximate the inner product as:
#
# .. math::
#     (\psi , \vec{u} \cdot \nabla \phi) = \mathbf{\psi_c^T M_c \tilde A_{fc}} \text{diag} (\mathbf{u_f} ) \mathbf{G_{cf} \, \phi_c}
#
# **With the operators constructed below, you can compute all of the
# inner products.**

# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])

# Inner product matricies
Mc = sdiag(mesh.vol * sig)  # Inner product matrix (centers)

# Differential operators
mesh.setCellGradBC(["neumann", "dirichlet", "neumann"])  # Set boundary conditions
Gcf = mesh.cellGrad  # Cells to faces gradient
Dfc = mesh.faceDiv  # Faces to centers divergence

# Averaging and summing matrix
Afc = mesh.dim * mesh.aveF2CC
import matplotlib.pyplot as plt


#####################################################