"""hump mesh module (single block)"""
import os
from typing import Any
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from enum import Enum
from g3d.defs import G3DWorld
from g3d.defs import G3DDisplay
from g3d.defs import AngleType
from g3d.defs import Axis
from g3d.defs import NodeDistribution
from g3d.defs import Direction
from g3d.defs import BlockFace
from g3d.defs import BCTypes
from g3d.defs import DIMMode
from g3d.defs import FaceOrientation
from g3d.mesh.modules import MeshModule
from g3d.mesh.geo import Coord
from g3d.mesh.geo import Edge
from g3d.mesh.geo import Block
DomainType = Enum('DomainType', ['HILL_LEFT', 'HILL_CENTER', 'HILL_RIGHT', 'UNDEFINED'])
[docs]
class HumpModule(MeshModule):
def __init__(self, *args: Any, **kwargs: Any) -> None:
self._ni = 1 # number of nodes in I-direction
self._nj = 1 # number of nodes in J-direction
self._nk = 2 # number of nodes in K-direction (set to 2 for 2D-meshes)
self._x1 = 0.0 # left-end x-coordinate
self._x2 = 1.0 # right-end x-coordinate
self._y1 = 0.0 # lower y-coordinate
self._y2 = 1.0 # upper y-coordinate
self._z1 = 0.0 # lower z-coordinate
self._z2 = 1.0 # upper z-coordinate
self._dy_wall = 0.001 # height of first cell in wall normal direction
self._nj_wall_offset = 0 # j-index of wall offset edge
self._stretch_rate_x = 1.0 # cell stretch rate in the x-direction
self._stretch_rate_y = 2.0 # cell stretch rate in the y-direction
self._shape_control = 0.5 # shape factor for inlet and outlet faces
self._x_control = 0. # location of mesh control (x-coordinate)
self._dxf_control = 1. # spacing factor at control point
self._top_spacing_factor = 0.5 # spacing factor at top edge
self._n_control_lines = 2 # number of control lines
self.__ni_control = 1
self.__domain_type = DomainType.UNDEFINED
super().__init__(*args, **kwargs)
def _initialize_module(self, **kwargs: Any) -> None:
super()._initialize_module(**kwargs)
def _set_code_type(self) -> None:
super()._set_code_type()
g3d_world = G3DWorld()
g3d_world.info.code_sub_description = 'hill mesh module'
@property
def ni(self) -> int:
"""
:return:
"""
return self._ni
@ni.setter
def ni(self, val: int) -> None:
self._ni = val
@property
def nj(self) -> int:
"""
:return:
"""
return self._nj
@nj.setter
def nj(self, val: int) -> None:
self._nj = val
@property
def nk(self) -> int:
"""
:return:
"""
return self._nk
@nk.setter
def nk(self, val: int) -> None:
self._nk = val
@property
def x1(self) -> float:
"""
:return:
"""
return self._x1
@x1.setter
def x1(self, val: float) -> None:
self._x1 = val
@property
def x2(self) -> float:
"""
:return:
"""
return self._x2
@x2.setter
def x2(self, val: float) -> None:
self._x2 = val
@property
def y1(self) -> float:
"""
:return:
"""
return self._y1
@y1.setter
def y1(self, val: float) -> None:
self._y1 = val
@property
def y2(self) -> float:
"""
:return:
"""
return self._y2
@y2.setter
def y2(self, val: float) -> None:
self._y2 = val
@property
def z1(self) -> float:
"""
:return:
"""
return self._z1
@z1.setter
def z1(self, val: float) -> None:
self._z1 = val
@property
def z2(self) -> float:
"""
:return:
"""
return self._z2
@z2.setter
def z2(self, val: float) -> None:
self._z2 = val
@property
def dy_wall(self) -> float:
"""
:return:
"""
return self._dy_wall
@dy_wall.setter
def dy_wall(self, val: float) -> None:
self._dy_wall = val
@property
def x_control(self) -> float:
"""
:return:
"""
return self._x_control
@x_control.setter
def x_control(self, val: float) -> None:
self._x_control = val
@property
def dxf_control(self) -> float:
"""
:return:
"""
return self._dxf_control
@dxf_control.setter
def dxf_control(self, val: float) -> None:
self._dxf_control = val
@property
def nj_wall_offset(self) -> int:
"""
:return:
"""
return self._nj_wall_offset
@nj_wall_offset.setter
def nj_wall_offset(self, val: int) -> None:
self._nj_wall_offset = val
@property
def stretch_rate_x(self) -> float:
"""
:return:
"""
return self._stretch_rate_x
@stretch_rate_x.setter
def stretch_rate_x(self, val: float) -> None:
self._stretch_rate_x = val
@property
def stretch_rate_y(self) -> float:
"""
:return:
"""
return self._stretch_rate_y
@stretch_rate_y.setter
def stretch_rate_y(self, val: float) -> None:
self._stretch_rate_y = val
@property
def shape_control(self) -> float:
"""
:return:
"""
return self._shape_control
@shape_control.setter
def shape_control(self, val: float) -> None:
self._shape_control = val
@property
def top_spacing_factor(self) -> float:
"""
:return:
"""
return self._top_spacing_factor
@top_spacing_factor.setter
def top_spacing_factor(self, val: float) -> None:
self._top_spacing_factor = val
@property
def n_control_lines(self) -> int:
"""
:return:
"""
return self._n_control_lines
@n_control_lines.setter
def n_control_lines(self, val: int) -> None:
self._n_control_lines = val
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
g3d_world = G3DWorld()
if super().set_from_key_value(key, value):
if key == 'dimension' or key == 'dim':
if g3d_world.options.dim_mode == DIMMode.G2D:
self._nk = 2
self._z1 = 0.
self._z2 = 1.
return True
return False
def _generate_mesh(self) -> None:
# =====================================================================
# Allocate block
# =====================================================================
b1 = Block(self._ni, self._nj, self._nk)
self.add_Block(b1)
# =====================================================================
# Generate mesh-defining node coordinates
# =====================================================================
c11 = Coord(self._x1, self._data[0, 1], self._z1)
c21 = Coord(self._x2, self._data[-1, 1], self._z1)
c12 = Coord(self._x1, self._y2, self._z1)
c22 = Coord(self._x2, self._y2, self._z1)
b1.add_Coord(c11)
b1.add_Coord(c21)
b1.add_Coord(c12)
b1.add_Coord(c22)
# =====================================================================
# Generate block boundary edges
# =====================================================================
# define domain bottom edge
e1 = Edge(self._ni, c11, c21, distribute=NodeDistribution.CURVE_1,
direction=Direction.X_DIRECTION, data=self._data)
b1.add_Edge(e1, [1, 1, 1, 1, self._ni])
# define domain top edge
if self._dxf_control != 1.0:
self.__ni_control = e1.find_index(self._x_control, axis=Axis.X_AXIS) + 1
self._x_control = e1.get_X(self.__ni_control - 1)
c12m = Coord(c12, spacing=e1.h0)
c22m = Coord(c22, spacing=e1.h1)
c0 = Coord(self._x_control, self._y2, self._z1, spacing=e1.h0 * self._dxf_control)
b1.add_Coord(c0)
if np.abs((self._x2 - self._x_control) - (self._x_control - self._x1)) < (2 * e1.h1):
# one CURVE_2 in each direction
self.__domain_type = DomainType.HILL_CENTER
e21 = Edge(self.__ni_control, c12m, c0, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e21, [1, self._nj, 1, 1, self.__ni_control])
e22 = Edge(self._ni - self.__ni_control + 1, c0, c22m, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e22, [self.__ni_control, self._nj, 1, 1, self._ni - self.__ni_control + 1])
elif (self._x2 - self._x_control) > (self._x_control - self._x1):
# one CURVE_2 on left side and two CURVE_2 on right side
self.__domain_type = DomainType.HILL_LEFT
c1 = Coord(2. * self._x_control - self._x1, self._y2, self._z1, spacing=e1.h0)
b1.add_Coord(c1)
e21 = Edge(self.__ni_control, c12m, c0, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e21, [1, self._nj, 1, 1, self.__ni_control])
e22 = Edge(self.__ni_control, c0, c1, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e22, [self.__ni_control, self._nj, 1, 1, self.__ni_control])
e23 = Edge(self._ni - 2 * self.__ni_control + 2, c1, c22m, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e23, [2 * self.__ni_control - 1, self._nj, 1, 1,
self._ni - 2 * self.__ni_control + 2])
else:
# two CURVE_2 on left side and one CURVE_2 on right side
self.__domain_type = DomainType.HILL_RIGHT
c1 = Coord(2. * self._x_control - self._x2, self._y2, self._z1, spacing=e1.h1)
b1.add_Coord(c1)
self.__ni_control = self._ni - self.__ni_control + 1
e21 = Edge(self._ni - 2 * self.__ni_control + 2, c12m, c1, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e21, [1, self._nj, 1, 1, self._ni - 2 * self.__ni_control + 2])
e22 = Edge(self.__ni_control, c1, c0, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e22, [self._ni - 2 * self.__ni_control + 2, self._nj, 1, 1, self.__ni_control])
e23 = Edge(self.__ni_control, c0, c22m, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e23, [self._ni - self.__ni_control + 1, self._nj, 1, 1, self.__ni_control])
else:
# no mesh control defined (one CURVE_2)
e2 = Edge(self._ni, c12, c22, distribute=NodeDistribution.CURVE_2)
b1.add_Edge(e2, [1, self._nj, 1, 1, self._ni])
# define left and right domain edge (CURVE_3)
c11m = Coord(c11, ialt=AngleType.IALT_2, alpha=0., beta=0.,
shape_factor=self._shape_control, spacing=self._dy_wall)
e3 = Edge(self._nj, c11m, c12, distribute=NodeDistribution.CURVE_3,
stretch_rate=self._stretch_rate_y)
b1.add_Edge(e3, [1, 1, 1, 3, self._nj])
c21m = Coord(c21, ialt=AngleType.IALT_2, alpha=0., beta=0.,
shape_factor=self._shape_control, spacing=self._dy_wall)
e4 = Edge(self._nj, c21m, c22, distribute=NodeDistribution.CURVE_3,
stretch_rate=self._stretch_rate_y)
b1.add_Edge(e4, [self._ni, 1, 1, 3, self._nj])
self._fill_klow_plane(e1)
def _fill_klow_plane(self, e1):
b1 = self._blocks[0]
b1.fill_Face(BlockFace.KLOW)
if self._nj_wall_offset > 1:
nj2 = self._nj - self._nj_wall_offset + 1
# =====================================================================
# replace inlet and outlet edges
# =====================================================================
e0 = b1.get_Edge(self.__ni_control, self._nj_wall_offset, 1, 3, nj2)
spacing_top = e0.h1 * self._top_spacing_factor
# =====================================================================
# Generate offset edge with a node distribution where each node is
# located in the wall normal direction
# =====================================================================
c11 = e1.start
c12 = b1.get_Coord_from_IJK(1, self._nj_wall_offset, 1)
c22 = b1.get_Coord_from_IJK(self._ni, self._nj_wall_offset, 1)
d = c11.get_dist(c12)
e2 = e1.generate_offset_Edge(d)
e2.set_start_end_Coord(c12, c22)
b1.add_Edge(e2, [1, self._nj_wall_offset, 1, 1, self._ni])
# =====================================================================
# Fill 2D plane with nodes
# =====================================================================
b1.fill_Face([1, 1, 1, 1, 3, self._ni, self._nj_wall_offset])
b1.fill_Face([1, self._nj_wall_offset, 1, 1, 3, self._ni, nj2])
if self.__domain_type != DomainType.UNDEFINED and self._n_control_lines > 0:
# re-define domain left edge (CURV4 from offset out)
e1 = b1.get_Edge(1, 1, 1, 3, self._nj_wall_offset)
e2 = b1.get_Edge(1, self._nj_wall_offset, 1, 3, nj2)
shape_factor = 0.5 * (e2.end.y - e2.start.y)
c1 = Coord(e2.start,
ialt=AngleType.IALT_1,
alpha=90.,
beta=0.,
shape_factor=shape_factor,
spacing=e1.h1,
angle_dim=AngleType.DEGREES)
c2 = Coord(e2.end,
ialt=AngleType.IALT_1,
alpha=90.,
beta=0.,
shape_factor=shape_factor,
spacing=spacing_top,
angle_dim=AngleType.DEGREES)
e3 = Edge(nj2, c1, c2, distribute=NodeDistribution.CURVE_4, stretch_rate=2.0)
b1.add_Edge(e3, [1, self._nj_wall_offset, 1, 3, nj2])
# re-define domain right edge (CURV4 from offset out)
e1 = b1.get_Edge(self._ni, 1, 1, 3, self._nj_wall_offset)
e2 = b1.get_Edge(self._ni, self._nj_wall_offset, 1, 3, nj2)
shape_factor = 0.5 * (e2.end.y - e2.start.y)
c1 = Coord(e2.start,
ialt=AngleType.IALT_1,
alpha=90.,
beta=0.,
shape_factor=shape_factor,
spacing=e1.h1,
angle_dim=AngleType.DEGREES)
c2 = Coord(e2.end,
ialt=AngleType.IALT_1,
alpha=90.,
beta=0.,
shape_factor=shape_factor,
spacing=spacing_top,
angle_dim=AngleType.DEGREES)
e3 = Edge(nj2, c1, c2, distribute=NodeDistribution.CURVE_4, stretch_rate=2.0)
b1.add_Edge(e3, [self._ni, self._nj_wall_offset, 1, 3, nj2])
# add help curves and refill KLOW plane
nid = int(self.__ni_control / (self._n_control_lines + 1))
ic0 = (self.__ni_control - 1) - self._n_control_lines * (nid - 1) + 1
if self.__domain_type == DomainType.HILL_RIGHT:
ic0 = (self._ni - 1) - (self.__ni_control - 1) - self._n_control_lines * (nid - 1) + 1
i_fill = 1
for i in range(2 * self._n_control_lines + 1):
ic = ic0 + i * (nid - 1)
e1 = b1.get_Edge(ic, 1, 1, 3, self._nj_wall_offset)
e2 = b1.get_Edge(ic, self._nj_wall_offset, 1, 3, nj2)
shape_factor = 0.5 * (e2.end.y - e2.start.y)
alpha, beta = e1.get_angle(ialt=AngleType.IALT_1,
direction=Direction.START_TO_END,
angle_dim=AngleType.DEGREES)
c1 = Coord(e2.start,
ialt=AngleType.IALT_1,
alpha=alpha,
beta=beta,
shape_factor=shape_factor,
spacing=e1.h1,
angle_dim=AngleType.DEGREES)
c2 = Coord(e2.end,
ialt=AngleType.IALT_1,
alpha=90.,
beta=0.,
shape_factor=shape_factor,
spacing=spacing_top,
angle_dim=AngleType.DEGREES)
b1.add_Coord(c1)
b1.add_Coord(c2)
e3 = Edge(nj2, c1, c2, distribute=NodeDistribution.CURVE_4, stretch_rate=2.0)
b1.add_Edge(e3, [ic, self._nj_wall_offset, 1, 3, nj2])
b1.fill_Face([i_fill, self._nj_wall_offset, 1, 1, 3, ic - i_fill + 1, nj2])
i_fill = ic
b1.fill_Face([i_fill, self._nj_wall_offset, 1, 1, 3, self._ni - i_fill + 1, nj2])
def _generate_mesh_post(self) -> int:
self._smooth()
self._extrude_mesh()
return super()._generate_mesh_post()
def _smooth(self) -> None:
# =====================================================================
# Smooth 2D plane
# =====================================================================
if self._smooth_faces:
g3d_disp = G3DDisplay()
g3d_disp.header('smoothing mesh (surface smoother)')
b1 = self._blocks[0]
self.add_smooth_face(b1.get_Face(BlockFace.KLOW))
self.smooth(self._n_smooth)
def _extrude_mesh(self) -> None:
# =====================================================================
# Extrude 2D face in K direction
# =====================================================================
g3d_world = G3DWorld()
g3d_world.display.header('extruding 2D mesh plane in K-direction')
if g3d_world.options.dim_mode == DIMMode.G2D:
self._nk = 2
if self._nk == 2:
self._z2 = 1.
b1 = self._blocks[0]
c11 = b1.get_Coord_from_IJK(1, 1, 1)
c_ext = Coord(c11.x, c11.y, self._z2)
e_ext = Edge(self._nk, c11, c_ext, distribute=NodeDistribution.CURVE_2)
b1.extrude_Face(b1.get_Face(BlockFace.KLOW, FaceOrientation.INWARD_FACING_NORMAL), e_ext)
def _add_boundaries(self) -> None:
# =====================================================================
# Populate boundary list placeholders
# =====================================================================
b1 = self._blocks[0]
if self._nk > 2:
self.add_BC(b1.get_Face(BlockFace.KLOW), bc_type=BCTypes.SYMM, tag='symm')
self.add_BC(b1.get_Face(BlockFace.KHIGH), bc_type=BCTypes.SYMM, tag='symm')
self.add_IO(b1.get_Face(BlockFace.ILOW), tag='inlet')
self.add_IO(b1.get_Face(BlockFace.IHIGH), tag='outlet')
self.add_BC(b1.get_Face(BlockFace.JHIGH), bc_type=BCTypes.SYMM, tag='symm')
self.add_BC(b1.get_Face(BlockFace.JLOW), tag='walls')