"""
In this script, we define a mesh in a unit cube,
.. math::
\\Omega = [0,1]^3.
The mesh is of :math:`K^3` elements,
.. math::
\\Omega_{m}=\\Omega_{i+(j-1)K+(k-1)K^2}=\\Omega_{i,j,k},\\quad
i,j,k \\in\\left\\lbrace1,2,\\cdots,K\\right\\rbrace.
The mapping :math:`\\Phi_{m}=\\Phi_{i,j,k}:\\Omega_{\\mathrm{ref}}\\to
\\Omega_{i,j,k}` is given as
.. math::
\\Phi_{i,j,k} = \\mathring{\\Phi}\\circ\\Xi_{i,j,k},
where :math:`\\Xi_{i,j,k}` is a linear mapping,
:math:`\\Xi_{i,j,k}:\\Omega_{\\mathrm{ref}}\\to\\left(
\\left[\\dfrac{i-1}{K},\\dfrac{i}{K}\\right],
\\left[\\dfrac{j-1}{K},\\dfrac{j}{K}\\right],
\\left[\\dfrac{k-1}{K},\\dfrac{k}{K}\\right]
\\right)`, i.e.,
.. math::
\\begin{pmatrix}
r\\\\s\\\\t
\\end{pmatrix} = \\Xi_{i,j,k}(\\xi,\\eta,\\varsigma)
=
\\dfrac{1}{K}
\\begin{pmatrix}
i-1 + (\\xi+1)/2\\\\
j-1 + (\\eta+1)/2\\\\
k-1 + (\\varsigma+1)/2
\\end{pmatrix},
and :math:`\\mathring{\\Phi}` is a mapping,
.. math::
\\begin{pmatrix}
x\\\\y\\\\z
\\end{pmatrix} = \\mathring{\\Phi}(r,s,t) =
\\begin{pmatrix}
r + \\frac{1}{2}c \\sin(2\\pi r)\\sin(2\\pi s)\\sin(2\\pi t)\\\\
s + \\frac{1}{2}c \\sin(2\\pi r)\\sin(2\\pi s)\\sin(2\\pi t)\\\\
t + \\frac{1}{2}c \\sin(2\\pi r)\\sin(2\\pi s)\\sin(2\\pi t)
\\end{pmatrix},
where :math:`0 \leq c \leq 0.25` is a deformation factor. When
:math:`c=0`, :math:`\mathring{\Phi}` is also a linear mapping. Thus we
have a uniform orthogonal mesh. When :math:`c>0`, the mesh is
curvilinear. Two examples (left: :math:`c=0`, right: :math:`c=0.25`) of
this mesh for :math:`K=3` are shown below.
.. image:: workshop/crazy_mesh.png
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[crazy_mesh.py]</contents/LIBRARY/ptc/mathischeap_ptc/crazy_mesh.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, cos, pi
from coordinate_transformation import CoordinateTransformation
[docs]
class CrazyMesh(object):
"""The crazy mesh.
:param c: The deformation factor, :math:`0 \leq c \leq 0.25`.
:type c: float
:param K: The crazy mesh is of :math:`K^3` elements.
:type K: int
:example:
>>> cm = CrazyMesh(0.1, 2)
>>> e0 = cm.CT_of_element_number(0)
>>> e000 = cm.CT_of_element_index(0, 0, 0)
>>> e0 is e000
True
>>> e7 = cm.CT_of_element_number(7)
>>> e111 = cm.CT_of_element_index(1, 1, 1)
>>> e7 is e111
True
>>> e000 # doctest: +ELLIPSIS
<coordinate_transformation.CoordinateTransformation object at...
>>> e7 # doctest: +ELLIPSIS
<coordinate_transformation.CoordinateTransformation object at...
"""
def __init__(self, c, K):
assert 0 <= c <= 0.25, \
"The deformation factor must be in [0, 0.25]."
self.c = c
self.K = K
self._cache_ = dict() # we cache the output for repeated usage.
[docs]
def CT_of_element_index(self, i, j, k):
"""Return a :class:`CoordinateTransformation` instance for
element :math:`\\Omega_{i,j,k}`.
Note that Python index starts from :math:`0`. So for a
:class:`CrazyMesh` of :math:`K^3` elements, its indices,
:math:`i,j,k\\in\\left\\lbrace 0, 1, K-1\\right\\rbrace`.
:param i: Element index ``i``.
:param j: Element index ``j``.
:param k: Element index ``k``.
:type i: int
:type j: int
:type k: int
:return: A :class:`CoordinateTransformation` instance.
"""
c, K = self.c, self.K
assert 0 <= i < K and i % 1 == 0, \
f"i={i} is wrong, must be an integer between 0 and {K}."
assert 0 <= j < K and j % 1 == 0, \
f"j={j} is wrong, must be an integer between 0 and {K}."
assert 0 <= k < K and k % 1 == 0, \
f"k={k} is wrong, must be an integer between 0 and {K}."
indicator = str(i)+'-'+str(j)+'-'+str(k)
if indicator in self._cache_: # to see if this element is cached
# if yes, then return it from the cache.
return self._cache_[indicator]
# otherwise, we make it.
def Phi(xi, et, sg):
r = ( 1 / K ) * ( i + (xi + 1) / 2 )
s = ( 1 / K ) * ( j + (et + 1) / 2 )
t = ( 1 / K ) * ( k + (sg + 1) / 2 )
x = r + 0.5* c * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
y = s + 0.5* c * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
z = t + 0.5* c * sin(2 * pi * r) * sin(2 * pi * s) * sin(2 * pi * t)
return x, y, z
def d_Phi(xi, et, sg):
r = ( 1 / K ) * ( i + (xi + 1) / 2 )
s = ( 1 / K ) * ( j + (et + 1) / 2 )
t = ( 1 / K ) * ( k + (sg + 1) / 2 )
x_xi = (1 + c * pi * cos(2 * pi * r) * sin(
2 * pi * s) * sin(2 * pi * t)) / (2 * K)
x_et = (c * pi * sin(2 * pi * r) * cos(2 * pi * s) * sin(
2 * pi * t)) / (2 * K)
x_sg = (c * pi * sin(2 * pi * r) * sin(2 * pi * s) * cos(
2 * pi * t)) / (2 * K)
y_xi = (c * pi * cos(2 * pi * r) * sin(2 * pi * s) * sin(
2 * pi * t)) / (2 * K)
y_et = (1 + c * pi * sin(2 * pi * r) * cos(
2 * pi * s) * sin(2 * pi * t)) / (2 * K)
y_sg = (c * pi * sin(2 * pi * r) * sin(2 * pi * s) * cos(
2 * pi * t)) / (2 * K)
z_xi = (c * pi * cos(2 * pi * r) * sin(2 * pi * s) * sin(
2 * pi * t)) / (2 * K)
z_et = (c * pi * sin(2 * pi * r) * cos(2 * pi * s) * sin(
2 * pi * t)) / (2 * K)
z_sg = (1 + c * pi * sin(2 * pi * r) * sin(
2 * pi * s) * cos(2 * pi * t)) / (2 * K)
return (x_xi, x_et, x_sg), (y_xi, y_et, y_sg), (
z_xi, z_et, z_sg)
ct = CoordinateTransformation(Phi, d_Phi)
self._cache_[indicator] = ct # before return, we cache it.
return ct
[docs]
def CT_of_element_number(self, m):
"""Return a :class:`CoordinateTransformation` instance for
element :math:`\\Omega_{m}`.
Note that Python index starts from :math:`0`. So for a
:class:`CrazyMesh` of :math:`K^3` elements,
:math:`m\\in\\left\\lbrace 0,1,\\cdots,K^3-1\\right\\rbrace`.
:param m: Element No. ``m``.
:type m: int
:return: A :class:`CoordinateTransformation` instance.
"""
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.CT_of_element_index(i, j, k)
[docs]
class CrazyMeshGlobalNumbering(object):
"""A wrapper of global numberings for discrete variables in the
crazy mesh.
: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
:example:
>>> K = 2
>>> N = 1
>>> GM = CrazyMeshGlobalNumbering(K, N)
>>> GM.FP
array([[ 0, 1, 12, 14, 24, 28],
[ 1, 2, 13, 15, 25, 29],
[ 3, 4, 14, 16, 26, 30],
[ 4, 5, 15, 17, 27, 31],
[ 6, 7, 18, 20, 28, 32],
[ 7, 8, 19, 21, 29, 33],
[ 9, 10, 20, 22, 30, 34],
[10, 11, 21, 23, 31, 35]])
>>> GM.VP
array([[0],
[1],
[2],
[3],
[4],
[5],
[6],
[7]])
"""
def __init__(self, K, N):
self.K, self.N = K, N
@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.
"""
raise NotImplementedError("Could you complete this property?")
@property
def EP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{EP}_{N-1}(\\Omega)` on a crazy mesh of
:math:`K^3` elements.
"""
raise NotImplementedError("Could you complete this property?")
@property
def FP(self):
"""Generate a global numbering for the dofs of an element in
:math:`\\text{FP}_{N-1}(\Omega)` on a crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
KN = K * N
FP_dofs_3d_x = np.arange((KN+1)*KN**2, dtype='int').reshape(
(KN+1, KN, KN), order='F')
accumulated = (KN + 1) * KN**2
FP_dofs_3d_y = np.arange(accumulated,
accumulated + (KN+1) * KN**2,
dtype='int').reshape((KN, KN+1, KN),
order='F')
accumulated = 2 * (KN + 1) * KN**2
FP_dofs_3d_z = np.arange(accumulated,
accumulated + (KN+1) * KN**2,
dtype='int').reshape((KN, KN, KN+1),
order='F')
GM_FP_x = np.zeros((K**3, (N+1) * N**2), dtype='int')
GM_FP_y = np.zeros((K**3, (N+1) * N**2), dtype='int')
GM_FP_z = np.zeros((K**3, (N+1) * N**2), dtype='int')
for k in range(K):
for j in range(K):
for i in range(K):
m = i + j * K + k * K ** 2
GM_FP_x[m] = FP_dofs_3d_x[
i*N:(i+1)*N+1,
j*N:(j+1)*N,
k*N:(k+1)*N].ravel('F')
GM_FP_y[m] = FP_dofs_3d_y[
i*N:(i+1)*N,
j*N:(j+1)*N+1,
k*N:(k+1)*N].ravel('F')
GM_FP_z[m] = FP_dofs_3d_z[
i*N:(i+1)*N,
j*N:(j+1)*N,
k*N:(k+1)*N+1].ravel('F')
GM_FP = np.hstack((GM_FP_x, GM_FP_y, GM_FP_z))
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 crazy mesh of
:math:`K^3` elements.
"""
K, N = self.K, self.N
KN = K * N
VP_dofs_3d = np.arange(KN ** 3, dtype='int'
).reshape((KN, KN, KN), order='F')
GM_VP = np.zeros((K ** 3, N ** 3), dtype='int')
for k in range(K):
for j in range(K):
for i in range(K):
m = i + j * K + k * K ** 2
GM_VP[m] = VP_dofs_3d[i*N:(i+1)*N,
j*N:(j+1)*N,
k*N:(k+1)*N].ravel('F')
return GM_VP
[docs]
class CrazyMeshGlobalBoundaryDOFs(object):
"""We find the global numbering of the dofs on each boundary of the
crazy mesh.
: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
:example:
>>> K = 2
>>> N = 1
>>> B_DOFs = CrazyMeshGlobalBoundaryDOFs(K, N)
>>> FB_dofs = B_DOFs.FP
>>> FB_dofs['x_minus'] # x=0 face
array([0, 3, 6, 9])
>>> FB_dofs['x_plus'] # x=1 face
array([ 2, 5, 8, 11])
>>> FB_dofs['y_minus'] # y=0 face
array([12, 13, 18, 19])
>>> FB_dofs['y_plus'] # y=1 face
array([16, 17, 22, 23])
>>> FB_dofs['z_minus'] # z=0 face
array([24, 25, 26, 27])
>>> FB_dofs['z_plus'] # z=1 face
array([32, 33, 34, 35])
"""
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.
"""
raise NotImplementedError("Could you complete this property?")
@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.
"""
raise NotImplementedError("Could you complete this property?")
@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
KN = K * N
FP_dofs_3d_x = np.arange((KN+1)*KN**2).reshape(
(KN+1, KN, KN), order='F')
accumulated = (KN + 1) * KN**2
FP_dofs_3d_y = np.arange(accumulated,
accumulated + (KN+1) * KN**2).reshape(
(KN, KN+1, KN), order='F')
accumulated = 2 * (KN + 1) * KN**2
FP_dofs_3d_z = np.arange(accumulated,
accumulated + (KN+1) * KN**2).reshape(
(KN, KN, KN+1), order='F')
DOFs_on_boundary = dict()
DOFs_on_boundary['x_minus'] = FP_dofs_3d_x[0, :, :].ravel('F')
DOFs_on_boundary['x_plus'] = FP_dofs_3d_x[-1, :, :].ravel('F')
DOFs_on_boundary['y_minus'] = FP_dofs_3d_y[:, 0, :].ravel('F')
DOFs_on_boundary['y_plus'] = FP_dofs_3d_y[:, -1, :].ravel('F')
DOFs_on_boundary['z_minus'] = FP_dofs_3d_z[:, :, 0].ravel('F')
DOFs_on_boundary['z_plus'] = FP_dofs_3d_z[:, :, -1].ravel('F')
return DOFs_on_boundary
if __name__ == '__main__':
import doctest
doctest.testmod()
def phi_exact(x, y, z):
return sin(2 * pi * x) * sin(2 * pi * y) * sin(2 * pi * z)
def u_exact(x, y, z):
return 2 * pi * cos(2 * pi * x) * sin(2 * pi * y) * sin(
2 * pi * z)
def v_exact(x, y, z):
return 2 * pi * sin(2 * pi * x) * cos(2 * pi * y) * sin(
2 * pi * z)
def w_exact(x, y, z):
return 2 * pi * sin(2 * pi * x) * sin(2 * pi * y) * cos(
2 * pi * z)
def f_exact(x, y, z):
return 12 * pi ** 2 * sin(2 * pi * x) * sin(2 * pi * y) * sin(
2 * pi * z)
from projection import Reduction
from L2_error import L2Error
from mimetic_basis_polynomials import MimeticBasisPolynomials
N = 2
KK = np.array([i for i in range(1,11)])
c=0.1
Error_p = list()
Error_u = list()
Error_w = list()
Error_f = list()
for K in KK:
_bfN_ = 'Lobatto-' + str(N)
mbf = MimeticBasisPolynomials(_bfN_, _bfN_, _bfN_)
crazy_mesh = CrazyMesh(c, K)
p_exact_local = list()
p_L2 = list()
w_exact_local = list()
w_L2 = list()
u_exact_local = list()
u_L2 = list()
f_exact_local = list()
f_L2 = list()
for k in range(K):
for j in range(K):
for i in range(K):
ct = crazy_mesh.CT_of_element_index(i, j, k)
RD = Reduction(mbf, ct)
L2e = L2Error(mbf, ct)
p_dofs_local = RD.NP(phi_exact)
p_L2_local = L2e.NP(p_dofs_local, phi_exact)
p_L2.append(p_L2_local ** 2)
w_dofs_local = RD.EP((u_exact, v_exact, w_exact))
w_L2_local = L2e.EP(w_dofs_local,
(u_exact, v_exact, w_exact))
w_L2.append(w_L2_local ** 2)
u_dofs_local = RD.FP((u_exact, v_exact, w_exact))
u_L2_local = L2e.FP(u_dofs_local,
(u_exact, v_exact, w_exact))
u_L2.append(u_L2_local ** 2)
f_dofs_local = RD.VP(f_exact)
f_L2_local = L2e.VP(f_dofs_local, f_exact)
f_L2.append(f_L2_local**2)
Error_p.append(np.sum(p_L2)**0.5)
Error_w.append(np.sum(w_L2)**0.5)
Error_u.append(np.sum(u_L2)**0.5)
Error_f.append(np.sum(f_L2)**0.5)
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(12, 12))
plt.subplot2grid((2, 2), (0, 0))
x = 1/KK
y = Error_p
order = (np.log10(y[-1])- np.log10(y[-2])) /\
(np.log10(x[-1])-np.log10(x[-2]))
plt.loglog(x, y,
label=f'Node, h-convergence, c={c}, order=%.2f'%order)
plt.legend()
plt.subplot2grid((2, 2), (0, 1))
x = 1/KK
y = Error_w
order = (np.log10(y[-1])- np.log10(y[-2])) /\
(np.log10(x[-1])-np.log10(x[-2]))
plt.loglog(x, y,
label=f'Edge, h-convergence, c={c}, order=%.2f'%order)
plt.legend()
plt.subplot2grid((2, 2), (1, 0))
x = 1/KK
y = Error_u
order = (np.log10(y[-1])- np.log10(y[-2])) /\
(np.log10(x[-1])-np.log10(x[-2]))
plt.loglog(x, y,
label=f'Face, h-convergence, c={c}, order=%.2f'%order)
plt.legend()
plt.subplot2grid((2, 2), (1, 1))
x = 1/KK
y = Error_f
order = (np.log10(y[-1])- np.log10(y[-2])) /\
(np.log10(x[-1])-np.log10(x[-2]))
plt.loglog(x, y,
label=f'Volume form, h-convergence, c={c}, order=%.2f'%order)
plt.legend()
plt.show()
