"""
`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
packages or download required scripts. © mathischeap.com
"""
import numpy as np
from quadrature import Lobatto
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()
