"""Mesh class definition (FVM mesh placeholder)"""
from typing import Any, Union, Tuple
import os
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
from cupy import ndarray
else:
import numpy as np
from numpy import ndarray
from g3d.defs import G3DWorld
from g3d.defs import DIMMode
from g3d.defs import FaceType
from g3d.flow.geo import Metrics
from g3d.mesh.geo import Coordinates
from g3d.mesh.geo import Grid
from g3d.mesh.geo import Face
from g3d.flow.geo import Q1DPipe
[docs]
class Mesh(Grid):
"""Mesh class definition, extends class Grid"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
self.__mtr = Metrics()
super().__init__(*args, **kwargs)
def _set_code_type(self):
""" """
# code type is set by solver env
pass
@property
def mtr(self) -> Metrics:
"""Grid metrics object
Returns
-------
Metrics
allocated instance of grid metrics object
"""
return self.__mtr
[docs]
def allocate_arrays(self) -> None:
super().allocate_arrays()
self.__mtr.allocate()
[docs]
def print(self) -> None:
# print mesh topology information and extreme values
super().print()
self.print_extreme_values()
[docs]
def print_extreme_values(self) -> None:
super().print_extreme_values()
self.__mtr.print_extreme_values()
print()
[docs]
def get_volume(self, idx: int = None) -> Union[ndarray, float]:
"""Gets volume for specific cell if cell index is specified otherwise all cell volumes are returned as an array
Parameters
----------
idx : int :
cell index, optional (Default value = None)
Returns
-------
float or ndarray
cell volume or array of cell volumes
"""
return self.__mtr.get_volume(idx)
[docs]
def get_inverse_volume(self, idx: int = None) -> Union[ndarray, float]:
"""Gets inverse volume for specific cell if cell index is specified otherwise all cell inverse volumes are
returned as an array
Parameters
----------
idx : int :
cell index, optional (Default value = None)
Returns
-------
float or ndarray
inverse cell volume or array of inverse cell volumes
"""
return self.__mtr.get_inverse_volume(idx)
[docs]
def get_area(self, idx: int = None) -> Union[ndarray, float]:
"""Gets average area for specific cell if cell index is specified otherwise all average cell areas are
returned as an array
Parameters
----------
idx : int :
cell index, optional (Default value = None)
Returns
-------
float or ndarray
average cell area or array of average cell areas
"""
return self.__mtr.get_area(idx)
[docs]
def get_streamtube_normal(self, idx: int = None) -> Union[ndarray, float]:
"""Gets stream-tube normal for specific cell if cell index is specified otherwise all stream-tube normals are
returned as an array
Parameters
----------
idx : int :
cell index, optional (Default value = None)
Returns
-------
float or ndarray
stream-tube normal or array of stream-tube normals
"""
return self.__mtr.get_streamtube_normal(idx)
[docs]
def check_face_def(self, f: Face, face_type: Any = FaceType.A_SIDE) -> None:
"""Checks definition of provided face
Parameters
----------
f : Face
mesh face
face_type : FaceType
face type (Default value = FaceType.A_SIDE)
"""
im, i0, j0, k0, idir, jdir, npi, npj = f.get_def()
if im < 1 or im > len(self._blocks):
f.check = 'block No. out of range'
return
ni, nj, nk = self._blocks[im - 1].get_block_sizes()
if i0 < 1 or i0 > ni or j0 < 1 or j0 > nj or k0 < 1 or k0 > nk:
f.check = 'boundary not contained in block'
return
if 1 < i0 < ni and 1 < j0 < nj and f.k0 > 1 and k0 < nk:
f.check = 'ref. corner point not on block boundary'
return
if idir == jdir:
f.check = 'local i, j-directions are identical'
return
icp = i0
jcp = j0
kcp = k0
if icp < 1 or icp > ni or jcp < 1 or jcp > nj or kcp < 1 or kcp > nk:
f.check = 'boundary not contained in block'
return
icp = i0 + self._id[idir - 1] * (npi - 1)
jcp = j0 + self._jd[idir - 1] * (npi - 1)
kcp = k0 + self._kd[idir - 1] * (npi - 1)
if icp < 1 or icp > ni or jcp < 1 or jcp > nj or kcp < 1 or kcp > nk:
f.check = 'boundary not contained in block'
return
icp = i0 + self._id[jdir - 1] * (npj - 1)
jcp = j0 + self._jd[jdir - 1] * (npj - 1)
kcp = k0 + self._kd[jdir - 1] * (npj - 1)
if icp < 1 or icp > ni or jcp < 1 or jcp > nj or kcp < 1 or kcp > nk:
f.check = 'boundary not contained in block'
return
icp = i0 + self._id[idir - 1] * (npi - 1) + self._id[jdir - 1] * (npj - 1)
jcp = j0 + self._jd[idir - 1] * (npi - 1) + self._jd[jdir - 1] * (npj - 1)
kcp = k0 + self._kd[idir - 1] * (npi - 1) + self._kd[jdir - 1] * (npj - 1)
if icp < 1 or icp > ni or jcp < 1 or jcp > nj or kcp < 1 or kcp > nk:
f.check = 'boundary not contained in block'
return
idkd = self._jd[idir - 1] * self._kd[jdir - 1] - self._kd[idir - 1] * self._jd[jdir - 1]
jdkd = self._kd[idir - 1] * self._id[jdir - 1] - self._id[idir - 1] * self._kd[jdir - 1]
kdkd = self._id[idir - 1] * self._jd[jdir - 1] - self._jd[idir - 1] * self._id[jdir - 1]
m = 1
if face_type is FaceType.B_SIDE:
m = -1
icp = i0 + idkd * m
jcp = j0 + jdkd * m
kcp = k0 + kdkd * m
if 1 <= icp <= ni and 1 <= jcp <= nj and 1 <= kcp <= nk:
f.check = 'face-orientation error detected'
return
f.check = 'face definition ok'
return
[docs]
def get_mesh_coordinates_from_pipe_data(self, q1d_pipe: Q1DPipe) -> None:
"""Generates mesh coordinates from provided quasi-one-dimensional pipe object
Parameters
----------
q1d_pipe : Q1DPipe :
instance of quasi-one-dimensional pipe object
Returns
-------
"""
x, y, z = self._nodes.get_XYZ()
q1d_pipe.extract_coordinates()
PX, PR = q1d_pipe.get_arrays()
sf = np.sqrt(np.pi)
for im in range(len(self._blocks)):
ni, nj, nk = self._blocks[im - 1].get_block_sizes()
for i in range(ni):
xp = PX[i]
rp = PR[i] * sf
x[i] = xp
y[i] = -0.5 * rp
z[i] = -0.5 * rp
x[i + ni] = xp
y[i + ni] = 0.5 * rp
z[i + ni] = -0.5 * rp
x[i + 2 * ni] = xp
y[i + 2 * ni] = -0.5 * rp
z[i + 2 * ni] = 0.5 * rp
x[i + 3 * ni] = xp
y[i + 3 * ni] = 0.5 * rp
z[i + 3 * ni] = 0.5 * rp
self.update_metrics()
def _import_coordinates_post(self) -> None:
# copy block coordinates to node array
x, y, z = self._nodes.get_XYZ()
for b in self._blocks:
l1 = b.nlow
l2 = b.nhigh
x[l1:l2], y[l1:l2], z[l1:l2] = b.get_XYZ()
# update metrics data
self.update_metrics()
[docs]
def update_metrics(self) -> None:
"""Updates metric data - cell volumes and face areas"""
g3d_world = G3DWorld()
if g3d_world.sizes.n_cells > 0:
x, y, z = self._nodes.get_XYZ()
NP = len(x)
if NP > 0:
self.__mtr.allocate()
vol = self.get_volume()
voli = self.get_inverse_volume()
area = self.get_area()
for im in range(len(self._blocks)):
ni, nj, nk = self._blocks[im].get_block_sizes()
i = np.arange(ni - 1)
j = np.arange(nj - 1)
k = np.arange(nk - 1)
i, j, k = np.meshgrid(i, j, k, indexing='ij')
l0 = self._blocks[im].get_node(i, j, k)
ic = self._blocks[im].get_cell(i, j, k)
if g3d_world.options.dim_mode is DIMMode.G2D_AXISYMM:
xp_0 = x[l0]
xp_1 = x[l0 + 1]
xp_2 = x[l0 + ni]
xp_3 = x[l0 + ni + 1]
yp_0 = y[l0]
yp_1 = y[l0 + 1]
yp_2 = y[l0 + ni]
yp_3 = y[l0 + ni + 1]
vol[ic] = ((xp_3 - xp_0) * (np.power(yp_2, 2) - np.power(yp_1, 2)) -
(xp_2 - xp_1) * (np.power(yp_3, 2) - np.power(yp_0, 2)) +
(xp_3 * yp_3 - xp_0 * yp_0) * (yp_2 - yp_1) -
(xp_2 * yp_2 - xp_1 * yp_1) * (yp_3 - yp_0)) / 6.
voli[ic] = 1. / vol[ic]
area[ic] = 0.5 * ((xp_3 - xp_0) * (yp_2 - yp_1) - (yp_3 - yp_0) * (xp_2 - xp_1))
else:
xp_1 = x[l0 + 1] - x[l0]
xp_2 = x[l0 + ni] - x[l0]
xp_3 = x[l0 + ni + 1] - x[l0]
xp_4 = x[l0 + ni * nj] - x[l0]
xp_5 = x[l0 + ni * nj + 1] - x[l0]
xp_6 = x[l0 + ni * nj + ni] - x[l0]
xp_7 = x[l0 + ni * nj + ni + 1] - x[l0]
yp_1 = y[l0 + 1] - y[l0]
yp_2 = y[l0 + ni] - y[l0]
yp_3 = y[l0 + ni + 1] - y[l0]
yp_4 = y[l0 + ni * nj] - y[l0]
yp_5 = y[l0 + ni * nj + 1] - y[l0]
yp_6 = y[l0 + ni * nj + ni] - y[l0]
yp_7 = y[l0 + ni * nj + ni + 1] - y[l0]
zp_1 = z[l0 + 1] - z[l0]
zp_2 = z[l0 + ni] - z[l0]
zp_3 = z[l0 + ni + 1] - z[l0]
zp_4 = z[l0 + ni * nj] - z[l0]
zp_5 = z[l0 + ni * nj + 1] - z[l0]
zp_6 = z[l0 + ni * nj + ni] - z[l0]
zp_7 = z[l0 + ni * nj + ni + 1] - z[l0]
vol[ic] = \
(xp_7 * yp_6 - xp_6 * yp_7) * (zp_4 + zp_5 - zp_2 - zp_3) + \
(xp_7 * yp_5 - xp_5 * yp_7) * (zp_1 + zp_3 - zp_4 - zp_6) + \
(xp_7 * yp_3 - xp_3 * yp_7) * (zp_2 + zp_6 - zp_1 - zp_5) + \
(xp_1 * yp_3 - xp_3 * yp_1 + xp_6 * yp_4 - xp_4 * yp_6) * (zp_5 + zp_7 - zp_2) + \
(xp_1 * yp_5 - xp_5 * yp_1 + xp_6 * yp_2 - xp_2 * yp_6) * (zp_4 - zp_3 - zp_7) + \
(xp_2 * yp_3 - xp_3 * yp_2 + xp_5 * yp_4 - xp_4 * yp_5) * (zp_1 - zp_6 - zp_7) + \
(xp_1 * yp_2 - xp_2 * yp_1) * zp_3 + \
(xp_1 * yp_4 - xp_4 * yp_1) * (-zp_5) + \
(xp_1 * yp_7 - xp_7 * yp_1) * (zp_5 - zp_3) + \
(xp_2 * yp_4 - xp_4 * yp_2) * zp_6 + \
(xp_2 * yp_7 - xp_7 * yp_2) * (zp_3 - zp_6) + \
(xp_3 * yp_5 - xp_5 * yp_3) * (zp_1 - zp_7) + \
(xp_3 * yp_6 - xp_6 * yp_3) * (zp_7 - zp_2) + \
(xp_4 * yp_7 - xp_7 * yp_4) * (zp_6 - zp_5) + \
(xp_5 * yp_6 - xp_6 * yp_5) * (zp_4 - zp_7)
vol[ic] /= 12.
voli[ic] = 1. / vol[ic]
if g3d_world.options.dim_mode is DIMMode.Q1D:
self.update_Q1D_metrics()
if g3d_world.options.dim_mode is DIMMode.Q3D:
self.update_Q3D_metrics()
self.__mtr.update_extreme_values()
self.print()
@staticmethod
def __calculate_surface_area(x: ndarray, y: ndarray, z: ndarray, pointers: list) -> float:
sx = 0.5 * ((y[pointers[3]] - y[pointers[0]]) * (z[pointers[2]] - z[pointers[1]])
- (z[pointers[3]] - z[pointers[0]]) * (y[pointers[2]] - y[pointers[1]]))
sy = 0.5 * ((z[pointers[3]] - z[pointers[0]]) * (x[pointers[2]] - x[pointers[1]])
- (x[pointers[3]] - x[pointers[0]]) * (z[pointers[2]] - z[pointers[1]]))
sz = 0.5 * ((x[pointers[3]] - x[pointers[0]]) * (y[pointers[2]] - y[pointers[1]])
- (y[pointers[3]] - y[pointers[0]]) * (x[pointers[2]] - x[pointers[1]]))
return (sx ** 2 + sy ** 2 + sz ** 2) ** 0.5
[docs]
def update_Q1D_metrics(self) -> None:
"""Updates quasi-one-dimensional metrics"""
g3d_world = G3DWorld()
if g3d_world.sizes.n_cells > 0:
x, y, z = self._nodes.get_XYZ()
NP = len(x)
if NP > 0:
area = self.get_area()
for im in range(len(self._blocks)):
ni, nj, nk = self._blocks[im].get_block_sizes()
for i in range(ni - 1):
pointers_1 = [i + 1, i + 1 + ni, i + 1 + 2 * ni, i + 1 + 3 * ni]
pointers_2 = [i, i + ni, i + 2 * ni, i + 3 * ni]
area[i] = (self.__calculate_surface_area(x, y, z, pointers_1) -
self.__calculate_surface_area(x, y, z, pointers_2))
[docs]
def update_Q3D_metrics(self) -> None:
"""Updates quasi-three-dimensional metrics"""
g3d_world = G3DWorld()
if g3d_world.sizes.n_cells > 0:
kp = [0] * 12
x, y, z = self._nodes.get_XYZ()
stn = self.get_streamtube_normal()
for im in range(len(self._blocks)):
ni, nj, nk = self._blocks[im].get_block_sizes()
if ni == 2:
i = 0
for j in range(nj - 1):
for k in range(nk - 1):
kc = self._blocks[im].get_cell(i, j, k)
lp0 = self._blocks[im].get_node(i, j, k)
kp[0] = lp0
kp[1] = lp0 + ni
kp[2] = lp0 + ni * nj
kp[3] = lp0 + ni + ni * nj
kp[4] = lp0 + 1
kp[5] = lp0 + 1 + ni
kp[6] = lp0 + 1 + ni * nj
kp[7] = lp0 + 1 + ni + ni * nj
self.calculate_streamtube_normal(stn, x, y, z, kp, kc)
elif nj == 2:
j = 0
for i in range(ni - 1):
for k in range(nk - 1):
kc = self._blocks[im].get_cell(i, j, k)
lp0 = self._blocks[im].get_node(i, j, k)
kp[0] = lp0
kp[1] = lp0 + ni * nj
kp[2] = lp0 + 1
kp[3] = lp0 + ni * nj + 1
kp[4] = lp0 + ni
kp[5] = lp0 + ni + ni * nj
kp[6] = lp0 + ni + 1
kp[7] = lp0 + ni + ni * nj + 1
self.calculate_streamtube_normal(stn, x, y, z, kp, kc)
elif nk == 2:
k = 0
for i in range(ni - 1):
for j in range(nj - 1):
kc = self._blocks[im].get_cell(i, j, k)
lp0 = self._blocks[im].get_node(i, j, k)
kp[0] = lp0
kp[1] = lp0 + 1
kp[2] = lp0 + ni
kp[3] = lp0 + ni + 1
kp[4] = lp0 + ni * nj
kp[5] = lp0 + ni * nj + 1
kp[6] = lp0 + ni + ni * nj
kp[7] = lp0 + ni + ni * nj + 1
self.calculate_streamtube_normal(stn, x, y, z, kp, kc)
[docs]
@staticmethod
def calculate_streamtube_normal(stn: ndarray, x: ndarray, y: ndarray, z: ndarray,
kp: list, kc: list) -> None:
"""Calculates stream-tube normals (used in quasi-three-dimensional mode)
Parameters
----------
stn : ndarray :
stream-tube normal array
x : ndarray :
x-coordinate array
y : ndarray :
y-coordinate array
z : ndarray :
z-coordinate array
kp : list :
node indices (list of ints)
kc : list :
cell indices (list of ints)
"""
sx = 0.25 * (((y[kp[3]] - y[kp[0]])
* (z[kp[2]] - z[kp[1]])
- (z[kp[3]] - z[kp[0]])
* (y[kp[2]] - y[kp[1]]))
+ ((y[kp[7]] - y[kp[4]])
* (z[kp[6]] - z[kp[5]])
- (z[kp[7]] - z[kp[4]])
* (y[kp[6]] - y[kp[5]])))
sy = 0.25 * (((z[kp[3]] - z[kp[0]])
* (x[kp[2]] - x[kp[1]])
- (x[kp[3]] - x[kp[0]])
* (z[kp[2]] - z[kp[1]]))
+ ((z[kp[7]] - z[kp[4]])
* (x[kp[6]] - x[kp[5]])
- (x[kp[7]] - x[kp[4]])
* (z[kp[6]] - z[kp[5]])))
sz = 0.25 * (((x[kp[3]] - x[kp[0]])
* (y[kp[2]] - y[kp[1]])
- (y[kp[3]] - y[kp[0]])
* (x[kp[2]] - x[kp[1]]))
+ ((x[kp[7]] - x[kp[4]])
* (y[kp[6]] - y[kp[5]])
- (y[kp[7]] - y[kp[4]])
* (x[kp[6]] - x[kp[5]])))
sl = np.sqrt(sx * sx + sy * sy + sz * sz)
stn[kc, 0] = sx / sl
stn[kc, 1] = sy / sl
stn[kc, 2] = sz / sl
# def check_bounds_G3D_classic(self, ip: int, jp: int, npi: int, npj: int, fm: FMP) -> Tuple[int, int, int, int]:
#
# mim = -1
# mip = 1
# mjm = -1
# mjp = 1
#
# if ip == 0:
# mim = 0
# if ip == npi - 2:
# mip = 0
# if jp == 0:
# mjm = 0
# if jp == npj - 2:
# mjp = 0
#
# fm.update_factors(mim, mip, mjm, mjp)
#
# return mim, mip, mjm, mjp
[docs]
def check_bounds(self, kc: int, ldi: int, ldj: int) -> Tuple[int, int, int, int]:
"""Check cell face bounds
Parameters
----------
kc : int :
cell index
ldi : int :
array index increment for one step in local i-direction
ldj : int :
array index increment for one step in local j-direction
Returns
-------
tuple
mim, mip, mjm, mjp
"""
mim = -1
mip = 1
mjm = -1
mjp = 1
for b in self._blocks:
if b.clow <= kc < b.chigh:
mim, mip = b.check_face_block_bounds(kc, ldi)
mjm, mjp = b.check_face_block_bounds(kc, ldj)
return mim, mip, mjm, mjp
[docs]
def check_bounds_block_interface(self, kca: int,
ldia: int,
ldja: int,
kcb: int,
ldib: int,
ldjb: int) -> Tuple[int, int, int, int]:
"""Checks block interface bounds
Parameters
----------
kca : int
cell index A-side
ldia : int
index increment for one step in local i-direction (A-side)
ldja : int
index increment for one step in local j-direction (A-side)
kcb : int
cell index B-side
ldib : int
index increment for one step in local i-direction (B-side)
ldjb : int
index increment for one step in local j-direction (B-side)
Returns
-------
tuple
mim, mip, mjm, mjp
"""
mim_a, mip_a, mjm_a, mjp_a = self.check_bounds(kca, ldia, ldja)
mim_b, mip_b, mjm_b, mjp_b = self.check_bounds(kcb, ldib, ldjb)
mim = -mim_a * mim_b
mip = mip_a * mip_b
mjm = -mjm_a * mjm_b
mjp = mjp_a * mjp_b
return mim, mip, mjm, mjp
[docs]
def get_cell_centers(self) -> Tuple[ndarray, ...]:
"""Get cell centers for all cells
Returns
-------
tuple
cell center coordinate vectors (x, y, z)
"""
g3d_world = G3DWorld()
xp, yp, zp = self._nodes.get_XYZ()
cell_centers = Coordinates(tag='cell_coordinates', size=g3d_world.sizes.n_cells)
xc, yc, zc = cell_centers.get_XYZ()
for b in self._blocks:
ni, nj, nk = b.get_block_sizes()
for k in range(nk - 1):
for j in range(nj - 1):
for i in range(ni - 1):
lc = b.get_cell(i, j, k)
lp1 = b.get_node(i, j, k)
lp2 = b.get_node(i + 1, j, k)
lp3 = b.get_node(i, j + 1, k)
lp4 = b.get_node(i + 1, j + 1, k)
lp5 = b.get_node(i, j, k + 1)
lp6 = b.get_node(i + 1, j, k + 1)
lp7 = b.get_node(i, j + 1, k + 1)
lp8 = b.get_node(i + 1, j + 1, k + 1)
xc[lc] = 0.125 * (xp[lp1] + xp[lp2] + xp[lp3] + xp[lp4] + xp[lp5] + xp[lp6] + xp[lp7] + xp[lp8])
yc[lc] = 0.125 * (yp[lp1] + yp[lp2] + yp[lp3] + yp[lp4] + yp[lp5] + yp[lp6] + yp[lp7] + yp[lp8])
zc[lc] = 0.125 * (zp[lp1] + zp[lp2] + zp[lp3] + zp[lp4] + zp[lp5] + zp[lp6] + zp[lp7] + zp[lp8])
return xc, yc, zc