import os
import copy
from datetime import date
import locale
import itertools
import numpy as np
from scipy import interpolate
from PySide2 import QtGui, QtCore, QtWidgets
import PyAero
import GraphicsItemsCollection as gic
import GraphicsItem
import Connect
from Utils import Utils
from Settings import OUTPUTDATA
import logging
logger = logging.getLogger(__name__)
class Windtunnel:
"""docstring for Windtunnel"""
def __init__(self):
# contains list of BlockMesh objects
self.blocks = []
# get MainWindow instance (overcomes handling parents)
self.mainwindow = QtCore.QCoreApplication.instance().mainwindow
def AirfoilMesh(self, name='', contour=None, divisions=15, ratio=3.0,
thickness=0.04):
# get airfoil contour coordinates
x, y = contour
# make a list of point tuples
# [(x1, y1), (x2, y2), (x3, y3), ... , (xn, yn)]
line = list(zip(x, y))
# block mesh around airfoil contour
self.block_airfoil = BlockMesh(name=name)
self.block_airfoil.addLine(line)
# self.block_airfoil.extrudeLine(line, length=thickness, direction=3,
# divisions=divisions, ratio=ratio)
self.block_airfoil.extrudeLine_cell_thickness(line,
cell_thickness=thickness,
growth=ratio,
divisions=divisions,
direction=3)
self.blocks.append(self.block_airfoil)
def TrailingEdgeMesh(self, name='', te_divisions=3,
thickness=0.04, divisions=10, ratio=1.05):
# compile first line of trailing edge block
first = self.block_airfoil.getLine(number=0, direction='v')
last = self.block_airfoil.getLine(number=-1, direction='v')
last_reversed = copy.deepcopy(last)
last_reversed.reverse()
vec = np.array(first[0]) - np.array(last[0])
line = copy.deepcopy(last_reversed)
# in case of TE add the points from the TE
if self.mainwindow.airfoil.has_TE:
for i in range(1, te_divisions):
p = last_reversed[-1] + float(i) / te_divisions * vec
# p is type numpy.float, so convert it to float
line.append((float(p[0]), float(p[1])))
line += first
# handle case with sharp trailing edge
else:
line += first[1:]
# trailing edge block mesh
block_te = BlockMesh(name=name)
block_te.addLine(line)
# block_te.extrudeLine(line, length=length, direction=4,
# divisions=divisions, ratio=ratio)
block_te.extrudeLine_cell_thickness(line,
cell_thickness=thickness,
growth=ratio,
divisions=divisions,
direction=4)
# equidistant point distribution
block_te.distribute(direction='u', number=-1)
# make a transfinite interpolation
# i.e. recreate pooints inside the block
block_te.transfinite()
self.block_te = block_te
self.blocks.append(block_te)
def TunnelMesh(self, name='', tunnel_height=2.0, divisions_height=100,
ratio_height=10.0, dist='symmetric'):
block_tunnel = BlockMesh(name=name)
self.tunnel_height = tunnel_height
# line composed of trailing edge and airfoil meshes
line = self.block_te.getVLines()[-1]
line.reverse()
del line[-1]
line += self.block_airfoil.getULines()[-1]
del line[-1]
line += self.block_te.getVLines()[0]
block_tunnel.addLine(line)
# line composed of upper, lower and front line segments
p1 = np.array((block_tunnel.getULines()[0][0][0], tunnel_height))
p2 = np.array((0.0, tunnel_height))
p3 = np.array((0.0, -tunnel_height))
p4 = np.array((block_tunnel.getULines()[0][-1][0], -tunnel_height))
# upper line of wind tunnel
line = list()
vec = p2 - p1
for t in np.linspace(0.0, 1.0, 10):
p = p1 + t * vec
line.append(p.tolist())
del line[-1]
# front half circle of wind tunnel
for phi in np.linspace(90.0, 270.0, 200):
phir = np.radians(phi)
x = tunnel_height * np.cos(phir)
y = tunnel_height * np.sin(phir)
line.append((x, y))
del line[-1]
# lower line of wind tunnel
vec = p4 - p3
for t in np.linspace(0.0, 1.0, 10):
p = p3 + t * vec
line.append(p.tolist())
line = np.array(line)
tck, u = interpolate.splprep(line.T, s=0, k=1)
# point distribution on upper, front and lower part
if dist == 'symmetric':
ld = -1.3
ud = 1.3
if dist == 'lower':
ld = -1.2
ud = 1.5
if dist == 'upper':
ld = -1.5
ud = 1.2
xx = np.linspace(ld, ud, len(block_tunnel.getULines()[0]))
t = (np.tanh(xx) + 1.0) / 2.0
line = interpolate.splev(t, tck, der=0)
line = list(zip(line[0].tolist(), line[1].tolist()))
block_tunnel.addLine(line)
p5 = np.array(block_tunnel.getULines()[0][0])
p6 = np.array(block_tunnel.getULines()[0][-1])
# first vline
vline1 = BlockMesh.makeLine(p5, p1, divisions=divisions_height,
ratio=ratio_height)
# last vline
vline2 = BlockMesh.makeLine(p6, p4, divisions=divisions_height,
ratio=ratio_height)
boundary = [block_tunnel.getULines()[0],
block_tunnel.getULines()[-1],
vline1,
vline2]
block_tunnel.transfinite(boundary=boundary)
# blending between normals (inner lines) and transfinite (outer lines)
ulines = list()
old_ulines = block_tunnel.getULines()
for j, uline in enumerate(block_tunnel.getULines()):
# skip first and last line
if j == 0 or j == len(block_tunnel.getULines()) - 1:
ulines.append(uline)
continue
line = list()
xo, yo = list(zip(*old_ulines[0]))
xo = np.array(xo)
yo = np.array(yo)
normals = BlockMesh.curveNormals(xo, yo)
for i, point in enumerate(uline):
# skip first and last point
if i == 0 or i == len(uline) - 1:
line.append(point)
continue
pt = np.array(old_ulines[j][i])
pto = np.array(old_ulines[0][i])
vec = pt - pto
# projection of vec into normal
dist = np.dot(vec, normals[i]) / np.linalg.norm(normals[i])
pn = pto + dist * normals[i]
v = float(j) / float(len(block_tunnel.getULines()))
exp = 0.6
pnew = (1.0 - v**exp) * pn + v**exp * pt
line.append((pnew.tolist()[0], pnew.tolist()[1]))
ulines.append(line)
block_tunnel = BlockMesh(name=name)
for uline in ulines:
block_tunnel.addLine(uline)
ij = [0, 30, 0, len(block_tunnel.getULines()) - 1]
block_tunnel.transfinite(ij=ij)
ij = [len(block_tunnel.getVLines()) - 31,
len(block_tunnel.getVLines()) - 1,
0,
len(block_tunnel.getULines()) - 1]
block_tunnel.transfinite(ij=ij)
sm = 1
if sm == 1:
smooth = Smooth(block_tunnel)
nodes = smooth.selectNodes(domain='interior')
block_tunnel = smooth.smooth(nodes, iterations=1,
algorithm='laplace')
ij = [1, 30, 1, len(block_tunnel.getULines()) - 2]
nodes = smooth.selectNodes(domain='ij', ij=ij)
block_tunnel = smooth.smooth(nodes, iterations=2,
algorithm='laplace')
ij = [len(block_tunnel.getVLines()) - 31,
len(block_tunnel.getVLines()) - 2,
1,
len(block_tunnel.getULines()) - 2]
nodes = smooth.selectNodes(domain='ij', ij=ij)
block_tunnel = smooth.smooth(nodes, iterations=3,
algorithm='laplace')
self.block_tunnel = block_tunnel
self.blocks.append(block_tunnel)
def TunnelMeshWake(self, name='', tunnel_wake=2.0,
divisions=100, ratio=0.1, spread=0.4):
chord = 1.0
block_tunnel_wake = BlockMesh(name=name)
# line composed of trailing edge and block_tunnel meshes
line = self.block_tunnel.getVLines()[-1]
line.reverse()
del line[-1]
line += self.block_te.getULines()[-1]
del line[-1]
line += self.block_tunnel.getVLines()[0]
block_tunnel_wake.addLine(line)
#
p1 = np.array((self.block_te.getULines()[-1][0][0],
self.tunnel_height))
p4 = np.array((self.block_te.getULines()[-1][-1][0],
- self.tunnel_height))
p7 = np.array((tunnel_wake + chord, self.tunnel_height))
p8 = np.array((tunnel_wake + chord, -self.tunnel_height))
upper = BlockMesh.makeLine(p7, p1, divisions=divisions,
ratio=1.0 / ratio)
lower = BlockMesh.makeLine(p8, p4, divisions=divisions,
ratio=1.0 / ratio)
left = line
right = BlockMesh.makeLine(p8, p7, divisions=len(left) - 1, ratio=1.0)
boundary = [upper, lower, right, left]
block_tunnel_wake.transfinite(boundary=boundary)
# equalize division line in wake
for i, u in enumerate(block_tunnel_wake.getULines()[0]):
if u[0] < chord + tunnel_wake * spread:
ll = len(block_tunnel_wake.getULines()[0])
line_no = -ll + i
break
block_tunnel_wake.distribute(direction='v', number=line_no)
# transfinite left of division line
ij = [len(block_tunnel_wake.getVLines()) + line_no,
len(block_tunnel_wake.getVLines()) - 1,
0,
len(block_tunnel_wake.getULines()) - 1]
block_tunnel_wake.transfinite(ij=ij)
# transfinite right of division line
ij = [0,
len(block_tunnel_wake.getVLines()) + line_no,
0,
len(block_tunnel_wake.getULines()) - 1]
block_tunnel_wake.transfinite(ij=ij)
self.block_tunnel_wake = block_tunnel_wake
self.blocks.append(block_tunnel_wake)
def makeMesh(self):
toolbox = self.mainwindow.centralwidget.toolbox
if self.mainwindow.airfoil:
if not hasattr(self.mainwindow.airfoil, 'spline_data'):
message = 'Splining needs to be done first.'
self.mainwindow.slots.messageBox(message)
return
contour = self.mainwindow.airfoil.spline_data[0]
else:
self.mainwindow.slots.messageBox('No airfoil loaded.')
return
# delete blocks outline if existing
# because a new one will be generated
if hasattr(self.mainwindow.airfoil, 'mesh_blocks'):
self.mainwindow.scene.removeItem(
self.mainwindow.airfoil.mesh_blocks)
del self.mainwindow.airfoil.mesh_blocks
progdialog = QtWidgets.QProgressDialog(
"", "Cancel", 0, 100, self.mainwindow)
progdialog.setFixedWidth(300)
progdialog.setMinimumDuration(0)
progdialog.setWindowTitle('Generating the CFD mesh')
progdialog.setWindowModality(QtCore.Qt.WindowModal)
progdialog.show()
progdialog.setValue(10)
# progdialog.setLabelText('making blocks')
self.AirfoilMesh(name='block_airfoil',
contour=contour,
divisions=toolbox.points_n.value(),
ratio=toolbox.ratio.value(),
thickness=toolbox.normal_thickness.value())
progdialog.setValue(20)
if progdialog.wasCanceled():
return
self.TrailingEdgeMesh(name='block_TE',
te_divisions=toolbox.te_div.value(),
thickness=toolbox.length_te.value(),
divisions=toolbox.points_te.value(),
ratio=toolbox.ratio_te.value())
progdialog.setValue(30)
if progdialog.wasCanceled():
return
self.TunnelMesh(name='block_tunnel',
tunnel_height=toolbox.tunnel_height.value(),
divisions_height=toolbox.divisions_height.value(),
ratio_height=toolbox.ratio_height.value(),
dist=toolbox.dist.currentText())
progdialog.setValue(50)
if progdialog.wasCanceled():
return
self.TunnelMeshWake(name='block_tunnel_wake',
tunnel_wake=toolbox.tunnel_wake.value(),
divisions=toolbox.divisions_wake.value(),
ratio=toolbox.ratio_wake.value(),
spread=toolbox.spread.value() / 100.0)
progdialog.setValue(70)
if progdialog.wasCanceled():
return
# connect mesh blocks
connect = Connect.Connect(progdialog)
vertices, connectivity, progdialog = \
connect.connectAllBlocks(self.blocks)
# add mesh to Windtunnel instance
self.mesh = vertices, connectivity
# generate cell to vertex connectivity from mesh
self.makeLCV()
# generate cell to edge connectivity from mesh
self.makeLCE()
# generate boundaries from mesh connectivity
unique, seen, doubles, boundary_edges = self.makeBoundaries()
# find loops inside boundary_edges
self.boundary_loops = self.findLoops(boundary_edges)
logger.info('Mesh around {} created'.
format(self.mainwindow.airfoil.name))
logger.info('Mesh has {} vertices and {} elements'.
format(len(vertices), len(connectivity)))
self.drawMesh(self.mainwindow.airfoil)
self.drawBlockOutline(self.mainwindow.airfoil)
progdialog.setValue(100)
# enable mesh export and set filename and boundary definitions
toolbox.box_meshexport.setEnabled(True)
def makeLCV(self):
"""Make cell to vertex connectivity for the mesh
LCV is identical to connectivity
"""
_, connectivity = self.mesh
self.LCV = connectivity
def makeLCE(self):
"""Make cell to edge connectivity for the mesh"""
_, connectivity = self.mesh
self.LCE = dict()
self.edges = list()
for i, cell in enumerate(connectivity):
# example for Qudrilateral:
# cell: [0, 1, 5, 4]
# local_edges: [(0,1), (1,5), (5,4), (4,0)]
local_edges = [(cell[j], cell[(j + 1) % len(cell)])
for j in range(len(cell))]
# all edges for cell i
self.LCE[i] = local_edges
# all edges in one list
self.edges += [tuple(sorted(edge)) for edge in local_edges]
def makeLCC(self):
"""Make cell to cell connectivity for the mesh"""
pass
def makeBoundaries(self):
"""A boundary edge is an edge that belongs only to one cell"""
seen = set()
unique = list()
doubles = set()
for edge in self.edges:
if edge not in seen:
seen.add(edge)
unique.append(edge)
else:
doubles.add(edge)
boundary_edges = [edge for edge in unique if edge not in doubles]
return unique, seen, doubles, boundary_edges
def findLoops(self, edges):
"""Find loops in a list of edges which are stored in tuples
and return the "connected components", in the disjoint set.
In the case of boundary edges these are loops or
"cycles" in graph theory language
Args:
edges (list of tuples):
Returns:
TYPE: Description
"""
vertices, connectivity = self.mesh
# make disjoint set object
djs = DisjointSet()
# add all edges to the disjoint set
for edge in edges:
djs.add(edge[0], edge[1])
# get the boundary loops (airfoil, outer boundary)
# djs.group returns a dictionary containing all loops
# the key is an arbitrary node of the loop
# the values per key are a list of unordered nodes
# belonging to the loop
boundary_loops = djs.group
# in order to order the returned nodes again, their corresponding edges
# have to be found first
new_loops = dict()
for i, loop in enumerate(boundary_loops):
l_edges = list()
for node in boundary_loops[loop]:
l_edges += [sorted(edge) for edge in edges if node in edge]
# remove duplicate list elements from l_edges
loop_edges = [k for k, _ in itertools.groupby(sorted(l_edges))]
new_loops[i] = loop_edges
# new_loops[0] is airfoil
# new_loops[1] is complete outer boundary
# split outer boundary into inlet and outlet
self.is_outlet = list()
# edge is e.g.: (27, 53)
for i, edge in enumerate(new_loops[1]):
self.is_outlet.append(0)
# vertices[edge[0]] is e.g.: (1.0453006577029285, 3.5)
vector = Utils.vector(vertices[edge[0]], vertices[edge[1]])
# check angle against y-axis
angle = Utils.angle_between(vector, (0., 1.), degree=True)
# FIXME
# FIXME find better criterions or at leat refactor
# FIXME
# angle tolerance
tol = 0.5
# check only for edges downstream the airfoil
tol_wake = 1.1
if ((angle > - tol) and (angle < tol)) or \
((angle > 180. - tol) and (angle < 180. + tol)) and \
vertices[edge[0]][0] > tol_wake:
self.is_outlet[i] = 1
return new_loops
def drawMesh(self, airfoil):
"""Add the mesh as ItemGroup to the scene
Args:
airfoil (TYPE): object containing all airfoil properties and data
"""
# toggle spline points
self.mainwindow.centralwidget.cb3.click()
# delete old mesh if existing
if hasattr(airfoil, 'mesh'):
logger.debug('MESH item type: {}'.format(type(airfoil.mesh)))
self.mainwindow.scene.removeItem(airfoil.mesh)
mesh = list()
for block in self.blocks:
for lines in [block.getULines(),
block.getVLines()]:
for line in lines:
# instantiate a graphics item
contour = gic.GraphicsCollection()
# make it polygon type and populate its points
points = [QtCore.QPointF(x, y) for x, y in line]
contour.Polyline(QtGui.QPolygonF(points), '')
# set its properties
contour.pen.setColor(QtGui.QColor(0, 0, 0, 255))
contour.pen.setWidthF(0.8)
contour.pen.setCosmetic(True)
contour.brush.setStyle(QtCore.Qt.NoBrush)
# add contour as a GraphicsItem to the scene
# these are the objects which are drawn in the GraphicsView
meshline = GraphicsItem.GraphicsItem(contour)
mesh.append(meshline)
airfoil.mesh = self.mainwindow.scene.createItemGroup(mesh)
# activate viewing options if mesh is created and displayed
self.mainwindow.centralwidget.cb6.setChecked(True)
self.mainwindow.centralwidget.cb6.setEnabled(True)
def drawBlockOutline(self, airfoil):
"""Add the mesh block outlines to the scene
Args:
airfoil (TYPE): object containing all airfoil properties and data
"""
# FIXME
# FIXME Refactroing of code duplication here and in drawMesh
# FIXME
mesh_blocks = list()
for block in self.blocks:
for lines in [block.getULines()]:
for line in [lines[0], lines[-1]]:
# instantiate a graphics item
contour = gic.GraphicsCollection()
# make it polygon type and populate its points
points = [QtCore.QPointF(x, y) for x, y in line]
contour.Polyline(QtGui.QPolygonF(points), '')
# set its properties
contour.pen.setColor(QtGui.QColor(202, 31, 123, 255))
contour.pen.setWidthF(3.0)
contour.pen.setCosmetic(True)
contour.brush.setStyle(QtCore.Qt.NoBrush)
# add contour as a GraphicsItem to the scene
# these are the objects which are drawn in the GraphicsView
meshline = GraphicsItem.GraphicsItem(contour)
mesh_blocks.append(meshline)
for lines in [block.getVLines()]:
for line in [lines[0], lines[-1]]:
# instantiate a graphics item
contour = gic.GraphicsCollection()
# make it polygon type and populate its points
points = [QtCore.QPointF(x, y) for x, y in line]
contour.Polyline(QtGui.QPolygonF(points), '')
# set its properties
contour.pen.setColor(QtGui.QColor(202, 31, 123, 255))
contour.pen.setWidthF(3.0)
contour.pen.setCosmetic(True)
contour.brush.setStyle(QtCore.Qt.NoBrush)
# add contour as a GraphicsItem to the scene
# these are the objects which are drawn in the GraphicsView
meshline = GraphicsItem.GraphicsItem(contour)
mesh_blocks.append(meshline)
airfoil.mesh_blocks = self.mainwindow.scene \
.createItemGroup(mesh_blocks)
# activate viewing options if mesh is created and displayed
self.mainwindow.centralwidget.cb8.setChecked(True)
self.mainwindow.centralwidget.cb8.setEnabled(True)
# after instantiating everything above switch it off
# as blocks should not be shown as a default
# now visibility of blocks fits to checkbox setting
self.mainwindow.centralwidget.cb8.click()
class BlockMesh:
def __init__(self, name='block'):
self.name = name
self.ULines = list()
def addLine(self, line):
# line is a list of (x, y) tuples
self.ULines.append(line)
def getULines(self):
return self.ULines
def getVLines(self):
vlines = list()
U, V = self.getDivUV()
# loop over all u-lines
for i in range(U + 1):
# prepare new v-line
vline = list()
# collect i-th point on each u-line
for uline in self.getULines():
vline.append(uline[i])
vlines.append(vline)
return vlines
def getLine(self, number=0, direction='u'):
if direction.lower() == 'u':
lines = self.getULines()
if direction.lower() == 'v':
lines = self.getVLines()
return lines[number]
def getDivUV(self):
u = len(self.getULines()[0]) - 1
v = len(self.getULines()) - 1
return u, v
def getNodeCoo(self, node):
I, J = node[0], node[1]
uline = self.getULines()[J]
point = uline[I]
return np.array(point)
def setNodeCoo(self, node, new_pos):
I, J = node[0], node[1]
uline = self.getULines()[J]
uline[I] = new_pos
return
@staticmethod
def makeLine(p1, p2, divisions=1, ratio=1.0):
vec = p2 - p1
dist = np.linalg.norm(vec)
spacing = BlockMesh.spacing(divisions=divisions,
ratio=ratio, length=dist)
line = list()
line.append((p1.tolist()[0], p1.tolist()[1]))
for i in range(1, len(spacing)):
p = p1 + spacing[i] * Utils.unit_vector(vec)
line.append((p.tolist()[0], p.tolist()[1]))
del line[-1]
line.append((p2.tolist()[0], p2.tolist()[1]))
return line
def extrudeLine_cell_thickness(self, line, cell_thickness=0.04,
growth=1.05,
divisions=1,
direction=3):
x, y = list(zip(*line))
x = np.array(x)
y = np.array(y)
if direction == 3:
spacing, _ = self.spacing_cell_thickness(
cell_thickness=cell_thickness,
growth=growth,
divisions=divisions)
normals = self.curveNormals(x, y)
for i in range(1, len(spacing)):
xo = x + spacing[i] * normals[:, 0]
yo = y + spacing[i] * normals[:, 1]
line = list(zip(xo.tolist(), yo.tolist()))
self.addLine(line)
elif direction == 4:
spacing, _ = self.spacing_cell_thickness(
cell_thickness=cell_thickness,
growth=growth,
divisions=divisions)
normals = self.curveNormals(x, y)
normalx = normals[:, 0].mean()
normaly = normals[:, 1].mean()
for i in range(1, len(spacing)):
xo = x + spacing[i] * normalx
yo = y + spacing[i] * normaly
line = list(zip(xo.tolist(), yo.tolist()))
self.addLine(line)
def extrudeLine(self, line, direction=0, length=0.1, divisions=1,
ratio=1.00001, constant=False):
x, y = list(zip(*line))
x = np.array(x)
y = np.array(y)
if constant and direction == 0:
x.fill(length)
line = list(zip(x.tolist(), y.tolist()))
self.addLine(line)
elif constant and direction == 1:
y.fill(length)
line = list(zip(x.tolist(), y.tolist()))
self.addLine(line)
elif direction == 3:
spacing = self.spacing(divisions=divisions,
ratio=ratio,
length=length)
normals = self.curveNormals(x, y)
for i in range(1, len(spacing)):
xo = x + spacing[i] * normals[:, 0]
yo = y + spacing[i] * normals[:, 1]
line = list(zip(xo.tolist(), yo.tolist()))
self.addLine(line)
elif direction == 4:
spacing = self.spacing(divisions=divisions,
ratio=ratio,
length=length)
normals = self.curveNormals(x, y)
normalx = normals[:, 0].mean()
normaly = normals[:, 1].mean()
for i in range(1, len(spacing)):
xo = x + spacing[i] * normalx
yo = y + spacing[i] * normaly
line = list(zip(xo.tolist(), yo.tolist()))
self.addLine(line)
def distribute(self, direction='u', number=0, type='constant'):
if direction == 'u':
line = np.array(self.getULines()[number])
elif direction == 'v':
line = np.array(self.getVLines()[number])
# interpolate B-spline through data points
# here, a linear interpolant is derived "k=1"
# splprep returns:
# tck ... tuple (t,c,k) containing the vector of knots,
# the B-spline coefficients, and the degree of the spline.
# u ... array of the parameters for each given point (knot)
tck, u = interpolate.splprep(line.T, s=0, k=1)
if type == 'constant':
t = np.linspace(0.0, 1.0, num=len(line))
if type == 'transition':
first = np.array(self.getULines()[0])
last = np.array(self.getULines()[-1])
tck_first, u_first = interpolate.splprep(first.T, s=0, k=1)
tck_last, u_last = interpolate.splprep(last.T, s=0, k=1)
if number < 0.0:
number = len(self.getVLines())
v = float(number) / float(len(self.getVLines()))
t = (1.0 - v) * u_first + v * u_last
# evaluate function at any parameter "0<=t<=1"
line = interpolate.splev(t, tck, der=0)
line = list(zip(line[0].tolist(), line[1].tolist()))
if direction == 'u':
self.getULines()[number] = line
elif direction == 'v':
for i, uline in enumerate(self.getULines()):
self.getULines()[i][number] = line[i]
@staticmethod
def spacing_cell_thickness(cell_thickness=0.04, growth=1.1, divisions=10):
# add cell thickness of first layer
spacing = [cell_thickness]
for i in range(divisions - 1):
spacing.append(spacing[0] + spacing[-1] * growth)
spacing.insert(0, 0.0)
length = np.sum(spacing)
return spacing, length
@staticmethod
def spacing(divisions=10, ratio=1.0, length=1.0):
"""Calculate point distribution on a line
Args:
divisions (int, optional): Number of subdivisions
ratio (float, optional): Ratio of last to first subdivision size
length (float, optional): length of line
Returns:
array: individual line segment lengths
"""
if divisions == 1:
sp = [0.0, 1.0]
return np.array(sp)
growth = ratio**(1.0 / (float(divisions) - 1.0))
if growth == 1.0:
growth = 1.0 + 1.0e-10
s = [1.0]
for i in range(1, divisions + 1):
s.append(growth**i)
spacing = np.array(s)
spacing -= spacing[0]
spacing /= spacing[-1]
spacing *= length
return spacing
def mapLines(self, line_1, line_2):
"""Map the distribution of points from one line to another line
Args:
line_1 (LIST): Source line (will be mapped)
line_2 (LIST): Destination line (upon this line_1 is mapped)
"""
pass
@staticmethod
def curveNormals(x, y, closed=False):
istart = 0
iend = 0
n = list()
for i, _ in enumerate(x):
if closed:
if i == len(x) - 1:
iend = -i - 1
else:
if i == 0:
istart = 1
if i == len(x) - 1:
iend = -1
a = np.array([x[i + 1 + iend] - x[i - 1 + istart],
y[i + 1 + iend] - y[i - 1 + istart]])
e = Utils.unit_vector(a)
n.append([e[1], -e[0]])
istart = 0
iend = 0
return np.array(n)
def transfinite(self, boundary=[], ij=[]):
"""Make a transfinite interpolation.
http://en.wikipedia.org/wiki/Transfinite_interpolation
upper
--------------------
| |
| |
left | | right
| |
| |
--------------------
lower
Example input for the lower boundary:
lower = [(0.0, 0.0), (0.1, 0.3), (0.5, 0.4)]
"""
if boundary:
lower = boundary[0]
upper = boundary[1]
left = boundary[2]
right = boundary[3]
elif ij:
lower = self.getULines()[ij[2]][ij[0]:ij[1] + 1]
upper = self.getULines()[ij[3]][ij[0]:ij[1] + 1]
left = self.getVLines()[ij[0]][ij[2]:ij[3] + 1]
right = self.getVLines()[ij[1]][ij[2]:ij[3] + 1]
else:
lower = self.getULines()[0]
upper = self.getULines()[-1]
left = self.getVLines()[0]
right = self.getVLines()[-1]
# FIXME
# FIXME left and right need to swapped from input
# FIXME
# FIXME like: left, right = right, left
# FIXME
lower = np.array(lower)
upper = np.array(upper)
left = np.array(left)
right = np.array(right)
# convert the block boundary curves into parametric form
# as curves need to be between 0 and 1
# interpolate B-spline through data points
# here, a linear interpolant is derived "k=1"
# splprep returns:
# tck ... tuple (t,c,k) containing the vector of knots,
# the B-spline coefficients, and the degree of the spline.
# u ... array of the parameters for each given point (knot)
tck_lower, u_lower = interpolate.splprep(lower.T, s=0, k=1)
tck_upper, u_upper = interpolate.splprep(upper.T, s=0, k=1)
tck_left, u_left = interpolate.splprep(left.T, s=0, k=1)
tck_right, u_right = interpolate.splprep(right.T, s=0, k=1)
nodes = np.zeros((len(left) * len(lower), 2))
# corner points
c1 = lower[0]
c2 = upper[0]
c3 = lower[-1]
c4 = upper[-1]
for i, xi in enumerate(u_lower):
for j, eta in enumerate(u_left):
node = i * len(u_left) + j
point = (1.0 - xi) * left[j] + xi * right[j] + \
(1.0 - eta) * lower[i] + eta * upper[i] - \
((1.0 - xi) * (1.0 - eta) * c1 + (1.0 - xi) * eta * c2 +
xi * (1.0 - eta) * c3 + xi * eta * c4)
nodes[node, 0] = point[0]
nodes[node, 1] = point[1]
vlines = list()
vline = list()
i = 0
for node in nodes:
i += 1
vline.append(node)
if i % len(left) == 0:
vlines.append(vline)
vline = list()
vlines.reverse()
if ij:
ulines = self.makeUfromV(vlines)
n = -1
for k in range(ij[2], ij[3] + 1):
n += 1
self.ULines[k][ij[0]:ij[1] + 1] = ulines[n]
else:
self.ULines = self.makeUfromV(vlines)
return
@staticmethod
def makeUfromV(vlines):
ulines = list()
uline = list()
for i in range(len(vlines[0])):
for vline in vlines:
x, y = vline[i][0], vline[i][1]
uline.append((x, y))
ulines.append(uline[::-1])
uline = list()
return ulines
@staticmethod
def writeFLMA(wind_tunnel, name='', depth=0.3):
basename = os.path.basename(name)
nameroot, extension = os.path.splitext(basename)
mesh = wind_tunnel.mesh
vertices, connectivity = mesh
with open(name, 'w') as f:
number_of_vertices_2D = len(vertices)
numvertex = '8'
# write number of points to FLMA file (*2 for z-direction)
f.write(str(2 * number_of_vertices_2D) + '\n')
signum = -1.
# write x-, y- and z-coordinates to FLMA file
# loop 1D direction (symmetry)
for _ in range(2):
for vertex in vertices:
f.write(str(vertex[0]) + ' ' + str(vertex[1]) +
' ' + str(signum * depth / 2.0) + ' ')
signum = 1.
# write number of cells to FLMA file
cells = len(connectivity)
f.write('\n' + str(cells) + '\n')
# write cell connectivity to FLMA file
for cell in connectivity:
cell_connect = str(cell[0]) + ' ' + \
str(cell[1]) + ' ' + \
str(cell[2]) + ' ' + \
str(cell[3]) + ' ' + \
str(cell[0] + number_of_vertices_2D) + ' ' + \
str(cell[1] + number_of_vertices_2D) + ' ' + \
str(cell[2] + number_of_vertices_2D) + ' ' + \
str(cell[3] + number_of_vertices_2D) + '\n'
f.write(numvertex + '\n')
f.write(cell_connect)
# FIRE element type (FET) for HEX element
fetHEX = '5'
f.write('\n' + str(cells) + '\n')
for i in range(cells):
f.write(fetHEX + ' ')
f.write('\n\n')
# FIRE element type (FET) for Quad element
fetQuad = '3\n'
# write FIRE selections to FLMA file
f.write('6\n')
f.write('right\n')
f.write(fetQuad)
f.write(str(2 * len(connectivity)) + '\n')
for i in range(len(connectivity)):
f.write(' %s 0' % (i))
f.write('\n')
f.write('\n')
f.write('left\n')
f.write(fetQuad)
f.write(str(2 * len(connectivity)) + '\n')
for i in range(len(connectivity)):
f.write(' %s 1' % (i))
f.write('\n')
f.write('\n')
f.write('bottom\n')
f.write(fetQuad)
f.write('2\n')
f.write('0 2\n')
f.write('\n')
f.write('top\n')
f.write(fetQuad)
f.write('2\n')
f.write('0 3\n')
f.write('\n')
f.write('back\n')
f.write(fetQuad)
f.write('2\n')
f.write('0 4\n')
f.write('\n')
f.write('front\n')
f.write(fetQuad)
f.write('2\n')
f.write('0 5\n')
logger.info('FIRE type mesh saved as {}'.
format(os.path.join(OUTPUTDATA, basename)))
@staticmethod
def writeSU2(wind_tunnel, name=''):
basename = os.path.basename(name)
nameroot, extension = os.path.splitext(basename)
mesh = wind_tunnel.mesh
blocks = wind_tunnel.blocks
boundary_loops = wind_tunnel.boundary_loops
vertices, connectivity = mesh
airfoil_subdivisions, v = blocks[0].getDivUV()
trailing_edge_subdivisions, _ = blocks[1].getDivUV()
# SU2 element types
element_type_quadrilateral = '9'
_date = date.today().strftime("%A %d. %B %Y")
with open(name, 'w') as f:
f.write('%\n')
f.write('% Airfoil contour: ' + nameroot + ' \n')
f.write('%\n')
f.write('% File created with ' + PyAero.__appname__ + '\n')
f.write('% Version: ' + PyAero.__version__ + '\n')
f.write('% Author: ' + PyAero.__author__ + '\n')
f.write('% Date: ' + _date + '\n')
# dimension of the problem
f.write('%\n')
f.write('% Problem dimension\n')
f.write('%\n')
f.write('NDIME= 2\n')
# element connectivity
f.write('%\n')
f.write('% Inner element connectivity\n')
f.write('%\n')
# number of elements
f.write('NELEM= %s\n' % (len(connectivity)))
for cell_id, cell in enumerate(connectivity):
cell_connect = element_type_quadrilateral + ' ' + \
str(cell[0]) + ' ' + \
str(cell[1]) + ' ' + \
str(cell[2]) + ' ' + \
str(cell[3]) + ' ' + str(cell_id) + '\n'
f.write(cell_connect)
# comment for vertices
f.write('%\n')
f.write('% Node coordinates\n')
f.write('%\n')
f.write('NPOIN=%s\n' % (len(vertices)))
# x- and y-coordinates
for node, vertex in enumerate(vertices):
x, y = vertex[0], vertex[1]
f.write(' {:24.16e} {:24.16e} {}\n'.format(x, y, node))
# boundaries
f.write('%\n')
f.write('% Boundary elements\n')
f.write('%\n')
# number of vertices
# number of marks (Airfoil, Farfield, Symmetry)
# f.write('NMARK= 3\n')
f.write('NMARK= 2\n')
# boundary definition (tag) for the airfoil
f.write('MARKER_TAG= {}\n'.format(wind_tunnel.boundary_airfoil))
f.write('MARKER_ELEMS= {}\n'.format(len(boundary_loops[0])))
for edge in boundary_loops[0]:
f.write('3 {} {}\n'.format(edge[0], edge[1]))
# boundary definition (tag) for the farfield
# this loops the complete outer boundary
f.write('MARKER_TAG= {}\n'.format(wind_tunnel.boundary_inlet))
f.write('MARKER_ELEMS= {}\n'.format(len(boundary_loops[1])))
for edge in boundary_loops[1]:
f.write('3 {} {}\n'.format(edge[0], edge[1]))
# boundary definition (tag) for the symmetry
# f.write('MARKER_TAG= {}\n'.format(wind_tunnel.boundary_symmetry))
# f.write('MARKER_ELEMS= {}\n'.format(len(connectivity)))
# for cell_id, cell in enumerate(connectivity):
# cell_connect = element_type_quadrilateral + ' ' + \
# str(cell[0]) + ' ' + \
# str(cell[1]) + ' ' + \
# str(cell[2]) + ' ' + \
# str(cell[3])
# f.write('{}\n'.format(cell_connect))
logger.info('SU2 type mesh saved as {}'.
format(name))
@staticmethod
def writeGMSH(wind_tunnel, name=''):
"""export mesh in GMSH format 2
http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-_0028Legacy_0029
Args:
mesh (TYPE): Description
blocks (TYPE): Description
name (str, optional): Description
"""
basename = os.path.basename(name)
nameroot, extension = os.path.splitext(basename)
mesh = wind_tunnel.mesh
boundary_loops = wind_tunnel.boundary_loops
bnd_airfoil = wind_tunnel.lineedit_airfoil
bnd_inlet = wind_tunnel.lineedit_inlet
bnd_outlet = wind_tunnel.lineedit_outlet
bnd_symmetry = wind_tunnel.lineedit_symmetry
is_outlet = wind_tunnel.is_outlet
vertices, connectivity = mesh
# element type "1" is GMSH 2-node line
# element type "2" is GMSH 3-node triangle
# element type "3" is GMSH 4-node quadrangle
element_type_line = '1'
# element_type_triangle = '2'
element_type_quadrangle = '3'
# write date in English
locale.setlocale(locale.LC_ALL, 'en')
_date = date.today().strftime("%A %d. %B %Y")
with open(name, 'w') as f:
f.write('$MeshFormat\n')
f.write('2.2 0 8\n')
f.write('$EndMeshFormat\n')
f.write('$Comments\n')
f.write(' Airfoil contour: ' + nameroot + ' \n')
f.write(' File created with ' + PyAero.__appname__ + '\n')
f.write(' Version: ' + PyAero.__version__ + '\n')
f.write(' Author: ' + PyAero.__author__ + '\n')
f.write(' Date: ' + _date + '\n')
f.write('$EndComments\n')
'''
$PhysicalNames
number-of-names
physical-dimension physical-tag "physical-name"
$EndPhysicalNames
'''
f.write('$PhysicalNames\n')
f.write('4\n')
f.write('1 1 "{}"\n'.format(bnd_airfoil))
f.write('1 2 "{}"\n'.format(bnd_inlet))
f.write('1 3 "{}"\n'.format(bnd_outlet))
f.write('2 4 "{}"\n'.format(bnd_symmetry))
f.write('$EndPhysicalNames\n')
f.write('$Nodes\n')
f.write('%s\n' % (len(vertices)))
# x- and y-coordinates
for node, vertex in enumerate(vertices, start=1):
x, y = vertex[0], vertex[1]
f.write(' {:} {:16.8} {:16.8} 0.0\n'.format(node, x, y))
f.write('$EndNodes\n')
'''
$Elements
number-of-elements
elm-number elm-type number-of-tags < tag > … node-number-list
$EndElements
'''
f.write('$Elements\n')
# boundary_loops is a disjoint set groups element
# key for each loop is one arbitrary vertex of the loop
# one loop is made by the airfoil
# the other loop is made by the windtunnel outer boundary
keys = list(boundary_loops.keys())
# print('Number of boundary loops', len(keys))
elements_loop1 = len(list(boundary_loops[keys[0]]))
elements_loop2 = len(list(boundary_loops[keys[1]]))
number_of_cells = len(connectivity)
# number of elements
# compiled of airfoil, outer boundary and mesh itself
f.write('{}\n'.format(elements_loop1 + elements_loop2 +
number_of_cells))
element_id = 0
# FIXME
# FIXME refactor dicts and their usage
# FIXME
# write boundary elements (as per physical names)
physical = {0: '1', 1: '2'}
elementary_entities = {0: '8', 1: '7'}
for j, loop in enumerate(boundary_loops):
for i, node in enumerate(boundary_loops[loop]):
element_id += 1
# an element consists of:
# element_id
# element_type
#
if is_outlet[i]:
physical_l = '3'
elementary_entities_l = '9'
else:
physical_l = physical[j]
elementary_entities_l = elementary_entities[j]
element = ' ' + str(element_id) + ' ' + \
element_type_line + ' 3 ' + physical_l + ' ' + \
elementary_entities_l + ' 0 ' + str(node[0] + 1) + \
' ' + str(node[1] + 1) + '\n'
f.write(element)
# write mesh elements
# includes physical tag for symmetry "4"
for cell in connectivity:
element_id += 1
element = ' ' + str(element_id) + ' ' + \
element_type_quadrangle + ' 3 4 6 0 ' + \
str(cell[0] + 1) + ' ' + \
str(cell[1] + 1) + ' ' + \
str(cell[2] + 1) + ' ' + \
str(cell[3] + 1) + '\n'
f.write(element)
f.write('$EndElements')
logger.info('GMSH type mesh saved as {}'.
format(os.path.join(OUTPUTDATA, basename)))
class Smooth:
def __init__(self, block):
self.block = block
def getNeighbours(self, node):
"""Get a list of neighbours around a node """
i, j = node[0], node[1]
neighbours = {1: (i - 1, j - 1), 2: (i, j - 1), 3: (i + 1, j - 1),
4: (i + 1, j), 5: (i + 1, j + 1), 6: (i, j + 1),
7: (i - 1, j + 1), 8: (i - 1, j)}
return neighbours
def smooth(self, nodes, iterations=1, algorithm='laplace'):
"""Smoothing of a square lattice mesh
Algorithms:
- Angle based
Tian Zhou:
AN ANGLE-BASED APPROACH TO TWO-DIMENSIONAL MESH SMOOTHING
- Laplace
Mean of surrounding node coordinates
- Parallelogram smoothing
Sanjay Kumar Khattri:
A NEW SMOOTHING ALGORITHM FOR QUADRILATERAL AND HEXAHEDRAL MESHES
Args:
nodes (TYPE): List of (i, j) tuples for the nodes to be smoothed
iterations (int, optional): Number of smoothing iterations
algorithm (str, optional): Smoothing algorithm
"""
# loop number of smoothing iterations
for i in range(iterations):
new_pos = list()
# smooth a node (i, j)
for node in nodes:
nb = self.getNeighbours(node)
if algorithm == 'laplace':
new_pos = (self.block.getNodeCoo(nb[2]) +
self.block.getNodeCoo(nb[4]) +
self.block.getNodeCoo(nb[6]) +
self.block.getNodeCoo(nb[8])) / 4.0
if algorithm == 'parallelogram':
new_pos = (self.block.getNodeCoo(nb[1]) +
self.block.getNodeCoo(nb[3]) +
self.block.getNodeCoo(nb[5]) +
self.block.getNodeCoo(nb[7])) / 4.0 - \
(self.block.getNodeCoo(nb[2]) +
self.block.getNodeCoo(nb[4]) +
self.block.getNodeCoo(nb[6]) +
self.block.getNodeCoo(nb[8])) / 2.0
if algorithm == 'angle_based':
pass
self.block.setNodeCoo(node, new_pos.tolist())
return self.block
def selectNodes(self, domain='interior', ij=[]):
"""Generate a node index list
Args:
domain (str, optional): Defines the part of the domain where
nodes shall be selected
Returns:
List: Indices as (i, j) tuples
"""
U, V = self.block.getDivUV()
nodes = list()
# select all nodes except boundary nodes
if domain == 'interior':
istart = 1
iend = U
jstart = 1
jend = V
if domain == 'ij':
istart = ij[0]
iend = ij[1]
jstart = ij[2]
jend = ij[3]
for i in range(istart, iend):
for j in range(jstart, jend):
nodes.append((i, j))
return nodes
class DisjointSet:
"""Summary
Attributes:
group (dict): Description
leader (dict): Description
oldgroup (dict): Description
oldleader (dict): Description
from: https://stackoverflow.com/a/3067672/2264936
"""
def __init__(self, size=None):
if size is None:
# maps a member to the group's leader
self.leader = {}
# maps a group leader to the group (which is a set)
self.group = {}
self.oldgroup = {}
self.oldleader = {}
else:
self.group = {i: set([i]) for i in range(0, size)}
self.leader = {i: i for i in range(0, size)}
self.oldgroup = {i: set([i]) for i in range(0, size)}
self.oldleader = {i: i for i in range(0, size)}
def add(self, a, b):
self.oldgroup = self.group.copy()
self.oldleader = self.leader.copy()
leadera = self.leader.get(a)
leaderb = self.leader.get(b)
if leadera is not None:
if leaderb is not None:
if leadera == leaderb:
return # nothing to do
groupa = self.group[leadera]
groupb = self.group[leaderb]
if len(groupa) < len(groupb):
a, leadera, groupa, b, leaderb, groupb = \
b, leaderb, groupb, a, leadera, groupa
groupa |= groupb
del self.group[leaderb]
for k in groupb:
self.leader[k] = leadera
else:
self.group[leadera].add(b)
self.leader[b] = leadera
else:
if leaderb is not None:
self.group[leaderb].add(a)
self.leader[a] = leaderb
else:
self.leader[a] = self.leader[b] = a
self.group[a] = set([a, b])
def connected(self, a, b):
leadera = self.leader.get(a)
leaderb = self.leader.get(b)
if leadera is not None:
if leaderb is not None:
return leadera == leaderb
else:
return False
else:
return False
def undo(self):
self.group = self.oldgroup.copy()
self.leader = self.oldleader.copy()
