# Source code for coordinate_transformation_surface

"""
Compute the coordinate transformation related variables for a
surface mapping, :math:\\Psi,

.. math::

\\Psi : \\Pi_{\\mathrm{ref}}\\to\\Gamma

We get an instance of class CoordinateTransformationSurface. The
transformations are its methods.

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

"""

import numpy as np
from numpy import sin as sin
from numpy import cos as cos
from numpy import pi as pi

def Psi(rho, tau):  # define a surface mapping; for testing
x = rho + 0.1 * sin(pi * rho) * sin(pi * tau)
y =       0.1 * sin(pi * rho) * sin(pi * tau)
z = tau + 0.1 * sin(pi * rho) * sin(pi * tau)
return x, y, z

def d_Psi(rho, tau): # the Jacobian matrix of above surface mapping
x_rho = 1 + 0.1 * pi * cos(pi * rho) * sin(pi * tau)
x_tau =     0.1 * pi * sin(pi * rho) * cos(pi * tau)
y_rho =     0.1 * pi * cos(pi * rho) * sin(pi * tau)
y_tau =     0.1 * pi * sin(pi * rho) * cos(pi * tau)
z_rho =     0.1 * pi * cos(pi * rho) * sin(pi * tau)
z_tau = 1 + 0.1 * pi * sin(pi * rho) * cos(pi * tau)
return (x_rho, x_tau), (y_rho, y_tau), (z_rho, z_tau)

[docs]
class CoordinateTransformationSurface(object):
""" The surface coordinate transformation class.

:param Psi: The mapping :math:\\Psi. It should be a function which
returns three components of the mapping, i.e.,

.. math::

\\Psi = (\\Psi_x, \\Psi_y, \\Psi_z).

:type Psi: function
:param d_Psi: The Jacobian matrix of :math:\\Psi,
:math:\\mathcal{J}. It should be
a function return the 6 (3*2) components of the Jacobian matrix, i.e.,

(
:math:\\dfrac{\\partial x}{\\partial\\varrho},
:math:\\dfrac{\\partial x}{\\partial\\tau}),
(
:math:\\dfrac{\\partial y}{\\partial\\varrho},
:math:\\dfrac{\\partial y}{\\partial\\tau}),
(
:math:\\dfrac{\\partial z}{\\partial\\varrho},
:math:\\dfrac{\\partial z}{\\partial\\tau}).
:type d_Psi: function
:param iJM: The inverse Jacobian matrix.

:example:

>>> ctS = CoordinateTransformationSurface(Psi, d_Psi)
>>> rho = np.linspace(-1,1,20)
>>> tau = np.linspace(-1,1,20)
>>> rho, tau = np.meshgrid(rho, tau, indexing='ij')
>>> x, y, z = ctS.mapping(rho, tau)
>>> J = ctS.Jacobian_matrix(rho, tau)
>>> np.shape(J)
(3, 2, 20, 20)
>>> G = ctS.metric_matrix(rho, tau)
>>> np.shape(G)
(2, 2, 20, 20)
>>> g = ctS.metric(rho, tau)
>>> np.shape(g)
(20, 20)

"""
def __init__(self, Psi, d_Psi, iJM=None):
self.Psi = Psi
self.d_Psi = d_Psi
self.iJM = iJM

[docs]
def mapping(self, rho, tau):
"""A wrap of the input mapping, Psi, i.e., :math:\\Psi.

:param rho: :math:\\varrho.
:param tau: :math:\\tau.
:return: a tuple :math:(x,y,z) = \\Psi(\\varrho,\\tau).
"""
return self.Psi(rho, tau)

[docs]
def Jacobian_matrix(self, rho, tau):
"""A wrap of the input Jacobian matrix, d_Psi, i.e.,
:math:\\mathcal{J}.

:return: Return the Jacobian matrix :math:\\mathcal{J}:

.. math::
\\mathcal{J} =
\\begin{bmatrix}
\\dfrac{\\partial x}{\\partial \\varrho} &
\\dfrac{\\partial x}{\\partial \\tau} \\\\
\\dfrac{\\partial y}{\\partial \\varrho} &
\\dfrac{\\partial y}{\\partial \\tau} \\\\
\\dfrac{\\partial z}{\\partial \\varrho} &
\\dfrac{\\partial z}{\\partial \\tau}
\\end{bmatrix}\\,.
"""
return self.d_Psi(rho, tau)

[docs]
def inverse_Jacobian_matrix(self, rho, tau):
"""Compute the inverse Jacobian matrix,
:math:\\mathcal{J}^{-1}.

:return: Return the inverse Jacobian matrix
:math:\\mathcal{J}^-1:

.. math::
\\mathcal{J}^{-1} =
\\begin{bmatrix}
\\dfrac{\\partial \\varrho}{\\partial x} &
\\dfrac{\\partial \\varrho}{\\partial y} &
\\dfrac{\\partial \\varrho}{\\partial z} \\\\
\\dfrac{\\partial \\tau}{\\partial x} &
\\dfrac{\\partial \\tau}{\\partial y} &
\\dfrac{\\partial \\tau}{\\partial z}
\\end{bmatrix}\\,.
"""
return self.iJM(rho, tau)

[docs]
def metric_matrix(self, rho, tau):
"""Compute the metric matrix :math:\\mathcal{G}.

:return: Return the metric matrix :math:\\mathcal{G}:

.. math::
\\mathcal{G} =
\\begin{bmatrix}
g_{1,1} &
g_{1,2} \\\\
g_{2,1} &
g_{2,2}
\\end{bmatrix}\\,.

"""
J = self.Jacobian_matrix(rho, tau)
G = [[None for _ in range(2)] for __ in range(2)]
for i in range(2):
for j in range(i, 2):
G[i][j] = J[0][i] * J[0][j]
for l in range(1, 3):
G[i][j] += J[l][i] * J[l][j]
if i != j:
G[j][i] = G[i][j]
return G

[docs]
def metric(self, rho, tau):
"""The metric :math:g.
"""
G = self.metric_matrix(rho, tau)
return G[0][0] * G[1][1] - G[0][1]*G[1][0]

[docs]
def extract_surface_coordinate_transformations_of(ct, which_sides=None):
"""We extract the six boundary coordinate transformations from a
:class:CoordinateTransformation instance representing a 3D mapping
:math:\\Phi.

.. math::
\\Phi: \\Omega_{\\mathrm{ref}}\\to\\Omega

:param ct: A CoordinateTransformation instance that represents the
mapping :math:\\Phi:\\Omega_{\mathrm{ref}}\\to\\Omega.
:type ct: CoordinateTransformation
:param which_sides: (default: None) We want to extract
sub-mappings on which sides?
:return: A tuple of 6 :class:CoordinateTransformationSurface
instances representing the north (:math:\\xi^-),
south (:math:\\xi^+), west (:math:\\varsigma^-),
east (:math:\\eta^+), back (:math:\\varsigma^-) and
front (:math:\\varsigma^+) boundaries.
"""
assert ct.__class__.__name__ == 'CoordinateTransformation'

if which_sides is None:
which_sides = 'NSWEBF'

sub_mappings = list()

if 'N' in which_sides:# North side
def mapping_N(rho, tau):
return ct.mapping(-1, rho, tau)
def JM_N(rho, tau):
J = ct.Jacobian_matrix(-1, rho, tau)
return ((J[0][1], J[0][2]),
(J[1][1], J[1][2]),
(J[2][1], J[2][2]))
def iJM_N(rho, tau):
J = ct.inverse_Jacobian_matrix(-1, rho, tau)
return ((J[1][0], J[1][1], J[1][2]),
(J[2][0], J[2][1], J[2][2]))
cts_N = CoordinateTransformationSurface(mapping_N, JM_N, iJM=iJM_N)
sub_mappings.append(cts_N)

if 'S' in which_sides:# South side
def mapping_S(rho, tau):
return ct.mapping(1, rho, tau)
def JM_S(rho, tau):
J = ct.Jacobian_matrix(1, rho, tau)
return ((J[0][1], J[0][2]),
(J[1][1], J[1][2]),
(J[2][1], J[2][2]))
def iJM_S(rho, tau):
J = ct.inverse_Jacobian_matrix(1, rho, tau)
return ((J[1][0], J[1][1], J[1][2]),
(J[2][0], J[2][1], J[2][2]))
cts_S = CoordinateTransformationSurface(mapping_S, JM_S, iJM=iJM_S)
sub_mappings.append(cts_S)

if 'W' in which_sides:# West side
def mapping_W(rho, tau):
return ct.mapping(rho, -1, tau)
def JM_W(rho, tau):
J = ct.Jacobian_matrix(rho, -1, tau)
return ((J[0][2], J[0][0]),
(J[1][2], J[1][0]),
(J[2][2], J[2][0]))
def iJM_W(rho, tau):
J = ct.inverse_Jacobian_matrix(rho, -1, tau)
return ((J[2][0], J[2][1], J[2][2]),
(J[0][0], J[0][1], J[0][2]))
cts_W = CoordinateTransformationSurface(mapping_W, JM_W, iJM=iJM_W)
sub_mappings.append(cts_W)

if 'E' in which_sides:# Ease side
def mapping_E(rho, tau):
return ct.mapping(rho, 1, tau)
def JM_E(rho, tau):
J = ct.Jacobian_matrix(rho, 1, tau)
return ((J[0][2], J[0][0]),
(J[1][2], J[1][0]),
(J[2][2], J[2][0]))
def iJM_E(rho, tau):
J = ct.inverse_Jacobian_matrix(rho, 1, tau)
return ((J[2][0], J[2][1], J[2][2]),
(J[0][0], J[0][1], J[0][2]))
cts_E = CoordinateTransformationSurface(mapping_E, JM_E, iJM=iJM_E)
sub_mappings.append(cts_E)

if 'B' in which_sides:# Back side
def mapping_B(rho, tau):
return ct.mapping(rho, tau, -1)
def JM_B(rho, tau):
J = ct.Jacobian_matrix(rho, tau, -1)
return ((J[0][0], J[0][1]),
(J[1][0], J[1][1]),
(J[2][0], J[2][1]))
def iJM_B(rho, tau):
J = ct.inverse_Jacobian_matrix(rho, tau, -1)
return ((J[0][0], J[0][1], J[0][2]),
(J[1][0], J[1][1], J[1][2]))
cts_B = CoordinateTransformationSurface(mapping_B, JM_B, iJM=iJM_B)
sub_mappings.append(cts_B)

if 'F' in which_sides:# Front side
def mapping_F(rho, tau):
return ct.mapping(rho, tau, 1)
def JM_F(rho, tau):
J = ct.Jacobian_matrix(rho, tau, 1)
return ((J[0][0], J[0][1]),
(J[1][0], J[1][1]),
(J[2][0], J[2][1]))
def iJM_F(rho, tau):
J = ct.inverse_Jacobian_matrix(rho, tau, 1)
return ((J[0][0], J[0][1], J[0][2]),
(J[1][0], J[1][1], J[1][2]))
cts_F = CoordinateTransformationSurface(mapping_F, JM_F, iJM=iJM_F)
sub_mappings.append(cts_F)

return sub_mappings

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

from coordinate_transformation import CoordinateTransformation, Phi, d_Phi

ct = CoordinateTransformation(Phi, d_Phi)

N, S, W, E, B, F = extract_surface_coordinate_transformations_of(ct)
ctS = S
# ctS = CoordinateTransformationSurface(Psi, d_Psi)

rho = np.linspace(-1,1,20)
tau = np.linspace(-1,1,20)

rho, tau = np.meshgrid(rho, tau, indexing='ij')

x, y, z = ctS.mapping(rho, tau)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(8, 6))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
`