"""mesh building block - Block class definition"""
from typing import Union, Any, Tuple, BinaryIO
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
import struct
from g3d.defs import G3DWorld
from g3d.defs import G3DDisplay
from g3d.defs import FaceType
from g3d.defs import FaceOrientation
from g3d.defs import BlockFace
from g3d.defs import AngleType
from g3d.defs import Axis
from g3d.defs import NodeDistribution
from g3d.defs import FunctionByPass
from g3d.mesh.geo import Coordinates
from g3d.mesh.geo import Coord, Edge, Face
[docs]
class Block(Coordinates):
"""class Block (mesh block), extends class Coordinates"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
self.__block_is_defined = False
self.__reset_def()
self._pos = 0
self._filename = "msh.bin"
self.__edges = []
self.__coordinates = []
self.__faces = []
self.__faces_sn = []
self._id = [1, -1, 0, 0, 0, 0]
self._jd = [0, 0, 1, -1, 0, 0]
self._kd = [0, 0, 0, 0, 1, -1]
self._im = 0
self._ni = 0
self._nj = 0
self._nk = 0
self._nmp = 0
self._nmc = 0
self._np = 0
self._nc = 0
self._size = 0
self.__DX = np.zeros(0)
self.__DY = np.zeros(0)
self.__DZ = np.zeros(0)
super().__init__(*args, **kwargs)
if not self.__block_is_defined:
self.__allocate()
def _set_from_args(self, *args: Any, **kwargs: Any) -> ndarray:
bps = self._get_function_evaluate_status(**kwargs)
if bps is FunctionByPass.SELF:
args_stat = np.ones(len(args))
else:
args_stat = super()._set_from_args(*args, bps=FunctionByPass.BYPASS)
if bps in [FunctionByPass.EVALUATE, FunctionByPass.SELF]:
if np.sum(args_stat) > 0:
block_def = []
for i in range(len(args)):
if args_stat[i] == 1:
arg = args[i]
if type(arg) is Block:
# if caught by a super class, the child class instance of copy will be used
args_stat[i] = 0
self.copy(arg)
self.__block_is_defined = True
elif type(arg) in [int, np.int64, np.int32]:
args_stat[i] = 0
if type(arg) in [np.int64, np.int32]:
block_def.append(arg.item())
else:
block_def.append(arg)
elif type(arg) is list:
args_stat[i] = 0
if self._list_is_all_ints(arg):
self._nk = 2
if len(arg) > 0:
self._ni = arg[0]
if len(arg) > 1:
self._nj = arg[1]
if len(arg) > 2:
self._nk = arg[2]
if len(block_def) > 0:
# will be caught by the for-loop below in a sub call to _set_from_args
self._set_from_args(block_def, bps=FunctionByPass.SELF)
return args_stat
def _set_default(self) -> None:
super()._set_default()
self._tag = 'g3d_block'
def __reset_def(self) -> None:
self._im = 0
self._ni = 0
self._nj = 0
self._nk = 0
self._nmp = 0
self._nmc = 0
self.update_counts()
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
"""
sets code parameters from specified key => value data
input: key, value
valid keys for Block class:
'filename', 'pos', 'ni', 'nj', 'nk'
example: set_from_key_value(key, value)
"""
if super().set_from_key_value(key, value):
return True
elif key == 'filename':
self._filename = value
return True
elif key == 'pos':
self._pos = value
return True
elif key == 'ni':
self._ni = value
return True
elif key == 'nj':
self._nj = value
return True
elif key == 'nk':
self._nk = value
return True
return False
[docs]
def allocate_XYZ(self, N: int = None) -> None:
"""
allocate XYZ vectors
input: vector length (optional input, default value = np)
example: allocate_XYZ(N)
"""
self._size = self._np
super().allocate_XYZ(N)
@property
def im(self) -> int:
"""mesh-block number"""
return self._im
@im.setter
def im(self, i: int) -> None:
self._im = i
self.__update_im()
@property
def ni(self) -> int:
"""number of nodes in the i-direction"""
return self._ni
@ni.setter
def ni(self, i: int) -> None:
self._ni = i
self.update_counts()
@property
def nj(self) -> int:
"""number of nodes in the j-direction"""
return self._nj
@nj.setter
def nj(self, i: int) -> None:
self._nj = i
self.update_counts()
@property
def nk(self) -> int:
"""number of nodes in the k-direction"""
return self._nk
@nk.setter
def nk(self, i: int) -> None:
self._nk = i
self.update_counts()
@property
def nmp(self) -> int:
"""accumulated node count in blocks with lower block number"""
return self._nmp
@nmp.setter
def nmp(self, i: int) -> None:
self._nmp = i
self.__update_node_range()
@property
def nmc(self) -> int:
"""accumulated cell count in blocks with lower block number"""
return self._nmc
@nmc.setter
def nmc(self, i: int) -> None:
self._nmc = i
self.__update_cell_range()
@property
def pos(self) -> int:
"""file position (file record pointer)"""
return self._pos
@pos.setter
def pos(self, p: int) -> None:
self._pos = p
@property
def filename(self) -> str:
"""mesh file filename"""
return self._filename
@filename.setter
def filename(self, f: str) -> None:
self._filename = f
@property
def n_cells(self) -> int:
"""number of cells in block"""
return self._nc
@property
def n_nodes(self) -> int:
"""number of nodes in block"""
return self._np
@property
def clow(self) -> int:
"""global cell index of first cell in block"""
return self._clow
@property
def chigh(self) -> int:
"""global cell index of last cell in block"""
return self._chigh
@property
def nlow(self) -> int:
"""global node index of first node in block"""
return self._nlow
@property
def nhigh(self) -> int:
"""global node index of last node in block"""
return self._nhigh
@property
def edges(self) -> list:
"""list of block edges"""
return self.__edges
@property
def coordinates(self) -> list:
"""block node coordinates"""
return self.__coordinates
[docs]
def update_counts(self) -> None:
"""updates node and cell counts"""
self._np = self._ni * self._nj * self._nk
self._nc = (self._ni - 1) * (self._nj - 1) * (self._nk - 1)
self.__update_ranges()
def __update_ranges(self) -> None:
self.__update_cell_range()
self.__update_node_range()
def __update_cell_range(self) -> None:
self._clow = self._nmc
self._chigh = self._clow + self._nc
def __update_node_range(self) -> None:
self._nlow = self._nmp
self._nhigh = self._nlow + self._np
[docs]
def get_block_sizes(self) -> Tuple[int, int, int]:
"""
return block node counts
input: no input
example: ni, nj, nk = get_block_sizes()
"""
return self._ni, self._nj, self._nk
def __get_lp(self, i: int, j: int, k: int) -> int: # local 0-based indices
return i + self._ni * j + self._ni * self._nj * k
def __get_lc(self, i: int, j: int, k: int) -> int: # local 0-based indices
return i + (self._ni - 1) * j + (self._ni - 1) * (self._nj - 1) * k
[docs]
def get_node(self, i: int, j: int, k: int) -> int: # global 0-based indices
"""
returns global node index
input: i, j, k (0-based indices)
example: l = get_node(i, j, k)
"""
return self.__get_lp(i, j, k) + self._nmp
[docs]
def get_cell(self, i: int, j: int, k: int) -> int: # global 0-based indices
"""
returns global cell index
input: i, j, k (0-based indices)
example: l = get_cell(i, j, k)
"""
return self.__get_lc(i, j, k) + self._nmc
[docs]
def get_index(self, i: int, j: int, k: int) -> int: # local 1-based indices
"""
returns local node index
input: i, j, k (1-based indices)
example:
"""
return (i - 1) + self._ni * (j - 1) + self._ni * self._nj * (k - 1)
[docs]
def cell_get_ijk(self, kc: int) -> Tuple[int, int, int]: # local 0-based indices
"""
get 0-based indices from global cell index
input: cell index
example i, j, k = cell_get_ijk(kc)
"""
kcl = kc - self._nmc
k = int(kcl / ((self._ni - 1) * (self._nj - 1)))
kcl -= k * ((self._ni - 1) * (self._nj - 1))
j = int(kcl / (self._ni - 1))
kcl -= j * (self._ni - 1)
i = kcl
return i, j, k
[docs]
def node_get_ijk(self, kp: int) -> Tuple[int, int, int]: # local 0-based indices
"""
get 0-based indices from global node index
input: node index
example i, j, k = node_get_ijk(kp)
"""
kpl = kp - self._nmp
k = int(kpl / (self._ni * self._nj))
kpl -= k * (self._ni * self._nj)
j = int(kpl / self._ni)
kpl -= j * self._ni
i = kpl
return i, j, k
[docs]
def get_def(self) -> Tuple[int, int, int, int]:
"""
returns block definition
input: no input
example: im, ni, nj, nk = get_def()
"""
return self._im, self._ni, self._nj, self._nk
[docs]
def check_face_block_bounds(self, kc: int, ld: int) -> Tuple[int, int]:
"""
check if face increment is within block boundaries
input: cell index (kc), cell increment (ld)
example: check_face_block_bounds(kc, ld)
"""
mm = -1
mp = 1
ndi = 1
ndj = (self._ni - 1)
ndk = (self._ni - 1) * (self._nj - 1)
i, j, k = self.cell_get_ijk(kc)
if ld == ndi or ld == -ndi:
# i - direction
if ld > 0:
if i == 0:
mm = 0
if i == self._ni - 2:
mp = 0
else:
if i == 0:
mp = 0
if i == self._ni - 2:
mm = 0
elif ld == ndj or ld == -ndj:
# j - direction
if ld > 0:
if j == 0:
mm = 0
if j == self._nj - 2:
mp = 0
else:
if j == 0:
mp = 0
if j == self._nj - 2:
mm = 0
elif ld == ndk or ld == -ndk:
# k - direction
if ld > 0:
if k == 0:
mm = 0
if k == self._nk - 2:
mp = 0
else:
if k == 0:
mp = 0
if k == self._nk - 2:
mm = 0
return mm, mp
[docs]
def print(self) -> None:
"""
prints block definition, block cell and node count
"""
g3d_disp = G3DDisplay()
g3d_disp.header('block: %s' % str(self._im))
super().print()
print(' ni, nj, nk = [%d, %d, %d]' % (self._ni, self._nj, self._nk))
print(' number of nodes: ' + str(self._np))
print(' number of cells: ' + str(self._nc))
[docs]
def copy(self, b: 'Block') -> None:
"""
copies block definition from specified block
input: block
example: copy(Block b)
"""
self.allocate(b.ni, b.nj, b.nk)
super().copy(b)
[docs]
def allocate(self, *args: Any, **kwargs: Any) -> None:
"""
allocates coordinate vectors
input: args, kwargs
args: ni, nj, nk
kwargs: ni=, nj=, nk=
"""
self.set(*args, **kwargs)
self.__allocate()
def __allocate(self) -> None:
if self._ni > 0 and self._nj > 0:
if self._nk == 0:
self._nk = 2
self.update_counts()
self.allocate_XYZ()
self.init_Face_list()
self.reset_smooth_delta()
[docs]
def allocate_smooth_delta(self) -> None:
"""allocate coordinate vectors used for mesh smooth operator"""
self.__DX = np.zeros(self._np)
self.__DY = np.zeros(self._np)
self.__DZ = np.zeros(self._np)
[docs]
def reset_smooth_delta(self) -> None:
"""reset coordinate vectors used for mesh smooth operator"""
self.__DX = np.zeros(0)
self.__DY = np.zeros(0)
self.__DZ = np.zeros(0)
[docs]
def init_Face_list(self) -> None:
"""initialize face lists based on block definition"""
self.__faces = []
self.__faces_sn = []
f1 = Face([self._im, 1, 1, 1, 5, 3, self._nk, self._nj])
f2 = Face([self._im, self._ni, 1, 1, 3, 5, self._nj, self._nk])
f3 = Face([self._im, 1, 1, 1, 1, 5, self._ni, self._nk])
f4 = Face([self._im, 1, self._nj, 1, 5, 1, self._nk, self._ni])
f5 = Face([self._im, 1, 1, 1, 3, 1, self._nj, self._ni])
f6 = Face([self._im, 1, 1, self._nk, 1, 3, self._ni, self._nj])
self.__faces.append(f1)
self.__faces.append(f2)
self.__faces.append(f3)
self.__faces.append(f4)
self.__faces.append(f5)
self.__faces.append(f6)
f11 = Face(f1)
f12 = Face(f2)
f13 = Face(f3)
f14 = Face(f4)
f15 = Face(f5)
f16 = Face(f6)
f11.swap_normal()
f12.swap_normal()
f13.swap_normal()
f14.swap_normal()
f15.swap_normal()
f16.swap_normal()
self.__faces_sn.append(f11)
self.__faces_sn.append(f12)
self.__faces_sn.append(f13)
self.__faces_sn.append(f14)
self.__faces_sn.append(f15)
self.__faces_sn.append(f16)
def __nullify_deltas(self) -> None:
self.__DX[:] = 0.
self.__DY[:] = 0.
self.__DZ[:] = 0.
def __update_im(self) -> None:
for f in self.__faces:
f.im = self._im
for f in self.__faces_sn:
f.im = self._im
for e in self.__edges:
e.im = self._im
[docs]
def add_Edge(self, e: Edge, edge_def: list = None) -> None:
"""
add edge to block edge list
input: Edge (e), edge definition in block reference (optional input)
"""
e.im = self._im
e.applied = False
if edge_def is not None:
if len(edge_def) == 5:
e.define_Edge(edge_def)
self.__edges.append(e)
self.__apply_edges()
[docs]
def add_Coord(self, c: Coord) -> None:
"""
add coordinate to block coordinate list
input: Coordinate (c)
"""
self.__coordinates.append(c)
def __apply_edges(self, force_apply: bool = False) -> None:
if force_apply:
for e in self.__edges:
e.applied = False
for e in self.__edges:
if not e.applied:
i0A, j0A, k0A, idirA, npiA = e.get_def()
ip = np.arange(npiA, dtype=np.int32)
lp = (i0A + self._id[idirA - 1] * ip - 1 +
self._ni * (j0A + self._jd[idirA - 1] * ip - 1) +
self._ni * self._nj * (k0A + self._kd[idirA - 1] * ip - 1))
self._X[lp], self._Y[lp], self._Z[lp] = e.get_XYZ()
e.applied = True
[docs]
def scale(self, sf: float) -> None:
"""
scale block coordinates
input: scale factor (sf)
example: scale(sf)
"""
super().scale(sf)
for e in self.__edges:
e.scale(sf)
for c in self.__coordinates:
c.scale(sf)
[docs]
def get_cell_pointers(self) -> ndarray:
"""
returns block range of cell pointers
input: no input
example: lc = get_cell_pointers()
"""
ll = np.zeros(self._nc, dtype=np.int32)
ip = np.arange((self._ni - 1), dtype=np.int32)
for k in range(self._nk - 1):
for j in range(self._nj - 1):
l1 = (self._ni - 1) * j + (self._ni - 1) * (self._nj - 1) * k
l2 = l1 + (self._ni - 1)
ll[l1:l2] = ip + self._ni * j + self._ni * self._nj * k
return ll
[docs]
def get_node_pointers(self) -> ndarray:
"""
returns block range of node pointers
input: no input
example: lp = get_node_pointers()
"""
ll = np.zeros(self._np, dtype=np.int32)
ip = np.arange(self._ni, dtype=np.int32)
for k in range(self._nk):
for j in range(self._nj):
l1 = self._ni * j + self._ni * self._nj * k
l2 = l1 + self._ni
ll[l1:l2] = ip + self._ni * j + self._ni * self._nj * k
return ll
[docs]
def check_volumes(self) -> Tuple[int, float, float]:
"""
checks volumes for all cells in block
returns number of negative volumes, min volume, and max volume
warning message printed in case negative volumes are found
input: no input
example: n_neg, v_min, v_max = check_volumes()
"""
g3d_disp = G3DDisplay()
V = self.__get_volume(self.get_cell_pointers())
v_min = np.amin(V)
v_max = np.amax(V)
n_neg = 0
if v_min < 0.:
if not g3d_disp.silent:
self.print()
for ic in range(self._nc):
if V[ic] < 0.:
n_neg += 1
i, j, k = self.cell_get_ijk(ic)
if not g3d_disp.silent:
print(' Warning! negative volume found for cell index: %d (i=%d, j=%d, k=%d)' % (ic, i, j, k))
if not g3d_disp.silent:
g3d_disp.hline('-')
print(' cell volume check report:')
g3d_disp.hline('-')
if n_neg > 0:
print(' warning! %d negative volumes found' % n_neg)
else:
print(' all cells are ok - no negative volumes found')
print(' cell volume min/max: %e/%e' % (v_min, v_max))
g3d_disp.hline()
return n_neg, v_min, v_max
def __get_volume(self, lp: Union[int, ndarray, list[int]]) -> Union[float, ndarray, list[float]]:
X1 = self._X[lp]
Y1 = self._Y[lp]
Z1 = self._Z[lp]
l2 = lp + 1
X2 = self._X[l2] - X1
Y2 = self._Y[l2] - Y1
Z2 = self._Z[l2] - Z1
l3 = lp + self._ni
X3 = self._X[l3] - X1
Y3 = self._Y[l3] - Y1
Z3 = self._Z[l3] - Z1
l4 = l3 + 1
X4 = self._X[l4] - X1
Y4 = self._Y[l4] - Y1
Z4 = self._Z[l4] - Z1
l5 = lp + self._ni * self._nj
X5 = self._X[l5] - X1
Y5 = self._Y[l5] - Y1
Z5 = self._Z[l5] - Z1
l6 = l5 + 1
X6 = self._X[l6] - X1
Y6 = self._Y[l6] - Y1
Z6 = self._Z[l6] - Z1
l7 = l5 + self._ni
X7 = self._X[l7] - X1
Y7 = self._Y[l7] - Y1
Z7 = self._Z[l7] - Z1
l8 = l7 + 1
X8 = self._X[l8] - X1
Y8 = self._Y[l8] - Y1
Z8 = self._Z[l8] - Z1
X1 = 0.
Y1 = 0.
Z1 = 0.
return ((X1 * Y2 - X2 * Y1 + X8 * Y7 - X7 * Y8) * (Z5 + Z6 - Z3 - Z4)
+ (X1 * Y3 - X3 * Y1 + X8 * Y6 - X6 * Y8) * (Z2 + Z4 - Z5 - Z7)
+ (X1 * Y5 - X5 * Y1 + X8 * Y4 - X4 * Y8) * (Z3 + Z7 - Z2 - Z6)
+ (X2 * Y4 - X4 * Y2 + X7 * Y5 - X5 * Y7) * (Z6 + Z8 - Z1 - Z3)
+ (X2 * Y6 - X6 * Y2 + X7 * Y3 - X3 * Y7) * (Z1 + Z5 - Z4 - Z8)
+ (X3 * Y4 - X4 * Y3 + X6 * Y5 - X5 * Y6) * (Z1 + Z2 - Z7 - Z8)
+ (X1 * Y4 - X4 * Y1) * (Z2 - Z3)
+ (X1 * Y6 - X6 * Y1) * (Z5 - Z2)
+ (X1 * Y7 - X7 * Y1) * (Z3 - Z5)
+ (X2 * Y3 - X3 * Y2) * (Z4 - Z1)
+ (X2 * Y5 - X5 * Y2) * (Z1 - Z6)
+ (X2 * Y8 - X8 * Y2) * (Z6 - Z4)
+ (X3 * Y5 - X5 * Y3) * (Z7 - Z1)
+ (X3 * Y8 - X8 * Y3) * (Z4 - Z7)
+ (X4 * Y6 - X6 * Y4) * (Z2 - Z8)
+ (X4 * Y7 - X7 * Y4) * (Z8 - Z3)
+ (X5 * Y8 - X8 * Y5) * (Z7 - Z6)
+ (X6 * Y7 - X7 * Y6) * (Z5 - Z8)) / 12.
[docs]
def get_Coord_from_IJK(self, i: int, j: int, k: int, **kwargs: Any) -> Coord:
"""
get Coord from (i, j, k)
input: indices (i, j ,k)
example: Coord c = get_Coord_from_IJK(i, j, k)
"""
x, y, z = self.get_XYZ_from_IJK(i, j, k)
c = Coord(x, y, z, **kwargs)
return c
[docs]
def get_XYZ_from_IJK(self, i: int, j: int, k: int) -> Tuple[float, float, float]:
"""
get (x, y, z) from (i, j, k)
input: indices (i, j ,k)
example: x, y, z = get_XYZ_from_IJK(i, j, k)
"""
lp = self.get_index(i, j, k)
return self._X[lp], self._Y[lp], self._Z[lp]
[docs]
def get_X_from_IJK(self, i: int, j: int, k: int) -> float:
"""
get x from (i, j, k)
input: indices (i, j ,k)
example: x = get_X_from_IJK(i, j, k)
"""
lp = self.get_index(i, j, k)
return self._X[lp]
[docs]
def get_Y_from_IJK(self, i: int, j: int, k: int) -> float:
"""
get y from (i, j, k)
input: indices (i, j ,k)
example: y = get_Y_from_IJK(i, j, k)
"""
lp = self.get_index(i, j, k)
return self._Y[lp]
[docs]
def get_Z_from_IJK(self, i: int, j: int, k: int) -> float:
"""
get z from (i, j, k)
input: indices (i, j ,k)
example: z = get_Z_from_IJK(i, j, k)
"""
lp = self.get_index(i, j, k)
return self._Z[lp]
[docs]
def remove_Coord_by_tag(self, tag: str) -> None:
"""
remove Coord with specified tag from coordinate list
input: tag
example: remove_Coord_by_tag(tag)
"""
tmp = []
for c in self.__coordinates:
if c.tag != tag:
tmp.append(c)
self.__coordinates.clear()
for c in tmp:
self.__coordinates.append(c)
[docs]
def get_Coord_by_tag(self, tag: str) -> Union[Coord, None]:
"""
get Coord with specified tag from coordinate list
input: tag
example: Coord c = get_Coord_by_tag(tag)
"""
for c in self.__coordinates:
if c.tag == tag:
return c
return None
[docs]
def get_Coord(self, lp: int, **kwargs: Any) -> Coord:
"""
get Coord from coordinate vector index
input: vector index (lp)
optional input: kwargs (forwarded to Coord constructor)
example: Coord c = get_Coord(lp, kwargs)
"""
x, y, z = self.get_XYZ(lp)
c = Coord(x, y, z, **kwargs)
return c
[docs]
def get_delta_arrays(self) -> Tuple[ndarray, ...]:
"""
returns delta coordinates (grid smooth placeholders)
input: no input
example: dx, dy, dz = get_delta_arrays()
"""
return self.__DX, self.__DY, self.__DZ
[docs]
def import_XYZ(self, f: BinaryIO) -> None:
"""
reads (x, y, z) coordinates from file
input: file handle (f)
example: import_XYZ(f)
"""
for k in range(self._nk):
for j in range(self._nj):
for i in range(self._ni):
lp = self.__get_lp(i, j, k)
self._X[lp] = struct.unpack('d', f.read(8))[0]
self._Y[lp] = struct.unpack('d', f.read(8))[0]
self._Z[lp] = struct.unpack('d', f.read(8))[0]
[docs]
def import_XY(self, f: BinaryIO) -> None:
"""
reads (x, y) coordinates from file
input: file handle (f)
example: import_XY(f)
"""
for j in range(self._nj):
for i in range(self._ni):
lp1 = self.__get_lp(i, j, 0)
lp2 = self.__get_lp(i, j, 1)
self._X[lp1] = struct.unpack('d', f.read(8))[0]
self._Y[lp1] = struct.unpack('d', f.read(8))[0]
self._Z[lp1] = 0.
self._X[lp2] = self._X[lp1]
self._Y[lp2] = self._Y[lp1]
self._Z[lp2] = 1.
[docs]
def export_XYZ(self, f: BinaryIO) -> None:
"""
writes (x, y, z) coordinates to file
input: file handle (f)
example: export_XYZ(f)
"""
for k in range(self._nk):
for j in range(self._nj):
for i in range(self._ni):
lp = self.__get_lp(i, j, k)
f.write(self._X[lp].tobytes())
f.write(self._Y[lp].tobytes())
f.write(self._Z[lp].tobytes())
[docs]
def export_XY(self, f: BinaryIO) -> None:
"""
writes (x, y) coordinates to file
input: file handle (f)
example: export_XY(f)
"""
for j in range(self._nj):
for i in range(self._ni):
lp = self.__get_lp(i, j, 0)
f.write(self._X[lp].tobytes())
f.write(self._Y[lp].tobytes())
[docs]
def remove_Edge_by_tag(self, tag: str) -> None:
"""
remove Edge with specified tag from Edge list
input: tag
example: remove_Edge_by_tag(tag)
"""
tmp = []
for e in self.__edges:
if e.tag != tag:
tmp.append(e)
self.__edges.clear()
for e in tmp:
self.__edges.append(e)
[docs]
def get_Edge_by_tag(self, tag: str) -> Union[Edge, None]:
"""
get Edge with specified tag from Edge list
input: tag
example: EDge e = get_Edge_by_tag(tag)
"""
for e in self.__edges:
if e.tag == tag:
return e
return None
[docs]
def get_Edge(self, i0: int, j0: int, k0: int, idir: int, npi: int) -> Edge:
"""
get Bock Edge from specified Edge definition
input: Edge definition (i0, j0, k0, idir, npi)
(i0, j0, k0) = start index
idir = direction
npi = number of nodes
example: EDge e = get_Edge_by_tag(tag)
"""
X = np.zeros(npi)
Y = np.zeros(npi)
Z = np.zeros(npi)
ip = np.arange(npi, dtype=np.int32)
lp = (i0 + self._id[idir - 1] * ip - 1 +
self._ni * (j0 + self._jd[idir - 1] * ip - 1) +
self._ni * self._nj * (k0 + self._kd[idir - 1] * ip - 1))
X[:] = self._X[lp]
Y[:] = self._Y[lp]
Z[:] = self._Z[lp]
c1 = self.get_Coord(lp[0])
c2 = self.get_Coord(lp[-1])
e = Edge([i0, j0, k0, idir, npi], c1, c2)
e.set_XYZ(X, Y, Z)
return e
[docs]
def populate_face_coordinate_placeholders(self, face: Face) -> None:
"""
populates face node placeholders
copies (x, y, z) coordinate data from Block coordinate arrays
to face coordinate vectors
input: no input
"""
if type(face) is not Face:
print('the specified input is not a face')
return
face.allocate_XYZ()
X, Y, Z = face.get_XYZ()
im, i0, j0, k0, idir, jdir, npi, npj = face.get_def()
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lp = self.__get_lp(ip, jp, kp)
face_node = face.get_node(i, j)
X[face_node] = self._X[lp]
Y[face_node] = self._Y[lp]
Z[face_node] = self._Z[lp]
[docs]
def apply_face_coordinates(self, face: Face) -> None:
"""
copies face coordinate vector data to Block coordinate placeholders
input: Face f
"""
if type(face) is not Face:
print('the specified input is not a face')
return
X, Y, Z = face.get_XYZ()
im, i0, j0, k0, idir, jdir, npi, npj = face.get_def()
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lp = self.__get_lp(ip, jp, kp)
face_node = face.get_node(i, j)
self._X[lp] = X[face_node]
self._Y[lp] = Y[face_node]
self._Z[lp] = Z[face_node]
[docs]
def face_offset(self, original_face: Face, offset_dist: float, offset_dir: int, offset_n: int) -> None:
"""
generates a Face from specified original Face and offset information
applies offset Face to Block coordinate vectors
input: Face (original face), offset distance (distance in face-normal direction),
offset direction, number of nodes in offset direction
definition of offset direction:
1 (+i)
2 (-i)
3 (+j)
4 (-j)
5 (+k)
6 (-k)
example: face_offset(Face, dist, dir, n_nodes)
"""
self.populate_face_coordinate_placeholders(original_face)
offset_face = original_face.generate_offset_Face(offset_dist)
offset_face.i0 = offset_face.i0 + self._id[offset_dir - 1] * offset_n
offset_face.j0 = offset_face.j0 + self._jd[offset_dir - 1] * offset_n
offset_face.k0 = offset_face.k0 + self._kd[offset_dir - 1] * offset_n
self.apply_face_coordinates(offset_face)
[docs]
def get_Face(self, N: Any, normal: Any = FaceOrientation.OUTWARD_FACING_NORMAL) -> Union[Face, None]:
"""
gets Block Face:
input: BlockFace (Enum), normal (Enum, default value = FaceOrientation.OUTWARD_FACING_NORMAL)
valid options:
Blockface: ILOW, IHIGH, JLOW, JHIGH, KLOW, KHIGH
normal: OUTWARD_FACING_NORMAL, INWARD_FACING_NORMAL
example: Face f = get_Face(BlockFace, normal)
"""
if type(N) is BlockFace:
N = int(N.value - 1)
if normal == FaceOrientation.INWARD_FACING_NORMAL:
return self.__faces_sn[N]
else:
return self.__faces[N]
return None
def __process_face_def(self, face_def: Union[list, BlockFace, Face]) -> Tuple[int, ...]:
if type(face_def) is Face:
face_def = face_def.get_def()
elif type(face_def) is BlockFace:
face_def = (self.get_Face(face_def)).get_def()
if len(face_def) == 7:
i0, j0, k0, idir, jdir, npi, npj = face_def
else:
bim, i0, j0, k0, idir, jdir, npi, npj = face_def
return i0, j0, k0, idir, jdir, npi, npj
[docs]
def get_face_pointers(self, face_def: Union[list, BlockFace, Face]) -> Tuple[int, ...]:
"""
returns Face data from specified input
input: Face or BlockFace, or Face definition
example: i0, j0, k0, idir, jdir, npi, npj = get_face_pointers(f)
(i0, j0, k0) = start index
(idir, jdir) = local i-direction, local j-direction
(npi, npj) = number of nodes in local i-direction and local j-direction
"""
return self.__process_face_def(face_def)
[docs]
def print_XYZ(self, face_def: Union[list, BlockFace, Face]) -> None:
"""
prints (x, y, z) coordinates for nodes on specified face
input: Face (given as Face, BlockFace, or face definition)
"""
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
for j in range(npj):
print(' j = ' + str(j + 1))
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lp = self.__get_lp(ip, jp, kp)
print('(i, j) = (' + str(i + 1) + ', ' + str(j + 1) + ') : [x, y, z] = ['
+ str(self._X[lp]) + ', ' + str(self._Y[lp]) + ', ' + str(self._Z[lp]) + ']')
[docs]
def scatter_face_XYZ(self, plt: Any, face_def: Union[list, BlockFace, Face], m: str = 'o') -> None:
"""
generate scatter plot for (x, y, z) coordinates for nodes on specified face
input: Face (given as Face, BlockFace, or face definition)
optional input: scatter symbol (default value = 'o')
"""
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
XM = np.zeros(npi * npj)
YM = np.zeros(npi * npj)
ZM = np.zeros(npi * npj)
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lp = self.__get_lp(ip, jp, kp)
ll = i + npi * j
XM[ll] = self._X[lp]
YM[ll] = self._Y[lp]
ZM[ll] = self._Z[lp]
plt.scatter(XM, YM, marker=m)
[docs]
def shift_Face(self, face_def: Union[list, BlockFace, Face], **kwargs: Any) -> None:
"""
:param face_def:
:param kwargs:
:return:
"""
x_shift = 0.
y_shift = 0.
z_shift = 0.
if len(kwargs) != 0:
g3d_world = G3DWorld()
for key, value in kwargs.items():
if g3d_world.display.verbose:
print('%s == %s' % (key, str(value)))
if key == 'x_shift':
x_shift = value
elif key == 'y_shift':
y_shift = value
elif key == 'z_shift':
z_shift = value
else:
print('warning! parameter %s is not defined' % key)
if x_shift == 0. and y_shift == 0. and z_shift == 0.:
print('warning! no shift specified')
return
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lp = self.__get_lp(ip, jp, kp)
self._X[lp] += x_shift
self._Y[lp] += y_shift
self._Z[lp] += z_shift
[docs]
def shift(self, **kwargs: Any) -> None:
"""
:param kwargs:
:return:
"""
x_shift = 0.
y_shift = 0.
z_shift = 0.
if len(kwargs) != 0:
g3d_world = G3DWorld()
for key, value in kwargs.items():
if g3d_world.display.verbose:
print('%s == %s' % (key, str(value)))
if key == 'x_shift':
x_shift = value
elif key == 'y_shift':
y_shift = value
elif key == 'z_shift':
z_shift = value
else:
print('warning! parameter %s is not defined' % key)
if x_shift == 0. and y_shift == 0. and z_shift == 0.:
print('warning! no shift specified')
return
for lp in range(self._np):
self._X[lp] += x_shift
self._Y[lp] += y_shift
self._Z[lp] += z_shift
[docs]
def rotate(self, **kwargs: Any) -> None:
"""
:param kwargs:
:return:
"""
alpha = 0.
beta = 0.
gamma = 0.
angle_dim = AngleType.DEGREES
if len(kwargs) != 0:
g3d_world = G3DWorld()
for key, value in kwargs.items():
if g3d_world.display.verbose:
print('%s == %s' % (key, str(value)))
if key == 'alpha':
alpha = value
elif key == 'beta':
beta = value
elif key == 'gamma':
gamma = value
elif key == 'angle_dim':
angle_dim = value
else:
print('warning! parameter %s is not defined' % key)
if alpha == 0. and beta == 0. and gamma == 0.:
print('warning! no rotation angle specified')
return
if angle_dim == AngleType.DEGREES:
alpha *= np.pi / 180.
beta *= np.pi / 180.
gamma *= np.pi / 180.
eps = 1.e-6
csa = np.cos(alpha)
ssa = np.sin(alpha)
csb = np.cos(beta)
ssb = np.sin(beta)
csg = np.cos(gamma)
ssg = np.sin(gamma)
for lp in range(self._np):
if np.abs(alpha) > eps:
yy = self._Y[lp]
zz = self._Z[lp]
self._Y[lp] = csa * yy - ssa * zz
self._Z[lp] = ssa * yy + csa * zz
if np.abs(beta) > eps:
xx = self._X[lp]
zz = self._Z[lp]
self._Z[lp] = csb * zz + ssb * xx
self._X[lp] = -ssb * zz + csb * xx
if np.abs(gamma) > eps:
xx = self._X[lp]
yy = self._Y[lp]
self._X[lp] = csg * xx - ssg * yy
self._Y[lp] = ssg * xx + csg * yy
[docs]
def copy_rotate_Edge(self, face_def: Union[list, BlockFace, Face], alpha: float, **kwargs: Any) -> None:
"""
:param face_def:
:param alpha:
:param kwargs:
"""
angle_dim = AngleType.DEGREES
axis = Axis.X_AXIS
if len(kwargs) != 0:
g3d_world = G3DWorld()
for key, value in kwargs.items():
if g3d_world.display.verbose:
print('%s == %s' % (key, str(value)))
if key == 'axis':
axis = value
elif key == 'angle_dim':
angle_dim = value
else:
print('warning! parameter %s is not defined' % key)
if angle_dim == AngleType.DEGREES:
alpha *= np.pi / 180.
self.__apply_edges()
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
d_alfa = alpha / (npj - 1)
cos_a = np.cos(d_alfa)
sin_a = np.sin(d_alfa)
if axis == Axis.X_AXIS:
for j in range(1, npj):
for i in range(npi):
ip1 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * (j - 1)
jp1 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * (j - 1)
kp1 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * (j - 1)
ip2 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp2 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp2 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
l1 = self.__get_lp(ip1, jp1, kp1)
l2 = self.__get_lp(ip2, jp2, kp2)
self._X[l2] = self._X[l1]
yy = self._Y[l1]
zz = self._Z[l1]
self._Y[l2] = cos_a * yy - sin_a * zz
self._Z[l2] = sin_a * yy + cos_a * zz
elif axis == Axis.Y_AXIS:
for j in range(1, npj):
for i in range(npi):
ip1 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * (j - 1)
jp1 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * (j - 1)
kp1 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * (j - 1)
ip2 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp2 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp2 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
l1 = self.__get_lp(ip1, jp1, kp1)
l2 = self.__get_lp(ip2, jp2, kp2)
self._Y[l2] = self._Y[l1]
xx = self._X[l1]
zz = self._Z[l1]
self._Z[l2] = cos_a * zz + sin_a * xx
self._X[l2] = -sin_a * zz + cos_a * xx
elif axis == Axis.Z_AXIS:
for j in range(1, npj):
for i in range(npi):
ip1 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * (j - 1)
jp1 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * (j - 1)
kp1 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * (j - 1)
ip2 = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp2 = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp2 = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
l1 = self.__get_lp(ip1, jp1, kp1)
l2 = self.__get_lp(ip2, jp2, kp2)
self._Z[l2] = self._Z[l1]
xx = self._X[l1]
yy = self._Y[l1]
self._X[l2] = cos_a * xx - sin_a * yy
self._Y[l2] = sin_a * xx + cos_a * yy
[docs]
def copy_rotate_Face(self, face_def: Union[list, BlockFace, Face], npk: int, alpha: float, **kwargs: Any) -> None:
"""
:param face_def:
:param npk:
:param alpha:
:param kwargs:
"""
angle_dim = AngleType.DEGREES
axis = Axis.X_AXIS
if len(kwargs) != 0:
g3d_world = G3DWorld()
for key, value in kwargs.items():
if g3d_world.display.verbose:
print('%s == %s' % (key, str(value)))
if key == 'axis':
axis = value
elif key == 'angle_dim':
angle_dim = value
else:
print('warning! parameter %s is not defined' % key)
if angle_dim == AngleType.DEGREES:
alpha *= np.pi / 180.
self.__apply_edges()
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
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]
ldka = idkd + self._ni * jdkd + self._ni * self._nj * kdkd
# kdir = 1
# if (ldka == -1):
# kdir = 2
# elif (np.abs(ldka) == self._ni):
# kdir = 3
# if (ldka < 0):
# kdir = 4
# elif (np.abs(ldka) == self._ni * self._nj):
# kdir = 5
# if (ldka < 0):
# kdir = 6
lc_face = np.zeros(npi * npj, dtype=np.uint32)
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lc_face[i + j * npi] = self.__get_lp(ip, jp, kp)
d_alfa = alpha / (npk - 1)
cos_a = np.cos(d_alfa)
sin_a = np.sin(d_alfa)
if axis == Axis.X_AXIS:
for k in range(1, npk):
l1 = lc_face + ldka * (k - 1)
l2 = lc_face + ldka * k
self._X[l2] = self._X[l1]
yy = self._Y[l1]
zz = self._Z[l1]
self._Y[l2] = cos_a * yy - sin_a * zz
self._Z[l2] = sin_a * yy + cos_a * zz
elif axis == Axis.Y_AXIS:
for k in range(1, npk):
l1 = lc_face + ldka * (k - 1)
l2 = lc_face + ldka * k
self._Y[l2] = self._Y[l1]
xx = self._X[l1]
zz = self._Z[l1]
self._Z[l2] = cos_a * zz + sin_a * xx
self._X[l2] = -sin_a * zz + cos_a * xx
elif axis == Axis.Z_AXIS:
for k in range(1, npk):
l1 = lc_face + ldka * (k - 1)
l2 = lc_face + ldka * k
self._Z[l2] = self._Z[l1]
xx = self._X[l1]
yy = self._Y[l1]
self._X[l2] = cos_a * xx - sin_a * yy
self._Y[l2] = sin_a * xx + cos_a * yy
[docs]
def extrude_Face(self, face_def: Union[list, BlockFace, Face], edge: Edge) -> None:
"""
:param face_def:
:param edge:
"""
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
npk = edge.npi
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]
ldka = idkd + self._ni * jdkd + self._ni * self._nj * kdkd
kdir = 1
if ldka == -1:
kdir = 2
elif np.abs(ldka) == self._ni:
kdir = 3
if ldka < 0:
kdir = 4
elif np.abs(ldka) == self._ni * self._nj:
kdir = 5
if ldka < 0:
kdir = 6
lc_face = np.zeros(npi * npj, dtype=np.uint32)
for j in range(npj):
for i in range(npi):
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lc_face[i + j * npi] = self.__get_lp(ip, jp, kp)
for k in range(npk):
lc = lc_face + ldka * k
self.copy_XYZ(lc_face, lc)
if kdir < 3:
self._X[lc] = edge.X[k]
elif kdir < 5:
self._Y[lc] = edge.Y[k]
else:
self._Z[lc] = edge.Z[k]
[docs]
def fill_Face(self, face_def: Union[list, BlockFace, Face]) -> Tuple[ndarray, ...]:
"""
:param face_def:
:return:
"""
self.__apply_edges()
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
u = np.zeros(npi)
for i in range(1, npi):
lu = np.zeros(4, dtype=np.uint32)
for jjd in range(2):
for iid in range(2):
iu = (i0 - 1) + self._id[idir - 1] * (i - iid) + self._id[jdir - 1] * (npj - 1) * jjd
ju = (j0 - 1) + self._jd[idir - 1] * (i - iid) + self._jd[jdir - 1] * (npj - 1) * jjd
ku = (k0 - 1) + self._kd[idir - 1] * (i - iid) + self._kd[jdir - 1] * (npj - 1) * jjd
lu[iid + 2 * jjd] = self.__get_lp(iu, ju, ku)
u[i] = u[i - 1] + float(self.__get_weight_2D(lu))
u /= u[-1]
v = np.zeros(npj)
for j in range(1, npj):
lv = np.zeros(4, dtype=np.uint32)
for iid in range(2):
for jjd in range(2):
iv = (i0 - 1) + self._id[idir - 1] * (npi - 1) * iid + self._id[jdir - 1] * (j - jjd)
jv = (j0 - 1) + self._jd[idir - 1] * (npi - 1) * iid + self._jd[jdir - 1] * (j - jjd)
kv = (k0 - 1) + self._kd[idir - 1] * (npi - 1) * iid + self._kd[jdir - 1] * (j - jjd)
lv[jjd + 2 * iid] = self.__get_lp(iv, jv, kv)
v[j] = v[j - 1] + float(self.__get_weight_2D(lv))
v /= v[-1]
for j in range(1, npj - 1):
b1 = 1 - v[j]
b2 = v[j]
for i in range(1, npi - 1):
a1 = 1 - u[i]
a2 = u[i]
lc = np.zeros(9, dtype=np.uint32)
ip = (i0 - 1)
jp = (j0 - 1)
kp = (k0 - 1)
lc[0] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * i
jp = (j0 - 1) + self._jd[idir - 1] * i
kp = (k0 - 1) + self._kd[idir - 1] * i
lc[1] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * (npi - 1)
jp = (j0 - 1) + self._jd[idir - 1] * (npi - 1)
kp = (k0 - 1) + self._kd[idir - 1] * (npi - 1)
lc[2] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[jdir - 1] * j
lc[3] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * j
lc[4] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * (npi - 1) + self._id[jdir - 1] * j
jp = (j0 - 1) + self._jd[idir - 1] * (npi - 1) + self._jd[jdir - 1] * j
kp = (k0 - 1) + self._kd[idir - 1] * (npi - 1) + self._kd[jdir - 1] * j
lc[5] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[jdir - 1] * (npj - 1)
jp = (j0 - 1) + self._jd[jdir - 1] * (npj - 1)
kp = (k0 - 1) + self._kd[jdir - 1] * (npj - 1)
lc[6] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * i + self._id[jdir - 1] * (npj - 1)
jp = (j0 - 1) + self._jd[idir - 1] * i + self._jd[jdir - 1] * (npj - 1)
kp = (k0 - 1) + self._kd[idir - 1] * i + self._kd[jdir - 1] * (npj - 1)
lc[7] = self.__get_lp(ip, jp, kp)
ip = (i0 - 1) + self._id[idir - 1] * (npi - 1) + self._id[jdir - 1] * (npj - 1)
jp = (j0 - 1) + self._jd[idir - 1] * (npi - 1) + self._jd[jdir - 1] * (npj - 1)
kp = (k0 - 1) + self._kd[idir - 1] * (npi - 1) + self._kd[jdir - 1] * (npj - 1)
lc[8] = self.__get_lp(ip, jp, kp)
self.__interp_2D(lc, a1, a2, b1, b2)
return u, v
[docs]
def fill(self, vol_def: list = None) -> None:
"""
:param vol_def:
"""
i1, i2, j1, j2, k1, k2 = [0, self._ni - 1, 0, self._nj - 1, 0, self._nk - 1]
if vol_def is not None:
if len(vol_def) == 6:
vol_def = np.array(vol_def)
i1, i2, j1, j2, k1, k2 = vol_def - 1
npi = i2 - i1 + 1
u = np.zeros(npi)
for i in range(1, npi):
ip = i + i1
lu = np.zeros(8, dtype=np.uint32)
lu[0] = self.__get_lp(ip, j1, k1)
lu[1] = lu[0] - 1
lu[2] = self.__get_lp(ip, j2, k1)
lu[3] = lu[2] - 1
lu[4] = self.__get_lp(ip, j1, k2)
lu[5] = lu[4] - 1
lu[6] = self.__get_lp(ip, j2, k2)
lu[7] = lu[6] - 1
u[i] = u[i - 1] + float(self.__get_weight_3D(lu))
u /= u[-1]
npj = j2 - j1 + 1
v = np.zeros(npj)
for j in range(1, npj):
jp = j + j1
lv = np.zeros(8, dtype=np.uint32)
lv[0] = self.__get_lp(i1, jp, k1)
lv[1] = lv[0] - self._ni
lv[2] = self.__get_lp(i2, jp, k1)
lv[3] = lv[2] - self._ni
lv[4] = self.__get_lp(i1, jp, k2)
lv[5] = lv[4] - self._ni
lv[6] = self.__get_lp(i2, jp, k2)
lv[7] = lv[6] - self._ni
v[j] = v[j - 1] + float(self.__get_weight_3D(lv))
v /= v[-1]
npk = k2 - k1 + 1
w = np.zeros(npk)
for k in range(1, npk):
kp = k + k1
lw = np.zeros(8, dtype=np.uint32)
lw[0] = self.__get_lp(i1, j1, kp)
lw[1] = lw[0] - self._ni * self._nj
lw[2] = self.__get_lp(i2, j1, kp)
lw[3] = lw[2] - self._ni * self._nj
lw[4] = self.__get_lp(i1, j2, kp)
lw[5] = lw[4] - self._ni * self._nj
lw[6] = self.__get_lp(i2, j2, kp)
lw[7] = lw[6] - self._ni * self._nj
w[k] = w[k - 1] + float(self.__get_weight_3D(lw))
w /= w[-1]
for k in range(1, npk - 1):
c1 = 1 - w[k]
c2 = w[k]
kp = k + k1
for j in range(1, npj - 1):
b1 = 1 - v[j]
b2 = v[j]
jp = j + j1
for i in range(1, npi - 1):
a1 = 1 - u[i]
a2 = u[i]
ip = i + i1
lc = np.zeros(27, dtype=np.uint32)
lc[0] = self.__get_lp(i1, j1, k1)
lc[1] = self.__get_lp(ip, j1, k1)
lc[2] = self.__get_lp(i2, j1, k1)
lc[3] = self.__get_lp(i1, jp, k1)
lc[4] = self.__get_lp(ip, jp, k1)
lc[5] = self.__get_lp(i2, jp, k1)
lc[6] = self.__get_lp(i1, j2, k1)
lc[7] = self.__get_lp(ip, j2, k1)
lc[8] = self.__get_lp(i2, j2, k1)
lc[9] = self.__get_lp(i1, j1, kp)
lc[10] = self.__get_lp(ip, j1, kp)
lc[11] = self.__get_lp(i2, j1, kp)
lc[12] = self.__get_lp(i1, jp, kp)
lc[13] = self.__get_lp(ip, jp, kp)
lc[14] = self.__get_lp(i2, jp, kp)
lc[15] = self.__get_lp(i1, j2, kp)
lc[16] = self.__get_lp(ip, j2, kp)
lc[17] = self.__get_lp(i2, j2, kp)
lc[18] = self.__get_lp(i1, j1, k2)
lc[19] = self.__get_lp(ip, j1, k2)
lc[20] = self.__get_lp(i2, j1, k2)
lc[21] = self.__get_lp(i1, jp, k2)
lc[22] = self.__get_lp(ip, jp, k2)
lc[23] = self.__get_lp(i2, jp, k2)
lc[24] = self.__get_lp(i1, j2, k2)
lc[25] = self.__get_lp(ip, j2, k2)
lc[26] = self.__get_lp(i2, j2, k2)
self.__interp_3D(lc, a1, a2, b1, b2, c1, c2)
def __get_weight_2D(self, ll: ndarray) -> float:
a = np.array([self._X[ll[0]], self._Y[ll[0]], self._Z[ll[0]]])
b = np.array([self._X[ll[1]], self._Y[ll[1]], self._Z[ll[1]]])
c = np.array([self._X[ll[2]], self._Y[ll[2]], self._Z[ll[2]]])
d = np.array([self._X[ll[3]], self._Y[ll[3]], self._Z[ll[3]]])
ab = float(np.linalg.norm(a - b))
cd = float(np.linalg.norm(c - d))
return ab + cd
# return np.linalg.norm(np.array([self._X[ll[0]], self._Y[ll[0]], self._Z[ll[0]]]) -
# np.array([self._X[ll[1]], self._Y[ll[1]], self._Z[ll[1]]])) + \
# np.linalg.norm(np.array([self._X[ll[2]], self._Y[ll[2]], self._Z[ll[2]]]) -
# np.array([self._X[ll[3]], self._Y[ll[3]], self._Z[ll[3]]]))
def __interp_2D(self, lc: ndarray, a1: float, a2: float, b1: float, b2: float) -> None:
self._X[lc[4]] = (a1 * self._X[lc[3]] + a2 * self._X[lc[5]]
+ b1 * self._X[lc[1]] + b2 * self._X[lc[7]]
- a1 * b1 * self._X[lc[0]] - a2 * b1 * self._X[lc[2]]
- a1 * b2 * self._X[lc[6]] - a2 * b2 * self._X[lc[8]])
self._Y[lc[4]] = (a1 * self._Y[lc[3]] + a2 * self._Y[lc[5]]
+ b1 * self._Y[lc[1]] + b2 * self._Y[lc[7]]
- a1 * b1 * self._Y[lc[0]] - a2 * b1 * self._Y[lc[2]]
- a1 * b2 * self._Y[lc[6]] - a2 * b2 * self._Y[lc[8]])
self._Z[lc[4]] = (a1 * self._Z[lc[3]] + a2 * self._Z[lc[5]]
+ b1 * self._Z[lc[1]] + b2 * self._Z[lc[7]]
- a1 * b1 * self._Z[lc[0]] - a2 * b1 * self._Z[lc[2]]
- a1 * b2 * self._Z[lc[6]] - a2 * b2 * self._Z[lc[8]])
def __get_weight_3D(self, ll: ndarray) -> float:
a = np.array([self._X[ll[0]], self._Y[ll[0]], self._Z[ll[0]]])
b = np.array([self._X[ll[1]], self._Y[ll[1]], self._Z[ll[1]]])
c = np.array([self._X[ll[2]], self._Y[ll[2]], self._Z[ll[2]]])
d = np.array([self._X[ll[3]], self._Y[ll[3]], self._Z[ll[3]]])
e = np.array([self._X[ll[4]], self._Y[ll[4]], self._Z[ll[4]]])
f = np.array([self._X[ll[5]], self._Y[ll[5]], self._Z[ll[5]]])
g = np.array([self._X[ll[6]], self._Y[ll[6]], self._Z[ll[6]]])
h = np.array([self._X[ll[7]], self._Y[ll[7]], self._Z[ll[7]]])
ab = float(np.linalg.norm(a - b))
cd = float(np.linalg.norm(c - d))
ef = float(np.linalg.norm(e - f))
gh = float(np.linalg.norm(g - h))
return ab + cd + ef + gh
# return np.linalg.norm(np.array([self._X[ll[0]], self._Y[ll[0]], self._Z[ll[0]]]) -
# np.array([self._X[ll[1]], self._Y[ll[1]], self._Z[ll[1]]])) + \
# np.linalg.norm(np.array([self._X[ll[2]], self._Y[ll[2]], self._Z[ll[2]]]) -
# np.array([self._X[ll[3]], self._Y[ll[3]], self._Z[ll[3]]])) + \
# np.linalg.norm(np.array([self._X[ll[4]], self._Y[ll[4]], self._Z[ll[4]]]) -
# np.array([self._X[ll[5]], self._Y[ll[5]], self._Z[ll[5]]])) + \
# np.linalg.norm(np.array([self._X[ll[6]], self._Y[ll[6]], self._Z[ll[6]]]) -
# np.array([self._X[ll[7]], self._Y[ll[7]], self._Z[ll[7]]]))
def __interp_3D(self, lc: ndarray, a1: float, a2: float, b1: float, b2: float, c1: float, c2: float) -> None:
self._X[lc[13]] = (a1 * self._X[lc[12]] + a2 * self._X[lc[14]]
+ b1 * self._X[lc[10]] + b2 * self._X[lc[16]]
+ c1 * self._X[lc[4]] + c2 * self._X[lc[22]]
- a1 * b1 * self._X[lc[9]] - a2 * b1 * self._X[lc[11]]
- a1 * b2 * self._X[lc[15]] - a2 * b2 * self._X[lc[17]]
- a1 * c1 * self._X[lc[3]] - a2 * c1 * self._X[lc[5]]
- a1 * c2 * self._X[lc[21]] - a2 * c2 * self._X[lc[23]]
- b1 * c1 * self._X[lc[1]] - b2 * c1 * self._X[lc[7]]
- b1 * c2 * self._X[lc[19]] - b2 * c2 * self._X[lc[25]]
+ a1 * b1 * c1 * self._X[lc[0]] + a2 * b1 * c1 * self._X[lc[2]]
+ a1 * b2 * c1 * self._X[lc[6]] + a2 * b2 * c1 * self._X[lc[8]]
+ a1 * b1 * c2 * self._X[lc[18]] + a2 * b1 * c2 * self._X[lc[20]]
+ a1 * b2 * c2 * self._X[lc[24]] + a2 * b2 * c2 * self._X[lc[26]])
self._Y[lc[13]] = (a1 * self._Y[lc[12]] + a2 * self._Y[lc[14]]
+ b1 * self._Y[lc[10]] + b2 * self._Y[lc[16]]
+ c1 * self._Y[lc[4]] + c2 * self._Y[lc[22]]
- a1 * b1 * self._Y[lc[9]] - a2 * b1 * self._Y[lc[11]]
- a1 * b2 * self._Y[lc[15]] - a2 * b2 * self._Y[lc[17]]
- a1 * c1 * self._Y[lc[3]] - a2 * c1 * self._Y[lc[5]]
- a1 * c2 * self._Y[lc[21]] - a2 * c2 * self._Y[lc[23]]
- b1 * c1 * self._Y[lc[1]] - b2 * c1 * self._Y[lc[7]]
- b1 * c2 * self._Y[lc[19]] - b2 * c2 * self._Y[lc[25]]
+ a1 * b1 * c1 * self._Y[lc[0]] + a2 * b1 * c1 * self._Y[lc[2]]
+ a1 * b2 * c1 * self._Y[lc[6]] + a2 * b2 * c1 * self._Y[lc[8]]
+ a1 * b1 * c2 * self._Y[lc[18]] + a2 * b1 * c2 * self._Y[lc[20]]
+ a1 * b2 * c2 * self._Y[lc[24]] + a2 * b2 * c2 * self._Y[lc[26]])
self._Z[lc[13]] = (a1 * self._Z[lc[12]] + a2 * self._Z[lc[14]]
+ b1 * self._Z[lc[10]] + b2 * self._Z[lc[16]]
+ c1 * self._Z[lc[4]] + c2 * self._Z[lc[22]]
- a1 * b1 * self._Z[lc[9]] - a2 * b1 * self._Z[lc[11]]
- a1 * b2 * self._Z[lc[15]] - a2 * b2 * self._Z[lc[17]]
- a1 * c1 * self._Z[lc[3]] - a2 * c1 * self._Z[lc[5]]
- a1 * c2 * self._Z[lc[21]] - a2 * c2 * self._Z[lc[23]]
- b1 * c1 * self._Z[lc[1]] - b2 * c1 * self._Z[lc[7]]
- b1 * c2 * self._Z[lc[19]] - b2 * c2 * self._Z[lc[25]]
+ a1 * b1 * c1 * self._Z[lc[0]] + a2 * b1 * c1 * self._Z[lc[2]]
+ a1 * b2 * c1 * self._Z[lc[6]] + a2 * b2 * c1 * self._Z[lc[8]]
+ a1 * b1 * c2 * self._Z[lc[18]] + a2 * b1 * c2 * self._Z[lc[20]]
+ a1 * b2 * c2 * self._Z[lc[24]] + a2 * b2 * c2 * self._Z[lc[26]])
[docs]
def fill_2D(self, w: float = 1.0) -> None:
"""
:param w:
"""
c1 = self.get_Coord_from_IJK(1, 1, 1)
c2 = Coord(c1)
c1.z = 0.0
c2.z = w
e_ext = Edge(self._nk, c1, c2, distribute=NodeDistribution.CURVE_2)
self.extrude_Face(self.get_Face(BlockFace.KLOW, FaceOrientation.INWARD_FACING_NORMAL), e_ext)
[docs]
def smooth_Face(self, face_def: Union[list, BlockFace, Face]) -> None:
"""
:param face_def:
"""
self.__apply_edges()
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
l0 = self.__get_lp(i0 - 1, j0 - 1, k0 - 1)
ldi = self._id[idir - 1] + self._ni * self._jd[idir - 1] + self._ni * self._nj * self._kd[idir - 1]
ldj = self._id[jdir - 1] + self._ni * self._jd[jdir - 1] + self._ni * self._nj * self._kd[jdir - 1]
ii, jj = np.meshgrid(np.arange(1, npi - 1), np.arange(1, npj - 1))
l11 = l0 + ldi * (ii - 1) + ldj * (jj - 1)
l21 = l0 + ldi * ii + ldj * (jj - 1)
l31 = l0 + ldi * (ii + 1) + ldj * (jj - 1)
l12 = l0 + ldi * (ii - 1) + ldj * jj
l22 = l0 + ldi * ii + ldj * jj
l32 = l0 + ldi * (ii + 1) + ldj * jj
l13 = l0 + ldi * (ii - 1) + ldj * (jj + 1)
l23 = l0 + ldi * ii + ldj * (jj + 1)
l33 = l0 + ldi * (ii + 1) + ldj * (jj + 1)
XI = .5 * (self._X[l32] - self._X[l12])
YI = .5 * (self._Y[l32] - self._Y[l12])
ZI = .5 * (self._Z[l32] - self._Z[l12])
XJ = .5 * (self._X[l23] - self._X[l21])
YJ = .5 * (self._Y[l23] - self._Y[l21])
ZJ = .5 * (self._Z[l23] - self._Z[l21])
alpha = XJ ** 2 + YJ ** 2 + ZJ ** 2
beta = -2. * (XI * XJ + YI * YJ + ZI * ZJ)
gamma = XI ** 2 + YI ** 2 + ZI ** 2
s_max = 4. * alpha + np.abs(beta) + 4. * gamma
scf = np.where(s_max != 0, 1. / s_max, 1.)
XII = self._X[l32] - 2. * self._X[l22] + self._X[l12]
YII = self._Y[l32] - 2. * self._Y[l22] + self._Y[l12]
ZII = self._Z[l32] - 2. * self._Z[l22] + self._Z[l12]
XJJ = self._X[l23] - 2. * self._X[l22] + self._X[l21]
YJJ = self._Y[l23] - 2. * self._Y[l22] + self._Y[l21]
ZJJ = self._Z[l23] - 2. * self._Z[l22] + self._Z[l21]
XIJ = .25 * (self._X[l33] - self._X[l13] + self._X[l11] - self._X[l31])
YIJ = .25 * (self._Y[l33] - self._Y[l13] + self._Y[l11] - self._Y[l31])
ZIJ = .25 * (self._Z[l33] - self._Z[l13] + self._Z[l11] - self._Z[l31])
self.__DX[l22] = scf * (alpha * XII + beta * XIJ + gamma * XJJ)
self.__DY[l22] = scf * (alpha * YII + beta * YIJ + gamma * YJJ)
self.__DZ[l22] = scf * (alpha * ZII + beta * ZIJ + gamma * ZJJ)
# for j in range(1, npj - 1):
# for i in range(1, npi - 1):
# l11 = l0 + ldi * (i - 1) + ldj * (j - 1)
# l21 = l0 + ldi * (i) + ldj * (j - 1)
# l31 = l0 + ldi * (i + 1) + ldj * (j - 1)
# l12 = l0 + ldi * (i - 1) + ldj * (j)
# l22 = l0 + ldi * (i) + ldj * (j)
# l32 = l0 + ldi * (i + 1) + ldj * (j)
# l13 = l0 + ldi * (i - 1) + ldj * (j + 1)
# l23 = l0 + ldi * (i) + ldj * (j + 1)
# l33 = l0 + ldi * (i + 1) + ldj * (j + 1)
# XI = .5 * (self._X[l32] - self._X[l12])
# YI = .5 * (self._Y[l32] - self._Y[l12])
# ZI = .5 * (self._Z[l32] - self._Z[l12])
# XJ = .5 * (self._X[l23] - self._X[l21])
# YJ = .5 * (self._Y[l23] - self._Y[l21])
# ZJ = .5 * (self._Z[l23] - self._Z[l21])
# alpha = XJ**2 + YJ**2 + ZJ**2
# beta = -2. * (XI * XJ + YI * YJ + ZI * ZJ)
# gamma = XI**2 + YI**2 + ZI**2
# s_max = 4. * alpha + np.abs(beta) + 4. * gamma
# scf = 1.
# if (s_max != 0.):
# scf /= s_max
# XII = self._X[l32] - 2.*self._X[l22] + self._X[l12]
# YII = self._Y[l32] - 2.*self._Y[l22] + self._Y[l12]
# ZII = self._Z[l32] - 2.*self._Z[l22] + self._Z[l12]
# XJJ = self._X[l23] - 2.*self._X[l22] + self._X[l21]
# YJJ = self._Y[l23] - 2.*self._Y[l22] + self._Y[l21]
# ZJJ = self._Z[l23] - 2.*self._Z[l22] + self._Z[l21]
# XIJ = .25 * (self._X[l33] - self._X[l13] + self._X[l11] - self._X[l31])
# YIJ = .25 * (self._Y[l33] - self._Y[l13] + self._Y[l11] - self._Y[l31])
# ZIJ = .25 * (self._Z[l33] - self._Z[l13] + self._Z[l11] - self._Z[l31])
# self.__DX[l22] = scf * (alpha * XII + beta * XIJ + gamma * XJJ)
# self.__DY[l22] = scf * (alpha * YII + beta * YIJ + gamma * YJJ)
# self.__DZ[l22] = scf * (alpha * ZII + beta * ZIJ + gamma * ZJJ)
[docs]
def smooth_apply(self, face_def: Union[list, BlockFace, Face]) -> None:
"""
:param face_def:
"""
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
l0 = self.__get_lp(i0 - 1, j0 - 1, k0 - 1)
ldi = self._id[idir - 1] + self._ni * self._jd[idir - 1] + self._ni * self._nj * self._kd[idir - 1]
ldj = self._id[jdir - 1] + self._ni * self._jd[jdir - 1] + self._ni * self._nj * self._kd[jdir - 1]
ii, jj = np.meshgrid(np.arange(0, npi), np.arange(0, npj))
ll = l0 + ldi * ii + ldj * jj
self._X[ll] += self.__DX[ll]
self._Y[ll] += self.__DY[ll]
self._Z[ll] += self.__DZ[ll]
[docs]
def get_glue_pointers(self,
face_def: Union[list, BlockFace, Face],
face_type: Any = FaceType.A_SIDE) -> Tuple[ndarray, ...]:
"""
:param face_def:
:param face_type:
:return:
"""
i0, j0, k0, idir, jdir, npi, npj = self.__process_face_def(face_def)
l0 = self.__get_lp(i0 - 1, j0 - 1, k0 - 1)
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]
ldip = self._id[idir - 1] + self._ni * self._jd[idir - 1] + self._ni * self._nj * self._kd[idir - 1]
ldjp = self._id[jdir - 1] + self._ni * self._jd[jdir - 1] + self._ni * self._nj * self._kd[jdir - 1]
ldkp = idkd + self._ni * jdkd + self._ni * self._nj * kdkd
ll1 = []
ll2 = []
k = 1
if face_type == FaceType.A_SIDE:
k = -1
for j in range(npj):
for i in range(npi):
ll1.append(l0 + ldip * i + ldjp * j)
ll2.append(l0 + ldip * i + ldjp * j + k * ldkp)
return self._X[ll1], self._Y[ll1], self._Z[ll1], \
self._X[ll2], self._Y[ll2], self._Z[ll2], \
self.__DX[ll1], self.__DY[ll1], self.__DZ[ll1]