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()