# Source code for coordinate_transformation

"""
Compute the coordinate transformation related variables for a mapping,

.. math::

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

We get an instance of class CoordinateTransformation. 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.py]</contents/LIBRARY/ptc/mathischeap_ptc/coordinate_transformation.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 Phi(xi, et, sg):  # define a mapping; for testing
r, s, t = (xi + 1) / 2, (et + 1) / 2, (sg + 1) / 2
x = r + 0.1 * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
y = s + 0.1 * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
z = t + 0.1 * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
return x, y, z

def d_Phi(xi, et, sg):  # the Jacobian matrix of above mapping
r, s, t = (xi + 1) / 2, (et + 1) / 2, (sg + 1) / 2
x_xi = (1 + 0.2*pi*cos(2*pi * r) * sin(2*pi * s) * sin(2*pi * t))/2
x_et = (    0.2*pi*sin(2*pi * r) * cos(2*pi * s) * sin(2*pi * t))/2
x_sg = (    0.2*pi*sin(2*pi * r) * sin(2*pi * s) * cos(2*pi * t))/2
y_xi = (    0.2*pi*cos(2*pi * r) * sin(2*pi * s) * sin(2*pi * t))/2
y_et = (1 + 0.2*pi*sin(2*pi * r) * cos(2*pi * s) * sin(2*pi * t))/2
y_sg = (    0.2*pi*sin(2*pi * r) * sin(2*pi * s) * cos(2*pi * t))/2
z_xi = (    0.2*pi*cos(2*pi * r) * sin(2*pi * s) * sin(2*pi * t))/2
z_et = (    0.2*pi*sin(2*pi * r) * cos(2*pi * s) * sin(2*pi * t))/2
z_sg = (1 + 0.2*pi*sin(2*pi * r) * sin(2*pi * s) * cos(2*pi * t))/2
return (x_xi, x_et, x_sg), (y_xi, y_et, y_sg), (z_xi, z_et, z_sg)

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

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

.. math::

\\Phi = (\\Phi_x, \\Phi_y, \\Phi_z).

:type Phi: function
:param d_Phi: The Jacobian matrix of :math:\\Phi,
:math:\\mathcal{J}. It should be
a function return the 9 components of the Jacobian matrix, i.e.,

(
:math:\\dfrac{\\partial x}{\\partial\\xi},
:math:\\dfrac{\\partial x}{\\partial\\eta},
:math:\\dfrac{\\partial x}{\\partial\\varsigma}),
(
:math:\\dfrac{\\partial y}{\\partial\\xi},
:math:\\dfrac{\\partial y}{\\partial\\eta},
:math:\\dfrac{\\partial y}{\\partial\\varsigma}),
(
:math:\\dfrac{\\partial z}{\\partial\\xi},
:math:\\dfrac{\\partial z}{\\partial\\eta},
:math:\\dfrac{\\partial z}{\\partial\\varsigma}).
:type d_Phi: function

:example:

>>> ct = CoordinateTransformation(Phi, d_Phi) # generate a CT instance
>>> xi = np.linspace(-0.9, 0.6, 3)
>>> et = np.linspace(-0.7, 0.9, 3)
>>> sg = np.linspace(-0.8, 0.7, 3)
>>> xi, et, sg = np.meshgrid(xi, et, sg, indexing='ij')
>>> x, y, z = ct.mapping(xi, et, sg)
>>> x # doctest: +ELLIPSIS
array([[[0.06469463, 0.05391086, 0.02977458],...
>>> y # doctest: +ELLIPSIS
array([[[0.16469463, 0.15391086, 0.12977458],...
>>> z # doctest: +ELLIPSIS
array([[[0.11469463, 0.47891086, 0.82977458],...
>>> J = ct.Jacobian_matrix(xi, et, sg)
>>> J[0][0] # J_{1,1} # doctest: +ELLIPSIS
array([[[0.64207986, 0.53781345, 0.30444385],...
>>> J[1][1] # J_{2,2} # doctest: +ELLIPSIS
array([[[0.53354051, 0.50892654, 0.45383545],...
>>> detJ = ct.Jacobian(xi, et, sg)
>>> detJ # doctest: +ELLIPSIS
array([[[0.1847901 , 0.11729178, 0.07611096],...
>>> iJ = ct.inverse_Jacobian_matrix(xi, et, sg)
>>> iJ[2][2] # J^{-1}_{3,3} # doctest: +ELLIPSIS
array([[[1.82807508, 2.33068327, 1.69672867],...
>>> iJ[1][2] # J^{-1}_{2,3} # doctest: +ELLIPSIS
array([[[-0.17192492,  0.33068327, -0.30327133],...
>>> G = ct.metric_matrix(xi, et, sg)
>>> G[0][1] # g_{1,2} # doctest: +ELLIPSIS
array([[[ 0.10210648,  0.02438263, -0.09377707],...
>>> G[1][2] # g_{2,3} # doctest: +ELLIPSIS
array([[[ 0.05493377, -0.03640053, -0.0063935 ],...
>>> g = ct.metric(xi, et, sg)
>>> g # doctest: +ELLIPSIS
array([[[0.03414738, 0.01375736, 0.00579288],...
>>> iG = ct.inverse_metric_matrix(xi, et, sg)
>>> iG[1][2] # g^{2,3} # doctest: +ELLIPSIS
array([[[-0.33977063,  0.72204401,  1.83434454],...

"""
def __init__(self, Phi, d_Phi):
self.Phi = Phi
self.d_Phi = d_Phi

[docs]
def mapping(self, xi, et, sg):
"""A wrap of the input mapping, Phi, i.e., :math:\\Phi.

:param xi: :math:\\xi.
:param et: :math:\\eta.
:param sg: :math:\\varsigma.
:return: a tuple :math:(x,y,z) = \\Phi(\\xi,\\eta,\\varsigma).
"""
return self.Phi(xi, et, sg)

[docs]
def Jacobian_matrix(self, xi, et, sg):
"""A wrap of the input Jacobian matrix, d_Phi, i.e.,
:math:\\mathcal{J}.

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

.. math::
\\mathcal{J} =
\\begin{bmatrix}
\\dfrac{\\partial x}{\\partial \\xi} &
\\dfrac{\\partial x}{\\partial \\eta} &
\\dfrac{\\partial x}{\\partial \\varsigma} \\\\
\\dfrac{\\partial y}{\\partial \\xi} &
\\dfrac{\\partial y}{\\partial \\eta} &
\\dfrac{\\partial y}{\\partial \\varsigma} \\\\
\\dfrac{\\partial z}{\\partial \\xi} &
\\dfrac{\\partial z}{\\partial \\eta} &
\\dfrac{\\partial z}{\\partial \\varsigma}
\\end{bmatrix}\\,.
"""
return self.d_Phi(xi, et, sg)

[docs]
def Jacobian(self, xi, et, sg):
"""Compute the Jacobian :math:\\mathrm{det}(\\mathcal{J}) of
the mapping.

:return: The Jacobian :math:\\mathrm{det}(\\mathcal{J}).
"""
J = self.Jacobian_matrix(xi, et, sg)
return + J[0][0]*J[1][1]*J[2][2] + J[0][1]*J[1][2]*J[2][0] \
+ J[0][2]*J[1][0]*J[2][1] - J[0][0]*J[1][2]*J[2][1] \
- J[0][1]*J[1][0]*J[2][2] - J[0][2]*J[1][1]*J[2][0]

[docs]
def metric_matrix(self, xi, et, sg):
"""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_{1,3} \\\\
g_{2,1} &
g_{2,2} &
g_{2,3} \\\\
g_{3,1} &
g_{3,2} &
g_{3,3}
\\end{bmatrix}\\,.

"""
J = self.Jacobian_matrix(xi, et, sg)
G = [[None for _ in range(3)] for __ in range(3)]
for i in range(3):
for j in range(i, 3):
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, xi, et, sg):
"""The metric
:math:g:= \\mathrm{det}(\\mathcal{G}):=(\mathrm{det}(\\mathcal{J}))^2.
"""
detJ = self.Jacobian(xi, et, sg)
return detJ ** 2

[docs]
def inverse_Jacobian_matrix(self, xi, et, sg):
"""Compute the Jacobian matrix of the inverse mapping
:math:\\Phi^{-1},

.. math::

(\\xi,\\eta,\\varsigma) = \\Phi^{-1}(x,y,z).

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

.. math::
\\mathcal{J}^{-1} =
\\begin{bmatrix}
\\dfrac{\\partial \\xi}{\\partial x} &
\\dfrac{\\partial \\xi}{\\partial y} &
\\dfrac{\\partial \\xi}{\\partial z} \\\\
\\dfrac{\\partial \\eta}{\\partial x} &
\\dfrac{\\partial \\eta}{\\partial y} &
\\dfrac{\\partial \\eta}{\\partial z} \\\\
\\dfrac{\\partial \\varsigma}{\\partial x} &
\\dfrac{\\partial \\varsigma}{\\partial y} &
\\dfrac{\\partial \\varsigma}{\\partial z}
\\end{bmatrix}\\,.

"""
J = self.Jacobian_matrix(xi, et, sg)
Jacobian = + J[0][0]*J[1][1]*J[2][2] + J[0][1]*J[1][2]*J[2][0] \
+ J[0][2]*J[1][0]*J[2][1] - J[0][0]*J[1][2]*J[2][1] \
- J[0][1]*J[1][0]*J[2][2] - J[0][2]*J[1][1]*J[2][0]
reciprocalJacobian = 1 / Jacobian
del Jacobian
iJ00 = reciprocalJacobian * (J[1][1] * J[2][2] - J[1][2] * J[2][1])
iJ01 = reciprocalJacobian * (J[2][1] * J[0][2] - J[2][2] * J[0][1])
iJ02 = reciprocalJacobian * (J[0][1] * J[1][2] - J[0][2] * J[1][1])
iJ10 = reciprocalJacobian * (J[1][2] * J[2][0] - J[1][0] * J[2][2])
iJ11 = reciprocalJacobian * (J[2][2] * J[0][0] - J[2][0] * J[0][2])
iJ12 = reciprocalJacobian * (J[0][2] * J[1][0] - J[0][0] * J[1][2])
iJ20 = reciprocalJacobian * (J[1][0] * J[2][1] - J[1][1] * J[2][0])
iJ21 = reciprocalJacobian * (J[2][0] * J[0][1] - J[2][1] * J[0][0])
iJ22 = reciprocalJacobian * (J[0][0] * J[1][1] - J[0][1] * J[1][0])
return [[iJ00, iJ01, iJ02],
[iJ10, iJ11, iJ12],
[iJ20, iJ21, iJ22]]

[docs]
def inverse_Jacobian(self, xi, et, sg):
"""Compute the inverse Jacobian, the determinant of the inverse
Jacobian matrix, :math:\\mathrm{det}(\\mathcal{J}^{-1})."""
iJ = self.inverse_Jacobian_matrix(xi, et, sg)
# noinspection PyUnresolvedReferences
return + iJ[0][0]*iJ[1][1]*iJ[2][2] + iJ[0][1]*iJ[1][2]*iJ[2][0]\
+ iJ[0][2]*iJ[1][0]*iJ[2][1] - iJ[0][0]*iJ[1][2]*iJ[2][1]\
- iJ[0][1]*iJ[1][0]*iJ[2][2] - iJ[0][2]*iJ[1][1]*iJ[2][0]

[docs]
def inverse_metric_matrix(self, xi, et, sg):
"""Compute the inverse metric matrix,
:math:\\mathcal{G}^{-1}.

:return: Return the inverse metric matrix :
:math:\\mathcal{G}^{-1}:

.. math::
\\mathcal{G}^{-1} =
\\begin{bmatrix}
g^{1,1} &
g^{1,2} &
g^{1,3} \\\\
g^{2,1} &
g^{2,2} &
g^{2,3} \\\\
g^{3,1} &
g^{3,2} &
g^{3,3}
\\end{bmatrix}\\,.
"""
iJ = self.inverse_Jacobian_matrix(xi, et, sg)
iG = [[None for _ in range(3)] for __ in range(3)]
for i in range(3):
for j in range(i, 3):
# noinspection PyTypeChecker
iG[i][j] = iJ[i][0] * iJ[j][0]
for l in range(1, 3):
iG[i][j] += iJ[i][l] * iJ[j][l]
if i != j:
iG[j][i] = iG[i][j]
return iG

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