"""
In this script, we define a mesh in a unit cube,
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[crazy_mesh_hybrid.py]</contents/LIBRARY/ptc/mathischeap_ptc/crazy_mesh_hybrid.py>`.
Dependence may exist. In case of error, check import and install required
packages or download required scripts. © mathischeap.com
"""
from abc import ABC
import numpy as np
from crazy_mesh import CrazyMesh, CrazyMeshGlobalNumbering, CrazyMeshGlobalBoundaryDOFs
from coordinate_transformation_surface import extract_surface_coordinate_transformations_of
[docs]
class CrazyMeshHybrid(CrazyMesh):
"""The hybrid version of class :class:`CrazyMesh`.
"""
[docs]
def surface_CT_of_element_index(self, surface_name, i, j, k):
"""For the ``surface_name`` side of element
:math:`\\Omega_{i,j,k}`.
:param surface_name: 'N' (:math:`\\xi^-`), 'S' (:math:`\\xi^+`),
'W'(:math:`\\eta^-`), 'E' (:math:`\\eta^+`),
'B' (:math:`\\varsigma^-`) or 'F' (:math:`\\varsigma^+`)?
:param i: Element index ``i``.
:param j: Element index ``j``.
:param k: Element index ``k``.
:return:
"""
assert len(surface_name) == 1 and surface_name in 'NSWEBF', \
f"surface_name={surface_name} wrong, be one of NSWEBF."
indicator = surface_name + ':' + str(i) + '-' + str(j) + '-' + str(k)
if indicator in self._cache_:
return self._cache_[indicator]
element_CT = self.CT_of_element_index(i, j, k)
surface_CT = extract_surface_coordinate_transformations_of(
element_CT, which_sides=surface_name)[0]
self._cache_[indicator] = surface_CT
return surface_CT
[docs]
def surface_CT_of_element_number(self, surface_name, m):
"""For the ``surface_name`` side of element
:math:`\\Omega_{m}`.
:param surface_name: 'N' (:math:`\\xi^-`), 'S' (:math:`\\xi^+`),
'W'(:math:`\\eta^-`), 'E' (:math:`\\eta^+`),
'B' (:math:`\\varsigma^-`) or 'F' (:math:`\\varsigma^+`)?
:param m: Element :math:`\\Omega_{m}`.
:return:
"""
K = self.K
assert 0 <= m < K**3 and m % 1 == 0, \
f"m={m} is wrong, must be an integer between 0 and {K**3}."
k = m // K**2
j = ( m - k * K**2 ) // K
i = m - k * K**2 - j * K
return self.surface_CT_of_element_index(surface_name, i, j, k)
[docs]
class CrazyMeshHybridGlobalNumbering(CrazyMeshGlobalNumbering):
"""The hybrid version of class :class:`CrazyMeshGlobalNumbering`.
"""
@property
def NP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{NP}_{N}(\\Omega)` on a crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
num_basis = (N + 1) ** 3
total_dofs = num_basis * K**3
GM_NP = np.arange(total_dofs, dtype='int').reshape(
(K ** 3, num_basis), order='C'
)
return GM_NP
@property
def EP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{EP}_{N-1}(\\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
num_basis = N * (N + 1)**2 * 3
total_dofs = num_basis * K**3
GM_EP = np.arange(total_dofs, dtype='int').reshape(
(K ** 3, num_basis), order='C'
)
return GM_EP
@property
def FP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{FP}_{N-1}(\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
num_basis = N**2 * (N + 1) * 3
total_dofs = num_basis * K**3
GM_FP = np.arange(total_dofs, dtype='int').reshape(
(K ** 3, num_basis), order='C'
)
return GM_FP
@property
def VP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{VP}_{N-1}(\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
Note that for :math:`\\text{VP}_{N-1}(\Omega)`, we obviously
can use the numbering for the non-hybrid crazy mesh such it is
discontinuous any way.
"""
K, N = self.K, self.N
num_basis = N**3
total_dofs = num_basis * K**3
GM_FP = np.arange(total_dofs, dtype='int').reshape(
(K ** 3, num_basis), order='C'
)
return GM_FP
@property
def TN(self):
"""Generate a global numbering for a discrete trace variable in
:math:`\\text{TN}_{N}(\\partial\\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
"""
raise NotImplementedError("Could you code it?")
@property
def TE(self):
"""Generate a global numbering for a discrete trace variable in
:math:`\\text{TE}_{N-1}(\\partial\\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
"""
raise NotImplementedError("Could you code it?")
@property
def TF(self):
"""Generate a global numbering for a discrete trace variable in
:math:`\\text{TF}_{N-1}(\\partial\\Omega)` on a hybrid crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
NN_KK = N**2 * K**2
NK = N * K
x_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((K+1, NK, NK), order='F')
y_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((NK, K+1, NK), order='F') + NN_KK * (K+1)
z_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((NK, NK, K+1), order='F') + 2 * NN_KK * (K+1)
GM_TF = np.zeros((K ** 3, 6 * N**2), dtype='int') # K^3 elements, on each element has 6 sides, and each sides have N^2 faces
for k in range(K):
for j in range(K):
for i in range(K):
m = i + j * K + k * K**2
GM_TF[m] = np.concatenate(
(x_face[i , j * N:(j + 1) * N, k * N:(k + 1) * N].ravel('F'),
x_face[i+1, j * N:(j + 1) * N, k * N:(k + 1) * N].ravel('F'),
y_face[i * N:(i + 1) * N, j , k * N:(k + 1) * N].ravel('F'),
y_face[i * N:(i + 1) * N, j+1, k * N:(k + 1) * N].ravel('F'),
z_face[i * N:(i + 1) * N, j * N:(j + 1) * N, k ].ravel('F'),
z_face[i * N:(i + 1) * N, j * N:(j + 1) * N, k+1].ravel('F')))
return GM_TF
[docs]
class CrazyMeshHybridGlobalBoundaryDOFs(CrazyMeshGlobalBoundaryDOFs, ABC):
"""The hybrid version of class :class:`CrazyMeshGlobalBoundaryDOFs`.
"""
@property
def TN(self):
"""Find the dofs of an element in
:math:`\\text{TN}_{N}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def TE(self):
"""Find the dofs of an element in
:math:`\\text{TE}_{N-1}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def TF(self):
"""Find the dofs of an element in
:math:`\\text{TF}_{N-1}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
K, N = self.K, self.N
NN_KK = N**2 * K**2
NK = N * K
x_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((K+1, NK, NK), order='F')
y_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((NK, K+1, NK), order='F') + NN_KK * (K+1)
z_face = np.arange(NN_KK * (K+1), dtype='int'
).reshape((NK, NK, K+1), order='F') + 2 * NN_KK * (K+1)
DOFs_on_boundary = dict()
DOFs_on_boundary['x_minus'] = x_face[0, :, :].ravel('F')
DOFs_on_boundary['x_plus'] = x_face[-1, :, :].ravel('F')
DOFs_on_boundary['y_minus'] = y_face[:, 0, :].ravel('F')
DOFs_on_boundary['y_plus'] = y_face[:, -1, :].ravel('F')
DOFs_on_boundary['z_minus'] = z_face[:, :, 0].ravel('F')
DOFs_on_boundary['z_plus'] = z_face[:, :, -1].ravel('F')
return DOFs_on_boundary
[docs]
class CrazyMeshHybridLocalBoundaryDOFs(object):
"""Like the :class:`CrazyMeshHybridGlobalBoundaryDOFs`, but we now
return the local numbering on each boundary.
:param K: The crazy mesh is of :math:`K^3` elements.
:param N: The degree :math:`N`. of the to be used mimetic polynomial
basis functions.
:type K: int
:type N: int
"""
def __init__(self, K, N):
self.K, self.N = K, N
@property
def NP(self):
"""Find the dofs of an element in
:math:`\\text{NP}_{N}(\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def EP(self):
"""Find the dofs of an element in
:math:`\\text{EP}_{N-1}(\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def FP(self):
"""Find the dofs of an element in
:math:`\\text{FP}_{N-1}(\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
K, N = self.K, self.N
local_DOFs_on_boundary = dict()
local_DOFs_on_boundary['x_minus'] = dict()
local_DOFs_on_boundary['x_plus'] = dict()
local_DOFs_on_boundary['y_minus'] = dict()
local_DOFs_on_boundary['y_plus'] = dict()
local_DOFs_on_boundary['z_minus'] = dict()
local_DOFs_on_boundary['z_plus'] = dict()
elements_global_numbering = np.arange(K ** 3, dtype='int').reshape(
(K, K, K), order='F'
)
local_numbering_x = np.arange(N ** 2 * (N+1), dtype='int').reshape(
(N+1, N, N), order='F'
)
local_numbering_y = np.arange(N ** 2 * (N+1), 2* N ** 2 * (N+1), dtype='int').reshape(
(N, N+1, N), order='F'
)
local_numbering_z = np.arange(2 *N ** 2 * (N+1), 3 *N ** 2 * (N+1), dtype='int').reshape(
(N, N, N+1), order='F'
)
# elements near North, x^-
elements = elements_global_numbering[0, :, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['x_minus'][i] = local_numbering_x[0, :, :].ravel('F')
# elements near South, x^+
elements = elements_global_numbering[-1, :, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['x_plus'][i] = local_numbering_x[-1, :, :].ravel('F')
# elements near North, y^-
elements = elements_global_numbering[:, 0, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['y_minus'][i] = local_numbering_y[:, 0, :].ravel('F')
# elements near South, y^+
elements = elements_global_numbering[:, -1, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['y_plus'][i] = local_numbering_y[:, -1, :].ravel('F')
# elements near North, z^-
elements = elements_global_numbering[:, :, 0].ravel('F')
for i in elements:
local_DOFs_on_boundary['z_minus'][i] = local_numbering_z[:, :, 0].ravel('F')
# elements near South, z^+
elements = elements_global_numbering[:, :, -1].ravel('F')
for i in elements:
local_DOFs_on_boundary['z_plus'][i] = local_numbering_z[:, :, -1].ravel('F')
return local_DOFs_on_boundary
@property
def TN(self):
"""Find the dofs of an element in
:math:`\\text{TN}_{N}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def TE(self):
"""Find the dofs of an element in
:math:`\\text{TE}_{N-1}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
raise NotImplementedError("Could you code it?")
@property
def TF(self):
"""Find the dofs of an element in
:math:`\\text{TF}_{N-1}(\\partial\\Omega)` which are on boundary of the
crazy mesh of :math:`K^3` elements.
:returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus',
'y_plus', 'z_minus', 'z_plus', and whose values are the
global numbering of the dofs on the boundaries indicated
by the keys.
"""
K, N = self.K, self.N
local_DOFs_on_boundary = dict()
local_DOFs_on_boundary['x_minus'] = dict()
local_DOFs_on_boundary['x_plus'] = dict()
local_DOFs_on_boundary['y_minus'] = dict()
local_DOFs_on_boundary['y_plus'] = dict()
local_DOFs_on_boundary['z_minus'] = dict()
local_DOFs_on_boundary['z_plus'] = dict()
elements_global_numbering = np.arange(K**3, dtype='int').reshape(
(K, K, K), order='F'
)
# elements near North, x^-
elements = elements_global_numbering[0, :, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['x_minus'][i] = [_ for _ in range(N**2)]
# elements near South, x^+
elements = elements_global_numbering[-1, :, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['x_plus'][i] = [N**2+_ for _ in range(N**2)]
# elements near North, y^-
elements = elements_global_numbering[:, 0, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['y_minus'][i] = [2*N**2+_ for _ in range(N**2)]
# elements near South, y^+
elements = elements_global_numbering[:, -1, :].ravel('F')
for i in elements:
local_DOFs_on_boundary['y_plus'][i] = [3*N**2+_ for _ in range(N**2)]
# elements near North, z^-
elements = elements_global_numbering[:, :, 0].ravel('F')
for i in elements:
local_DOFs_on_boundary['z_minus'][i] = [4*N**2+_ for _ in range(N**2)]
# elements near South, z^+
elements = elements_global_numbering[:, :, -1].ravel('F')
for i in elements:
local_DOFs_on_boundary['z_plus'][i] = [5*N**2+_ for _ in range(N**2)]
return local_DOFs_on_boundary
if __name__ == '__main__':
import doctest
doctest.testmod()
K = 2
N = 2
GN = CrazyMeshHybridGlobalNumbering(K, N)
BN = CrazyMeshHybridGlobalBoundaryDOFs(K, N)
LBN = CrazyMeshHybridLocalBoundaryDOFs(K, N)
tf = LBN.TF
print(tf)
tf = LBN.FP
print(tf)
