# Source code for mimetic_basis_polynomials

"""
Node, edge, face and volume polynomials in the 3d reference
domain :math:\Omega_{\mathrm{ref}}=[-1,1]^3.

.. _grid:

====
grid
====

Let A=(a1, a2), B=(b1,b2,b3), C=(c1,c2). Then
:ref:grid (A, B, C) refers to a sequence of coordinates,
(a1, b1, c1), (a2, b1, c1), (a1, b2, c1), (a2, b2, c1),
(a1, b3, c1), (a2, b3, c1), (a1, b1, c2), ...... Namely, we
first do a meshgrid (A, B, C), then put the coordinates
into a sequence one by one picking along A direction firstly, B
direction secondly and finally C direction. Also see :func:grid.

⭕ To access the source code, click on the [source] button at the right
side or click on
:download:[mimetic_basis_polynomials.py]</contents/LIBRARY/ptc/mathischeap_ptc/mimetic_basis_polynomials.py>.
Dependence may exist. In case of error, check import and install required

"""

import numpy as np
from Lagrange_and_edge_polynomials import Lagrange_polynomials as Lp1
from Lagrange_and_edge_polynomials import edge_polynomials as ep1

[docs]class MimeticBasisPolynomials(object):
"""A wrapper of basis node, edge, face and volume polynomials in the
3d reference domain :math:\Omega_{\mathrm{ref}}=[-1,1]^3.

:param nodes_xi: The nodes on which the 1D mimetic polynomials are
built along the first axis (:math:\\xi).
:param nodes_et: The nodes on which the 1D mimetic polynomials are
built along the second axis (:math:\\eta).
:param nodes_sg: The nodes on which the 1D mimetic polynomials
are built along the third axis (:math:\\varsigma).
:type nodes_xi: 1d np.array
:type nodes_et: 1d np.array
:type nodes_sg: 1d np.array

:example:

>>> bf = MimeticBasisPolynomials('Lobatto-3', 'Lobatto-3', 'Lobatto-3')
>>> bf.degree # N = N_xi = N_eta = N_sigma = 3
[3, 3, 3]
>>> xi = np.linspace(-1, 1, 5)
>>> et = np.linspace(-1, 1, 6)
>>> sg = np.linspace(-1, 1, 7)
>>> NP = bf.node_polynomials(xi, et, sg)
>>> NP.shape # 4^3=64 node polynomials evaluated at 5*6*7=210 points
(64, 210)
>>> EP_xi, EP_et, EP_sg = bf.edge_polynomials(xi, et, sg)
>>> EP_xi.shape # 3*4*4=48 edge polynomials
(48, 210)
>>> FP_xi, FP_et, FP_sg = bf.face_polynomials(xi, et, sg)
>>> FP_et.shape # 3*4*3=36 face polynomials
(36, 210)
>>> VP = bf.volume_polynomials(xi, et, sg)
>>> VP.shape # 3^3=27 volume polynomials
(27, 210)

"""
def __init__(self, nodes_xi, nodes_et, nodes_sg):
self.nodes = [self._parse_nodes_(nodes_xi),
self._parse_nodes_(nodes_et),
self._parse_nodes_(nodes_sg)]
self.degree = [len(self.nodes[i])-1 for i in range(3)]

@staticmethod
def _parse_nodes_(nodes):
"""Parse the nodes. You can customize your own input shortcut
through this method.
"""
if isinstance(nodes, str):
if nodes[:7] == 'Lobatto': # for example, nodes = "Lobatto-5"
_, nodes_degree = nodes.split('-')
nodes_degree = int(nodes_degree)
nodes, _ = Lobatto(nodes_degree)
else:
pass
return nodes

[docs]    def node_polynomials(self, xi, et, sg):
"""Evaluate the node polynomials at grid(xi, et, sg).

:param xi: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param et: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param sg: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:return: A 2d np.array whose 0-dimension refers to the number
of node polynomials and 1-dimension refers to the values
evaluated at :ref:grid (xi, eta, sigma).
"""
nodes_xi, nodes_et, nodes_sg = self.nodes
np_xi = Lp1(nodes_xi, xi)
np_eta = Lp1(nodes_et, et)
np_sigma = Lp1(nodes_sg, sg)
return np.kron(np_sigma, np.kron(np_eta, np_xi))

[docs]    def edge_polynomials(self, xi, et, sg):
"""Evaluate the edge polynomials at grid(xi, et, sg).

:param xi: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param et: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param sg: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:return: A 2d np.array whose 0-dimension refers to the number
of node polynomials and 1-dimension refers to the values
evaluated at :ref:grid (xi, eta, sigma).
"""
nodes_xi, nodes_et, nodes_sg = self.nodes
np_xi = Lp1(nodes_xi, xi)
np_eta = Lp1(nodes_et, et)
np_sigma = Lp1(nodes_sg, sg)
ep_xi = ep1(nodes_xi, xi)
ep_eta = ep1(nodes_et, et)
ep_sigma = ep1(nodes_sg, sg)
return np.kron(np_sigma, np.kron(np_eta, ep_xi)), \
np.kron(np_sigma, np.kron(ep_eta, np_xi)), \
np.kron(ep_sigma, np.kron(np_eta, np_xi)),

[docs]    def face_polynomials(self, xi, et, sg):
"""Evaluate the face polynomials at grid(xi, et, sg).

:param xi: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param et: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param sg: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:return: A 2d np.array whose 0-dimension refers to the number
of node polynomials and 1-dimension refers to the values
evaluated at :ref:grid (xi, eta, sigma).
"""
nodes_xi, nodes_et, nodes_sg = self.nodes
np_xi = Lp1(nodes_xi, xi)
np_eta = Lp1(nodes_et, et)
np_sigma = Lp1(nodes_sg, sg)
ep_xi = ep1(nodes_xi, xi)
ep_eta = ep1(nodes_et, et)
ep_sigma = ep1(nodes_sg, sg)
return np.kron(ep_sigma, np.kron(ep_eta, np_xi)), \
np.kron(ep_sigma, np.kron(np_eta, ep_xi)), \
np.kron(np_sigma, np.kron(ep_eta, ep_xi)),

[docs]    def volume_polynomials(self, xi, et, sg):
"""Evaluate the volume polynomials at grid(xi, et, sg).

:param xi: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param et: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:param sg: :ref:grid (xi, eta, sigma) is the grid
to evaluate the polynomials.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:return: A 2d np.array whose 0-dimension refers to the number
of node polynomials and 1-dimension refers to the values
evaluated at :ref:grid (xi, eta, sigma).
"""
nodes_xi, nodes_et, nodes_sg = self.nodes
ep_xi = ep1(nodes_xi, xi)
ep_eta = ep1(nodes_et, et)
ep_sigma = ep1(nodes_sg, sg)
return np.kron(ep_sigma, np.kron(ep_eta, ep_xi))

[docs]def grid(A, B, C):
"""The function to compute the grid of three sets of nodes.

:param A: The first (along :math:\\xi) set of nodes.
:param B: The second (along :math:\\eta) set of nodes.
:param C: The third (along :math:\\varsigma) set of nodes.
:type A: 1d data object
:type B: 1d data object
:type C: 1d data object

:return: A tuple of three outputs:

#. (1d np.array) The grided :math:\\xi coordinates.
#. (1d np.array) The grided :math:\\eta coordinates.
#. (1d np.array) The grided :math:\\varsigma coordinates.

:example:

>>> A = np.array([1, 2])
>>> B = np.array([3, 4, 5])
>>> C = np.array([6, 7])
>>> x, y, z = grid(A, B, C)
>>> D = np.vstack((x,y,z)).T
>>> D
array([[1, 3, 6],
[2, 3, 6],
[1, 4, 6],
[2, 4, 6],
[1, 5, 6],
[2, 5, 6],
[1, 3, 7],
[2, 3, 7],
[1, 4, 7],
[2, 4, 7],
[1, 5, 7],
[2, 5, 7]])
"""
x, y, z = np.meshgrid(A, B, C, indexing='ij')
x = x.ravel('F')
y = y.ravel('F')
z = z.ravel('F')
return x, y, z

if __name__ == '__main__':
import doctest
doctest.testmod()