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
required packages or download required scripts. © mathischeap.com
"""
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 = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
plt.show()
