"""mesh placeholder - Grid class definition"""
import os
import sys
import re
import struct
from typing import Any, Union, Tuple, BinaryIO
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
from cupy import ndarray, float64
else:
import numpy as np
from numpy import ndarray, float64
import matplotlib.pyplot as plt
from g3d.defs import XMLDataItem
from g3d.defs import G3DWorld
from g3d.defs import G3DCode
from g3d.defs import G3DDisplay
from g3d.defs import G3DSizes
from g3d.defs import assure_path_exists
from g3d.defs import clean_path
from g3d.defs import DIMMode
from g3d.defs import BCTypes
from g3d.defs import INTFTypes
from g3d.defs import CodeType
from g3d.defs import BinaryFormat
from g3d.mesh.geo import Coordinates
from g3d.mesh.geo import Block
from g3d.mesh.geo import Face
from g3d.mesh.geo import Glue
[docs]
class Grid(G3DCode):
def __init__(self, *args: Any, **kwargs: Any) -> None:
self._blocks = []
self._nodes = Coordinates(tag='node_coordinates')
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.__bc = {}
self.__intf = {}
self.__bc_types = {}
self.__intf_types = {}
self.__smooth_faces = []
self.__smooth_face_glue = []
self.__default_width = 1.0
super().__init__(*args, **kwargs)
def _set_code_type(self) -> None:
g3d_world = G3DWorld()
g3d_world.info.code_type = CodeType.MESHER
g3d_world.info.code_name = 'g3dmesh.py'
g3d_world.info.code_description = 'a python-based mesh tool for generation of structured multi-block meshes'
g3d_world.options.short_opt += 'p'
g3d_world.options.long_opt.append('plot_edges')
def _print_help_message(self) -> None:
super()._print_help_message()
def _evaluate_run_option(self, opt: str, arg: Any) -> None:
g3d_world = G3DWorld()
super()._evaluate_run_option(opt, arg)
if opt in ['-p', '--plot_edges']:
g3d_world.display.plot_edges = True
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
g3d_world = G3DWorld()
if super().set_from_key_value(key, value):
return True
elif key == 'plot_edges':
g3d_world.display.plot_edges = value
return True
elif key == 'width':
self.__default_width = value
return False
[docs]
def allocate_arrays(self) -> None:
"""
"""
g3d_sizes = G3DSizes()
self._nodes.allocate_XYZ(g3d_sizes.n_nodes + g3d_sizes.n_ghost_nodes)
[docs]
def redefine_ghosts(self) -> None:
"""
"""
g3d_sizes = G3DSizes()
self._nodes.backup()
self._nodes.allocate_XYZ(g3d_sizes.n_nodes + g3d_sizes.n_ghost_nodes)
self._nodes.get_from_backup(True)
[docs]
def print_extreme_values(self) -> None:
"""
"""
g3d_world = G3DWorld()
g3d_world.display.header('computational domain')
self._nodes.print_extreme_values()
@property
def n_blocks(self) -> int:
"""
:return:
"""
return len(self._blocks)
@property
def blocks(self) -> list:
"""
:return:
"""
return self._blocks
@property
def nodes(self) -> Coordinates:
"""
:return:
"""
return self._nodes
@property
def n_nodes(self) -> int:
"""
:return:
"""
return G3DSizes().n_nodes
@property
def n_cells(self) -> int:
"""
:return:
"""
return G3DSizes().n_cells
[docs]
def nullify_boundary_lists(self) -> None:
"""
"""
self.__bc = {}
self.__intf = {}
self.__bc_types = {}
self.__intf_types = {}
[docs]
def nullify_smooth_faces(self) -> None:
"""
"""
self.__smooth_faces = []
self.__smooth_face_glue = []
[docs]
def print(self) -> None:
g3d_world = G3DWorld()
g3d_world.display.header('grid report')
print(' number of mesh blocks = %d' % len(self._blocks))
for b in self._blocks:
b.print()
print()
g3d_world.display.hline()
print(' node count = %d' % g3d_world.sizes.n_nodes)
print(' cell count = %d' % g3d_world.sizes.n_cells)
g3d_world.display.hline()
[docs]
def get_extreme_coordinates(self) -> Tuple[float64, ...]:
"""
:return:
"""
return self._nodes.get_extreme_values()
[docs]
def get_x_min_x_max(self) -> Tuple[float64, float64]:
"""
:return:
"""
return self._nodes.get_x_min_x_max()
[docs]
def get_y_min_y_max(self) -> Tuple[float64, float64]:
"""
:return:
"""
return self._nodes.get_y_min_y_max()
[docs]
def get_z_min_z_max(self) -> Tuple[float64, float64]:
"""
:return:
"""
return self._nodes.get_z_min_z_max()
[docs]
def scale(self, sf: float) -> None:
"""
:param sf:
"""
for b in self._blocks:
b.scale(sf)
[docs]
def add_smooth_face(self, face_def: Union[Face, list]) -> None:
"""
:param face_def:
"""
if type(face_def) is Face:
face_def = face_def.get_def()
self.__smooth_faces.append(face_def)
[docs]
def add_smooth_face_glue(self, faceA_def: Union[Face, list], faceB_def: Union[Face, list]) -> None:
"""
:param faceA_def:
:param faceB_def:
"""
if type(faceA_def) is Face:
faceA_def = faceA_def.get_def()
if type(faceB_def) is Face:
faceB_def = faceB_def.get_def()
g = Glue(faceA_def, faceB_def)
self.__smooth_face_glue.append(g)
[docs]
def check_volumes(self) -> int:
"""
:return:
"""
n_neg_global = 0
v_min_global = 1.e20
v_max_global = -1.e20
for b in self._blocks:
n_neg, v_min, v_max = b.check_volumes()
n_neg_global += n_neg
v_min_global = min([v_min_global, v_min])
v_max_global = max([v_max_global, v_max])
g3d_disp = G3DDisplay()
if not g3d_disp.silent:
g3d_disp.hline()
print('cell volume check completed')
if n_neg_global > 0:
print('%d negative volumes found' % n_neg_global)
else:
print('no negative volumes found')
print('cell volume min/max: %e/%e' % (v_min_global, v_max_global))
g3d_disp.hline()
return n_neg_global
[docs]
def fill(self) -> None:
"""
"""
g3d_world = G3DWorld()
dm = g3d_world.options.dim_mode
if dm == DIMMode.G2D or dm == DIMMode.G2D_AXISYMM or dm == DIMMode.Q3D:
for b in self._blocks:
b.fill_2D(self.__default_width)
elif dm == DIMMode.G3D:
for b in self._blocks:
b.fill()
[docs]
def export_connections(self, bcf: str = '') -> None:
"""
:param bcf:
"""
f = sys.stdout
if bcf != 'screen':
g3d_world = G3DWorld()
if len(bcf) == 0:
bcf = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.bc')
elif '/' not in bcf:
bcf = clean_path(g3d_world.io.wd + '/' + bcf)
assure_path_exists(bcf)
if not g3d_world.display.silent:
print('exporting face definitions to file: ', bcf)
f = open(bcf, 'w')
if len(self.__intf) > 0:
for tag in self.__intf.keys():
print(self.__intf_types[tag].name, file=f)
if self.__intf_types[tag].name.lower() != tag:
print(tag, file=f)
print(len(self.__intf[tag]), file=f)
for intf in self.__intf[tag]:
intfA = intf[0]
intfB = intf[1]
deltaA = [0., 0., 0.]
alphaA = 0.
deltaB = [0., 0., 0.]
alphaB = 0.
if type(intfA) is Face:
deltaA = intfA.get_delta()
alphaA = intfA.alpha
intfA = intfA.get_def()
if type(intfB) is Face:
deltaB = intfB.get_delta()
alphaB = intfB.alpha
intfB = intfB.get_def()
if self.__intf_types[tag] is INTFTypes.ROTPERINT:
print(*intfA, *alphaA, file=f)
print(*intfB, *alphaB, file=f)
elif self.__intf_types[tag] is INTFTypes.TRANSPERINT:
print(*intfA, *deltaA, file=f)
print(*intfB, *deltaB, file=f)
else:
print(*intfA, file=f)
print(*intfB, file=f)
print('', file=f)
if len(self.__bc) > 0:
for tag in self.__bc.keys():
print(self.__bc_types[tag].name, file=f)
print(tag, file=f)
print(len(self.__bc[tag]), file=f)
for bc in self.__bc[tag]:
if type(bc) is Face:
bc = bc.get_def()
print(*bc, file=f)
print('', file=f)
if bcf != 'screen':
f.close()
[docs]
def add_Intf(self, faceA_def: Union[Face, list[int]], faceB_def: Union[Face, list[int]], **kwargs: Any) -> None:
"""
:param faceA_def:
:param faceB_def:
:param kwargs:
"""
# default tag will just be used for grouping of faces
# i.e. will not be added in the text output
tag = 'default'
intf_type = INTFTypes.INTERFACES
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 == 'tag' or key == 'surface_tag':
tag = value
elif key == 'intf_type' or key == 'type':
intf_type = value
else:
print('warning! parameter %s is not defined' % key)
if tag == 'default':
tag = intf_type.name.lower()
if tag not in self.__intf.keys():
self.__intf[tag] = []
self.__intf_types[tag] = intf_type
self.__intf[tag].append([faceA_def, faceB_def])
[docs]
def add_BC(self, bc: Union[Face, list[int]], **kwargs: Any) -> None:
"""
:param bc:
:param kwargs:
"""
tag = 'walls'
bc_type = BCTypes.WALLS
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 == 'tag' or key == 'surface_tag':
tag = value
elif key == 'bc_type' or key == 'type':
bc_type = value
else:
print('warning! parameter %s is not defined' % key)
if tag not in self.__bc.keys():
self.__bc[tag] = []
self.__bc_types[tag] = bc_type
self.__bc[tag].append(bc)
[docs]
def add_IO(self, io: Union[Face, list[int]], **kwargs: Any) -> None:
"""
:param io:
:param kwargs:
"""
# this is just a wrapper for addBC
tag = 'ioflow'
bc_type = BCTypes.IOFLOW
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 == 'tag' or key == 'surface_tag':
tag = value
elif key == 'bc_type':
bc_type = value
else:
print('warning! parameter %s is not defined' % key)
self.add_BC(io, tag=tag, bc_type=bc_type)
[docs]
def add_Block(self, ni: Union[int, Block], nj: int = 0, nk: int = 0) -> None:
"""
:param ni:
:param nj:
:param nk:
"""
g3d_sizes = G3DSizes()
if type(ni) is Block:
b = ni
else:
b = Block(ni, nj, nk)
im = len(self._blocks)
b.im = im + 1
b.nmp = g3d_sizes.n_nodes
b.nmc = g3d_sizes.n_cells
self._blocks.append(b)
g3d_sizes.add_nodes(b.n_nodes)
g3d_sizes.add_cells(b.n_cells)
[docs]
def get_XYZ(self, im: int, lp: int = None) -> Union[Tuple[ndarray, ...], Tuple[float64, ...]]:
"""
:param im:
:param lp:
:return:
"""
return self._blocks[im - 1].get_XYZ(lp)
[docs]
def get_X(self, im: int, lp: int = None) -> Union[ndarray, float64]:
"""
:param im:
:param lp:
:return:
"""
return self._blocks[im - 1].get_X(lp)
[docs]
def get_Y(self, im: int, lp: int = None) -> Union[ndarray, float64]:
"""
:param im:
:param lp:
:return:
"""
return self._blocks[im - 1].get_Y(lp)
[docs]
def get_Z(self, im: int, lp: int = None) -> Union[ndarray, float64]:
"""
:param im:
:param lp:
:return:
"""
return self._blocks[im - 1].get_Z(lp)
[docs]
def set_XYZ(self, im: int, X: Union[float, ndarray], Y: Union[float, ndarray],
Z: Union[float, ndarray], lp=None) -> None:
"""
:param im:
:param X:
:param Y:
:param Z:
:param lp:
"""
self._blocks[im - 1].set_XYZ(X, Y, Z, lp)
[docs]
def set_X(self, im: int, X: Union[float, ndarray], lp=None) -> None:
"""
:param im:
:param X:
:param lp:
"""
self._blocks[im - 1].set_X(X, lp)
[docs]
def set_Y(self, im: int, Y: Union[float, ndarray], lp=None) -> None:
"""
:param im:
:param Y:
:param lp:
"""
self._blocks[im - 1].set_Y(Y, lp)
[docs]
def set_Z(self, im: int, Z: Union[float, ndarray], lp=None) -> None:
"""
:param im:
:param Z:
:param lp:
"""
self._blocks[im - 1].set_Z(Z, lp)
[docs]
def copy_XYZ(self, im_to: int, l_to: Union[int, list[int], ndarray],
im_from: int, l_from: Union[int, list[int], ndarray]) -> None:
"""
:param im_to:
:param l_to:
:param im_from:
:param l_from:
"""
self._blocks[im_to - 1].X[l_to] = self._blocks[im_from - 1].X[l_from]
self._blocks[im_to - 1].Y[l_to] = self._blocks[im_from - 1].Y[l_from]
self._blocks[im_to - 1].Z[l_to] = self._blocks[im_from - 1].Z[l_from]
[docs]
def copy_X(self, im_to: int, l_to: Union[int, list[int], ndarray],
im_from: int, l_from: Union[int, list[int], ndarray]) -> None:
"""
:param im_to:
:param l_to:
:param im_from:
:param l_from:
"""
self._blocks[im_to - 1].X[l_to] = self._blocks[im_from - 1].X[l_from]
[docs]
def copy_Y(self, im_to: int, l_to: Union[int, list[int], ndarray],
im_from: int, l_from: Union[int, list[int], ndarray]) -> None:
"""
:param im_to:
:param l_to:
:param im_from:
:param l_from:
"""
self._blocks[im_to - 1].Y[l_to] = self._blocks[im_from - 1].Y[l_from]
[docs]
def copy_Z(self, im_to: int, l_to: Union[int, list[int], ndarray],
im_from: int, l_from: Union[int, list[int], ndarray]) -> None:
"""
:param im_to:
:param l_to:
:param im_from:
:param l_from:
"""
self._blocks[im_to - 1].Z[l_to] = self._blocks[im_from - 1].Z[l_from]
[docs]
def populate_face_coordinate_placeholders(self, im: int, b: Face) -> None:
"""
:param im:
:param b:
"""
self._blocks[im - 1].populate_face_coordinate_placeholders(b)
[docs]
def cell_get_block(self, kc: int) -> int:
"""
:param kc:
:return:
"""
for im in range(len(self._blocks)):
if self._blocks[im].clow <= kc < self._blocks[im].chigh:
return im
return -1
[docs]
def node_get_block(self, kp: int) -> int:
"""
:param kp:
:return:
"""
for im in range(len(self._blocks)):
if self._blocks[im].nlow <= kp < self._blocks[im].nhigh:
return im
return -1
[docs]
def cell_get_ijk(self, kc: int) -> Tuple[int, int, int]:
"""
:param kc:
:return:
"""
for b in self._blocks:
if b.clow <= kc < b.chigh:
return b.cell_get_ijk(kc)
return -1, -1, -1
[docs]
def node_get_ijk(self, kp: int) -> Tuple[int, int, int]:
"""
:param kp:
:return:
"""
for b in self._blocks:
if b.nlow <= kp < b.nhigh:
return b.node_get_ijk(kp)
return -1, -1, -1
[docs]
def copy_Face(self, faceA_def: Union[Face, list[int]], faceB_def: Union[Face, list[int]]) -> None:
"""
:param faceA_def:
:param faceB_def:
"""
if type(faceA_def) is Face:
faceA_def = faceA_def.get_def()
if type(faceB_def) is Face:
faceB_def = faceB_def.get_def()
imA, i0A, j0A, k0A, idirA, jdirA, npiA, npjA = faceA_def
imB, i0B, j0B, k0B, idirB, jdirB, npiB, npjB = faceB_def
niA, njA, nkA = self._blocks[imA - 1].get_block_sizes()
niB, njB, nkB = self._blocks[imB - 1].get_block_sizes()
i = np.arange(npiB)
for j in range(npjB):
iA = (i0A - 1) + self._id[idirA - 1] * i + self._id[jdirA - 1] * j
jA = (j0A - 1) + self._jd[idirA - 1] * i + self._jd[jdirA - 1] * j
kA = (k0A - 1) + self._kd[idirA - 1] * i + self._kd[jdirA - 1] * j
lA = iA + niA * jA + niA * njA * kA
iB = (i0B - 1) + self._id[idirB - 1] * i + self._id[jdirB - 1] * j
jB = (j0B - 1) + self._jd[idirB - 1] * i + self._jd[jdirB - 1] * j
kB = (k0B - 1) + self._kd[idirB - 1] * i + self._kd[jdirB - 1] * j
lB = iB + niB * jB + niB * njB * kB
self.copy_XYZ(imB, lB, imA, lA)
[docs]
def glue_blocks(self) -> None:
"""
"""
for tag in self.__intf.keys():
intf_type = self.__intf_types[tag]
for intf in self.__intf[tag]:
intfA = intf[0]
intfB = intf[1]
if intf_type is INTFTypes.ROTPERINT:
# self.copyFace(intfA, intfB)
print('rotperint not implemented')
elif intf_type is INTFTypes.TRANSPERINT:
dx = intfA.delta.x - intfB.delta.x
dy = intfA.delta.y - intfB.delta.y
dz = intfA.delta.z - intfB.delta.z
self.copy_Face(intfA, intfB)
self._blocks[intfB.im - 1].shift_Face(intfB, x_shift=dx, y_shift=dy, z_shift=dz)
else:
self.copy_Face(intfA, intfB)
[docs]
def plot_edges_XY(self) -> None:
"""
"""
g3d_world = G3DWorld()
if g3d_world.display.plot_edges:
for b in self._blocks:
for e in b.edges:
e.scatter_XY(plt)
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.show()
[docs]
def plot_edges_YZ(self) -> None:
"""
"""
g3d_world = G3DWorld()
if g3d_world.display.plot_edges:
for b in self._blocks:
for e in b.edges:
e.scatter_YZ(plt)
plt.xlabel('Y')
plt.ylabel('Z')
plt.axis('equal')
plt.show()
[docs]
def plot_edges_XZ(self) -> None:
"""
"""
g3d_world = G3DWorld()
if g3d_world.display.plot_edges:
for b in self._blocks:
for e in b.edges:
e.scatter_XZ(plt)
plt.xlabel('X')
plt.ylabel('Z')
plt.axis('equal')
plt.show()
[docs]
def import_mesh(self, *args, **kwargs: Any) -> None:
"""
imports mesh from file
:param args: arguments (BinaryFormat, file base name (str)), optional
:param kwargs: keyword value pairs, optional
valid key values: 'format' binary_format (defaults to BinaryFormat.NATIVE),
'msh_file' -> filename, 'topo_file' -> filename
"""
g3d_world = G3DWorld()
binary_format = BinaryFormat.NATIVE
msh_file = ''
topo_file = ''
if len(args) > 0:
for arg in args:
if type(arg) is BinaryFormat:
binary_format = arg
elif type(arg) is str:
msh_file = arg + '.geo'
topo_file = arg + '.topo'
for key, value in kwargs.items():
if key == 'format':
binary_format = value
elif key == 'msh_file':
msh_file = value
elif key == 'topo_file':
topo_file = value
if len(topo_file) == 0:
if len(msh_file) == 0:
msh_file = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.geo')
topo_file = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.topo')
else:
topo_file = msh_file.replace('.geo', '.topo')
g3d_world.io.msh = msh_file
g3d_world.io.topo = topo_file
if binary_format == BinaryFormat.NATIVE:
self.__import_native(**kwargs)
elif binary_format == BinaryFormat.G3D_NATIVE:
self.__import_G3D_native(**kwargs)
def __import_native(self, **kwargs):
g3d_world = G3DWorld()
file_name = ''
for key, value in kwargs.items():
if key == 'filename' or key == 'name':
file_name = value
topo = g3d_world.io.topo
if len(file_name) > 0:
topo = file_name + '.topo'
# read mesh topology and allocate mesh blocks
with open(topo, 'rb') as f:
D = struct.unpack('i', f.read(4))[0]
if D == 2:
if g3d_world.options.dim_mode == DIMMode.G3D:
g3d_world.options.dim_mode = DIMMode.G2D
elif D == 1:
g3d_world.options.dim_mode = DIMMode.Q1D
if g3d_world.options.dim_mode == DIMMode.Q1D:
ni = struct.unpack('i', f.read(4))[0]
nj = 2
nk = 2
self.add_Block(ni, nj, nk)
else:
N = struct.unpack('i', f.read(4))[0]
for im in range(N):
ni = struct.unpack('i', f.read(4))[0]
nj = struct.unpack('i', f.read(4))[0]
nk = struct.unpack('i', f.read(4))[0]
self.add_Block(ni, nj, nk)
# allocated mesh placeholders
self.allocate_arrays()
# read mesh coordinates
geo = g3d_world.io.msh
if len(file_name) > 0:
geo = file_name + '.geo'
with open(geo, 'rb') as f:
self._import_coordinates(f)
def __import_G3D_native(self, **kwargs: Any) -> None:
g3d_world = G3DWorld()
file_name = ''
for key, value in kwargs.items():
if key == 'filename' or key == 'name':
file_name = value
if len(file_name) == 0:
file_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.msh')
elif '/' not in file_name:
file_name = clean_path(g3d_world.io.wd + '/' + file_name)
with open(file_name, 'rb') as f:
line = f.readline()
print(line.strip())
N = struct.unpack('i', f.read(4))[0]
# read mesh topology and allocate mesh blocks
# var = struct.unpack('i', f.read(4))[0] # reading leadin
f.read(4) # reading leadin
for im in range(N):
# var = struct.unpack('i', f.read(4))[0] # reading N
f.read(4) # reading N
ni = struct.unpack('i', f.read(4))[0]
nj = struct.unpack('i', f.read(4))[0]
nk = struct.unpack('i', f.read(4))[0]
self.add_Block(ni, nj, nk)
# allocated mesh placeholders
self.allocate_arrays()
# read mesh coordinates
# var = struct.unpack('i', f.read(4))[0] # reading leadin
f.read(4) # reading leadin
self._import_coordinates(f)
def __import_XML(self, **kwargs: Any) -> None:
file_name = ''
for key, value in kwargs.items():
if key == 'filename' or key == 'name':
file_name = value
msh_path = os.path.dirname(os.path.abspath(file_name))
DataCollections = []
regex_reference = re.compile(
r'DataItem Reference="/Xdmf/DataItem\[(?P<data_ref>\d+)'
)
regex_mesh_data = re.compile(
r'DataItem'
r'\s+Dimensions="(?P<NK>\d+)'
r'\s+(?P<NJ>\d+)'
r'\s+(?P<NI>\d+)'
r'\s+(?P<ND>\d+)\s*"'
r'\s+NumberType="(?P<type>\w+)"'
r'\s+Precision="(?P<prec>\d+)"'
r'\s+Format="(?P<format>\w+)"'
r'\s+Endian="(?P<endian>\w+)"'
r'\s+Seek="(?P<pos>\d+)">\s*'
)
regex_collection = re.compile(
r'DataItem'
r'\s+Name="(?P<name>[\w\s\-.]+)"'
r'\s+ItemType="Collection')
regex_data_item_begin = re.compile(
r'\s*<DataItem'
)
regex_data_item_end = re.compile(
r'\s*</DataItem'
)
regex_block = re.compile(
r'Grid'
r'\s+Name="Block\s+(?P<blockno>\d+)"'
r'\s*GridType="Uniform"'
)
regex_topology = re.compile(
r'Topology'
r'\s+TopologyType="3DSMesh"'
r'\s+Dimensions="(?P<NK>\d+)'
r'\s+(?P<NJ>\d+)'
r'\s+(?P<NI>\d+)\s*"'
)
regex_filename = re.compile(
r'(?P<filename>[\w_\-.]+)'
)
with open(file_name, 'r') as f:
line = f.readline()
while line:
match = regex_collection.search(line)
if match: # collection definition found
data_item_level = 1
while data_item_level:
line = f.readline()
match = regex_mesh_data.search(line)
if match:
XMLData = XMLDataItem()
XMLData.ni = int(match.group('NI'))
XMLData.nj = int(match.group('NJ'))
XMLData.nk = int(match.group('NK'))
XMLData.n_dim = int(match.group('ND'))
XMLData.NumberType = match.group('type')
XMLData.precision = int(match.group('prec'))
XMLData.format = match.group('format')
XMLData.endian = match.group('endian')
XMLData.pos = int(match.group('pos'))
while True:
line = f.readline()
match = regex_filename.search(line)
if match:
XMLData.file = msh_path + "/" + match.group('filename')
break
DataCollections.append(XMLData)
match = regex_data_item_begin.search(line)
if match:
data_item_level += 1
else:
match = regex_data_item_end.search(line)
if match:
data_item_level -= 1
match = regex_block.search(line)
if match: # block definition found
while True:
line = f.readline()
match = regex_topology.search(line)
if match:
idx = len(self._blocks)
self.add_Block(int(match.group('NI')), int(match.group('NJ')), int(match.group('NK')))
while True:
line = f.readline()
match = regex_reference.search(line)
if match:
self._blocks[idx].filename = DataCollections[
int(match.group('data_ref')) - 1].GetFile()
self._blocks[idx].pos = DataCollections[int(match.group('data_ref')) - 1].get_pos()
break
break
line = f.readline()
# allocated mesh placeholders
self.allocate_arrays()
# read grid coordinates from file
for b in self._blocks:
with open(b.filename, 'rb') as f:
f.seek(b.pos)
b.import_XYZ(f)
self._import_coordinates_post()
def _import_coordinates(self, f: BinaryIO) -> None:
g3d_world = G3DWorld()
if g3d_world.options.dim_mode == DIMMode.G3D or g3d_world.options.dim_mode == DIMMode.Q3D:
for b in self._blocks:
b.import_XYZ(f)
elif g3d_world.options.dim_mode == DIMMode.G2D or g3d_world.options.dim_mode == DIMMode.G2D_AXISYMM:
for b in self._blocks:
b.import_XY(f)
elif g3d_world.options.dim_mode == DIMMode.Q1D:
ni, nj, nk = self._blocks[0].get_block_sizes()
X, Y, Z = self._blocks[0].get_XYZ()
for i in range(ni):
xp = struct.unpack('d', f.read(8))[0]
A = struct.unpack('d', f.read(8))[0]
rp = np.sqrt(A)
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._import_coordinates_post()
def _import_coordinates_post(self) -> None:
pass
[docs]
def export(self, *args: Any, **kwargs: Any) -> None:
"""
exports mesh to file
:param args: arguments (BinaryFormat), optional
:param kwargs: keyword value pairs, optional
valid key values: 'format' binary_format (defaults to BinaryFormat.NATIVE)
"""
binary_format = BinaryFormat.NATIVE
if len(args) > 0:
if type(args[0]) is BinaryFormat:
binary_format = args[0]
for key, value in kwargs.items():
if key == 'format':
binary_format = value
if binary_format == BinaryFormat.NATIVE:
self.__export_native(**kwargs)
elif binary_format == BinaryFormat.G3D_NATIVE:
self.__export_G3D_native(**kwargs)
else:
print('the specified export format is not supported')
def __export_native(self, *args: Any, **kwargs: Any) -> None:
file_name = ''
if len(args) > 0:
if type(args[0]) is str:
file_name = args[0]
for key, value in kwargs.items():
if key == 'filename' or key == 'name':
file_name = value
bo = sys.byteorder
g3d_world = G3DWorld()
if len(file_name) == 0:
file_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname)
file_name = clean_path(file_name)
msh = file_name + '.geo'
topo = file_name + '.topo'
xmf = file_name + '_msh.xmf'
assure_path_exists(file_name)
if not g3d_world.display.silent:
print('exporting mesh topology file: ', topo)
f = open(topo, 'w+b')
dim = 3
if g3d_world.options.dim_mode.value < DIMMode.G3D.value:
if g3d_world.options.dim_mode.value < DIMMode.G2D.value:
dim = 1
else:
dim = 2
f.write(dim.to_bytes(4, 'little'))
nb = len(self._blocks)
f.write(nb.to_bytes(4, 'little'))
for b in self._blocks:
im, ni, nj, nk = b.get_def()
f.write(ni.to_bytes(4, 'little'))
f.write(nj.to_bytes(4, 'little'))
f.write(nk.to_bytes(4, 'little'))
f.close()
msh = clean_path(msh)
if not g3d_world.display.silent:
print('exporting mesh coordinate file: ', msh)
f = open(msh, 'w+b')
fpos = [0]
if g3d_world.options.dim_mode == DIMMode.G3D or g3d_world.options.dim_mode == DIMMode.Q3D:
for b in self._blocks:
b.export_XYZ(f)
fpos.append(f.tell())
elif g3d_world.options.dim_mode == DIMMode.G2D or g3d_world.options.dim_mode == DIMMode.G2D_AXISYMM:
for b in self._blocks:
b.export_XY(f)
fpos.append(f.tell())
f.close()
if not g3d_world.display.silent:
print('generating xml meta data: ', xmf)
msh = msh.replace(os.path.dirname(xmf), '').strip('/')
f = open(xmf, 'w')
f.write('<?xml version=\"1.0\" ?>\n')
f.write('<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n')
f.write('<Xdmf Version=\"2.0\">\n')
for b in self._blocks:
im, ni, nj, nk = b.get_def()
f.write(' <DataItem Name=\"Collection ' + str(im) + '\" ItemType=\"Collection\">\n')
if g3d_world.options.dim_mode == DIMMode.G3D or g3d_world.options.dim_mode == DIMMode.Q3D:
f.write(' <DataItem Dimensions=\"' + str(nk) + ' ' + str(nj) + ' ' + str(ni)
+ ' 3\" NumberType=\"Float\" Precision=\"8\" Format=\"Binary\" Endian=\"'
+ bo + '\" Seek=\"' + str(fpos[im - 1]) + '\">\n')
elif g3d_world.options.dim_mode == DIMMode.G2D or g3d_world.options.dim_mode == DIMMode.G2D_AXISYMM:
f.write(' <DataItem Dimensions=\"1 ' + str(nj) + ' ' + str(ni)
+ ' 2\" NumberType=\"Float\" Precision=\"8\" Format=\"Binary\" Endian=\"'
+ bo + '\" Seek=\"' + str(fpos[im - 1]) + '\">\n')
f.write(' ' + msh + '\n')
f.write(' </DataItem>\n')
f.write(' </DataItem>\n')
f.write(' <Domain>\n')
for b in self._blocks:
im, ni, nj, nk = b.get_def()
f.write(' <Grid Name=\"Block ' + str(im) + '\" GridType=\"Uniform\">\n')
if g3d_world.options.dim_mode == DIMMode.G3D or g3d_world.options.dim_mode == DIMMode.Q3D:
f.write(' <Topology TopologyType=\"3DSMesh\" Dimensions=\"'
+ str(nk) + ' ' + str(nj) + ' ' + str(ni) + '\"/>\n')
f.write(' <Geometry GeometryType=\"XYZ\">\n')
elif g3d_world.options.dim_mode == DIMMode.G2D or g3d_world.options.dim_mode == DIMMode.G2D_AXISYMM:
f.write(' <Topology TopologyType=\"2DSMesh\" Dimensions=\"1 '
+ str(nj) + ' ' + str(ni) + '\"/>\n')
f.write(' <Geometry GeometryType=\"XY\">\n')
f.write(' <DataItem Reference=\"/Xdmf/DataItem[' + str(im) + ']/DataItem[1]\" />\n')
f.write(' </Geometry>\n')
f.write(' </Grid>\n')
f.write(' </Domain>\n')
f.write('</Xdmf>\n')
f.close()
def __export_G3D_native(self, *args: Any, **kwargs: Any) -> None:
file_name = ''
if len(args) > 0:
if type(args[0]) is str:
file_name = args[0]
for key, value in kwargs.items():
if key == 'filename' or key == 'name':
file_name = value
g3d_world = G3DWorld()
if len(file_name) == 0:
file_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname + '.msh')
elif '/' not in file_name:
file_name = clean_path(g3d_world.io.wd + '/' + file_name)
file_name = clean_path(file_name)
assure_path_exists(file_name)
print('exporting mesh file: ', file_name)
f = open(file_name, 'w+b')
header = str.encode("g3dflow binary grid (3D)\n")
if g3d_world.options.dim_mode is DIMMode.G2D:
header = str.encode("g3dflow binary grid (2D)\n")
f.write(header)
nb = len(self._blocks)
f.write(nb.to_bytes(4, 'little'))
f.write((nb * 4 * 4).to_bytes(4, 'little')) # leadin (topology bytes)
for b in self._blocks:
im, ni, nj, nk = b.get_def()
f.write(b.nmp.to_bytes(4, 'little'))
f.write(ni.to_bytes(4, 'little'))
f.write(nj.to_bytes(4, 'little'))
f.write(nk.to_bytes(4, 'little'))
# self.calculateSizes()
if g3d_world.options.dim_mode is DIMMode.G3D:
f.write((g3d_world.sizes.n_nodes * 3 * 8).to_bytes(4, 'little')) # leadin (coordinate bytes)
for b in self._blocks:
b.export_XYZ(f)
elif g3d_world.options.dim_mode is DIMMode.G2D:
f.write((g3d_world.sizes.n_nodes * 8).to_bytes(4, 'little')) # leadin (coordinate bytes)
for b in self._blocks:
b.export_XY(f)
f.close()
[docs]
def smooth(self, N: int = 1) -> None:
"""
:param N:
"""
for b in self._blocks:
b.allocate_smooth_delta()
self.__smooth_defined_faces(N)
self.__smooth_all_volumes(N)
for b in self._blocks:
b.reset_smooth_delta()
def __smooth_all_volumes(self, N: int = 1) -> None:
pass
def __smooth_defined_faces(self, N: int = 1) -> None:
for istep in range(N):
for face_def in self.__smooth_faces:
self._blocks[face_def[0] - 1].smooth_Face(face_def)
self.__smooth_glue()
for face_def in self.__smooth_faces:
self._blocks[face_def[0] - 1].smooth_apply(face_def)
def __smooth_glue(self) -> None:
for glue in self.__smooth_face_glue:
faceA_def, faceB_def = glue.get_face_defs()
bA = self._blocks[faceA_def.im - 1]
bB = self._blocks[faceB_def.im - 1]
NIA, NJA, NKA = bA.get_block_sizes()
NIB, NJB, NKB = bB.get_block_sizes()
i0A, j0A, k0A, idirA, jdirA, npiA, npjA = bA.get_face_pointers(faceA_def)
i0B, j0B, k0B, idirB, jdirB, npiB, npjB = bB.get_face_pointers(faceB_def)
l0A = bA.get_index(i0A, j0A, k0A)
idkd = self._jd[idirA - 1] * self._kd[jdirA - 1] - self._kd[idirA - 1] * self._jd[jdirA - 1]
jdkd = self._kd[idirA - 1] * self._id[jdirA - 1] - self._id[idirA - 1] * self._kd[jdirA - 1]
kdkd = self._id[idirA - 1] * self._jd[jdirA - 1] - self._jd[idirA - 1] * self._id[jdirA - 1]
ldiA = self._id[idirA - 1] + NIA * self._jd[idirA - 1] + NIA * NJA * self._kd[idirA - 1]
ldjA = self._id[jdirA - 1] + NIA * self._jd[jdirA - 1] + NIA * NJA * self._kd[jdirA - 1]
ldkA = idkd + NIA * jdkd + NIA * NJA * kdkd
l0B = bB.get_index(i0B, j0B, k0B)
idkd = self._jd[idirB - 1] * self._kd[jdirB - 1] - self._kd[idirB - 1] * self._jd[jdirB - 1]
jdkd = self._kd[idirB - 1] * self._id[jdirB - 1] - self._id[idirB - 1] * self._kd[jdirB - 1]
kdkd = self._id[idirB - 1] * self._jd[jdirB - 1] - self._jd[idirB - 1] * self._id[jdirB - 1]
ldiB = self._id[idirB - 1] + NIB * self._jd[idirB - 1] + NIB * NJB * self._kd[idirB - 1]
ldjB = self._id[jdirB - 1] + NIB * self._jd[jdirB - 1] + NIB * NJB * self._kd[jdirB - 1]
ldkB = idkd + NIB * jdkd + NIB * NJB * kdkd
XA, YA, ZA = bA.get_XYZ()
XB, YB, ZB = bB.get_XYZ()
if glue.type is INTFTypes.TRANSPERINT:
XA += faceA_def.delta.x
YA += faceA_def.delta.y
ZA += faceA_def.delta.z
XB += faceB_def.delta.x
YB += faceB_def.delta.y
ZB += faceB_def.delta.z
DXA, DYA, DZA = bA.get_delta_arrays()
DXB, DYB, DZB = bB.get_delta_arrays()
i1 = 1
i2 = npiA - 1
ldA = ldiA
ldB = ldiB
if npiA == 1:
i1 = 0
i2 = npiA
ldA = ldjA
ldB = ldjB
j1 = 1
j2 = npjA - 1
if npjA == 1:
j1 = 0
j2 = npjA
ii, jj = np.meshgrid(np.arange(i1, i2), np.arange(j1, j2))
lA = l0A + ldiA * ii + ldjA * jj
lB = l0B + ldiB * ii + ldjB * jj
XI = 0.5 * (XB[lB + ldkB] - XA[lA - ldkA])
YI = 0.5 * (YB[lB + ldkB] - YA[lA - ldkA])
ZI = 0.5 * (ZB[lB + ldkB] - ZA[lA - ldkA])
XJ = 0.5 * (XA[lA + ldA] - XA[lA - ldA])
YJ = 0.5 * (YA[lA + ldA] - YA[lA - ldA])
ZJ = 0.5 * (ZA[lA + ldA] - ZA[lA - ldA])
XII = XB[lB + ldkB] - 2. * XA[lA] + XA[lA - ldkA]
YII = YB[lB + ldkB] - 2. * YA[lA] + YA[lA - ldkA]
ZII = ZB[lB + ldkB] - 2. * ZA[lA] + ZA[lA - ldkA]
XJJ = XA[lA + ldA] - 2. * XA[lA] + XA[lA - ldA]
YJJ = YA[lA + ldA] - 2. * YA[lA] + YA[lA - ldA]
ZJJ = ZA[lA + ldA] - 2. * ZA[lA] + ZA[lA - ldA]
XIJ = 0.25 * (XB[lB + ldB + ldkB] - XA[lA + ldA - ldkA] + XA[lA - ldA - ldkA] - XB[lB - ldB + ldkB])
YIJ = 0.25 * (YB[lB + ldB + ldkB] - YA[lA + ldA - ldkA] + YA[lA - ldA - ldkA] - YB[lB - ldB + ldkB])
ZIJ = 0.25 * (ZB[lB + ldB + ldkB] - ZA[lA + ldA - ldkA] + ZA[lA - ldA - ldkA] - ZB[lB - ldB + ldkB])
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.)
DXA[lA] = DXB[lB] = scf * (alpha * XII + beta * XIJ + gamma * XJJ)
DYA[lA] = DYB[lB] = scf * (alpha * YII + beta * YIJ + gamma * YJJ)
DZA[lA] = DZB[lB] = scf * (alpha * ZII + beta * ZIJ + gamma * ZJJ)
if glue.type is INTFTypes.TRANSPERINT:
XA -= faceA_def.delta.x
YA -= faceA_def.delta.y
ZA -= faceA_def.delta.z
XB -= faceB_def.delta.x
YB -= faceB_def.delta.y
ZB -= faceB_def.delta.z