# Standard Python modules
import copy
import os
# External modules
from baseclasses.utils import Error
import numpy as np
from pyspline import Curve, Surface
from pyspline.utils import closeTecplot, openTecplot, writeTecplot2D
from scipy import sparse
from scipy.sparse.linalg import factorized
# Local modules
from . import geo_utils
from .topology import SurfaceTopology
[docs]
class pyGeo:
"""
pyGeo is a (fairly) complete geometry surfacing engine. It
performs multiple functions including producing surfaces from
cross sections and globally fitting surfaces. The actual b-spline
surfaces are of the pySpline Surface type.
The initialization type, initType, specifies what type of
initialization will be used. There are currently 3 initialization
types: plot3d, iges, and liftingSurface
Parameters
----------
initType : str, {'plot3d', 'iges', 'liftingSurface'}
A key word defining how this geo object will be defined.
fileName : str
Used for 'plot3d' and 'iges' only. Name of file to load.
xsections : list of filenames
List of the cross section coordinate files. Length = N
scale : list or array
List of the scaling factors (chord) for cross sections.
Length = N
offset : List or array
List of x-y offset to apply *before* scaling. Length = N
Xsec : List or array
List of spatial coordinates as to the placement of
cross sections. Size (N, 3)
x, y, z : list or array
If Xsec is not given, x,y,z arrays each of length N can
be given individually.
rot : List or array
List of x-y-z rotations to apply to cross sections. Length = N
rotX, rotY, rotZ : list or arrays
Individual lists of x,y, and z rotations. Each of length N
nCtl : int
Number of control points to use for fitting. If it is None, local
interpolation is performed which typically yields better results when
a small number of sections. When trying to match wing geometries, a
larger number of sections are often required. In this case, the
number of control points also needs to be increased to prevent the
file from becoming overly large.
kSpan : int
The spline order in span-wise direction. 2 for linear, 3 for quadratic
and 4 for cubic
ku : int
Spline order in u (for plot3d input files only)
kv : int
Spline order in v (for plot3d input files only)
nCtlu : int
Number of control points in u (for plot3d input files only)
nCtlv : int
Number of control points in v (for plot3d input files only)
"""
def __init__(self, initType, *args, **kwargs):
self.initType = initType
print("pyGeo Initialization Type is: %s" % (initType))
# ------------------- pyGeo Class Attributes -----------------
self.topo = None # The topology of the surfaces
self.surfs = [] # The list of surface (pySpline surf)
# objects
self.nSurf = None # The total number of surfaces
self.coef = None # The global (reduced) set of control
# points
if initType == "plot3d":
self._readPlot3D(*args, **kwargs)
elif initType == "iges":
self._readIges(*args, **kwargs)
elif initType == "liftingSurface":
self._init_lifting_surface(*args, **kwargs)
elif initType == "create": # Don't do anything
pass
else:
raise Error("Unknown init type. Valid Init types are 'plot3d', 'iges' and 'liftingSurface'")
# ----------------------------------------------------------------------------
# Initialization Type Functions
# ----------------------------------------------------------------------------
def _readPlot3D(self, fileName, order="f", ku=4, kv=4, nCtlu=4, nCtlv=4):
"""Load a plot3D file and create the splines to go with each patch
Parameters
----------
fileName : str
File name to load. Should end in .xyz
order : str
'f' for fortran ordering (usual), 'c' for c ordering
ku : int
Spline order in u
kv : int
Spline order in v
nCtlu : int
Number of control points in u
nCtlv : int
Number of control points in v
"""
f = open(fileName)
binary = False
nSurf = geo_utils.readNValues(f, 1, "int", binary)[0]
sizes = geo_utils.readNValues(f, nSurf * 3, "int", binary).reshape((nSurf, 3))
# ONE of Patch Sizes index must be one
nPts = 0
for i in range(nSurf):
if sizes[i, 0] == 1: # Compress back to indices 0 and 1
sizes[i, 0] = sizes[i, 1]
sizes[i, 1] = sizes[i, 2]
elif sizes[i, 1] == 1:
sizes[i, 1] = sizes[i, 2]
elif sizes[i, 2] == 1:
pass
else:
raise Error("One of the plot3d indices must be 1")
nPts += sizes[i, 0] * sizes[i, 1]
surfs = []
for i in range(nSurf):
curSize = sizes[i, 0] * sizes[i, 1]
surfs.append(np.zeros([sizes[i, 0], sizes[i, 1], 3]))
for idim in range(3):
surfs[-1][:, :, idim] = geo_utils.readNValues(f, curSize, "float", binary).reshape(
(sizes[i, 0], sizes[i, 1]), order=order
)
f.close()
# Now create a list of spline surface objects:
self.surfs = []
self.surfs0 = surfs
# Note This doesn't actually fit the surfaces...just produces
# the parametrization and knot vectors
self.nSurf = nSurf
for isurf in range(self.nSurf):
self.surfs.append(Surface(X=surfs[isurf], ku=ku, kv=kv, nCtlu=nCtlu, nCtlv=nCtlv))
def _readIges(self, fileName):
"""Load a Iges file and create the splines to go with each patch
Parameters
----------
fileName : str
Name of file to load.
"""
f = open(fileName)
Ifile = []
for line in f:
line = line.replace(";", ",") # This is a bit of a hack...
Ifile.append(line)
f.close()
start_lines = int(Ifile[-1][1:8])
general_lines = int(Ifile[-1][9:16])
directory_lines = int(Ifile[-1][17:24])
# parameter_lines = int((Ifile[-1][25:32]))
# Now we know how many lines we have to deal with
dir_offset = start_lines + general_lines
para_offset = dir_offset + directory_lines
surf_list = []
# Directory lines is ALWAYS a multiple of 2
for i in range(directory_lines // 2):
# 128 is bspline surface type
if int(Ifile[2 * i + dir_offset][0:8]) == 128:
start = int(Ifile[2 * i + dir_offset][8:16])
num_lines = int(Ifile[2 * i + 1 + dir_offset][24:32])
surf_list.append([start, num_lines])
self.nSurf = len(surf_list)
print("Found %d surfaces in Iges File." % (self.nSurf))
self.surfs = []
for isurf in range(self.nSurf): # Loop over our patches
data = []
# Create a list of all data
# -1 is for conversion from 1 based (iges) to python
para_offset = surf_list[isurf][0] + dir_offset + directory_lines - 1
for i in range(surf_list[isurf][1]):
aux = Ifile[i + para_offset][0:69].split(",")
for j in range(len(aux) - 1):
data.append(float(aux[j]))
# Now we extract what we need
Nctlu = int(data[1] + 1)
Nctlv = int(data[2] + 1)
ku = int(data[3] + 1)
kv = int(data[4] + 1)
counter = 10
tu = data[counter : counter + Nctlu + ku]
counter += Nctlu + ku
tv = data[counter : counter + Nctlv + kv]
counter += Nctlv + kv
weights = data[counter : counter + Nctlu * Nctlv]
weights = np.array(weights)
if weights.all() != 1:
print("WARNING: Not all weight in B-spline surface are 1. A NURBS surface CANNOT be replicated exactly")
counter += Nctlu * Nctlv
coef = np.zeros([Nctlu, Nctlv, 3])
for j in range(Nctlv):
for i in range(Nctlu):
coef[i, j, :] = data[counter : counter + 3]
counter += 3
# Last we need the ranges
prange = np.zeros(4)
prange[0] = data[counter]
prange[1] = data[counter + 1]
prange[2] = data[counter + 2]
prange[3] = data[counter + 3]
# Re-scale the knot vectors in case the upper bound is not 1
tu = np.array(tu)
tv = np.array(tv)
if not tu[-1] == 1.0:
tu /= tu[-1]
if not tv[-1] == 1.0:
tv /= tv[-1]
self.surfs.append(Surface(ku=ku, kv=kv, tu=tu, tv=tv, coef=coef))
# Generate dummy data for connectivity to work
u = np.linspace(0, 1, 3)
v = np.linspace(0, 1, 3)
[V, U] = np.meshgrid(v, u)
self.surfs[-1].X = self.surfs[-1](U, V)
self.surfs[-1].Nu = 3
self.surfs[-1].Nv = 3
self.surfs[-1].origData = True
def _init_lifting_surface(
self,
xsections,
X=None,
x=None,
y=None,
z=None,
rot=None,
rotX=None,
rotY=None,
rotZ=None,
scale=None,
offset=None,
nCtl=None,
kSpan=3,
teHeight=None,
teHeightScaled=None,
thickness=None,
bluntTe=False,
roundedTe=False,
bluntTaperRange=0.1,
squareTeTip=True,
teScale=0.75,
tip="rounded",
tipScale=0.25,
leOffset=0.001,
teOffset=0.001,
spanTang=0.5,
upTang=0.5,
):
"""Create a lifting surface by distributing the cross
sections. See pyGeo module documentation for information on
the most commonly used options."""
if X is not None:
Xsec = np.array(X)
else:
# We have to use x, y, z
Xsec = np.vstack([x, y, z]).T
N = len(Xsec)
if rot is not None:
rot = np.array(rot)
else:
if rotX is None:
rotX = np.zeros(N)
if rotY is None:
rotY = np.zeros(N)
if rotZ is None:
rotZ = np.zeros(N)
rot = np.vstack([rotX, rotY, rotZ]).T
if offset is None:
offset = np.zeros((N, 2))
if scale is None:
scale = np.ones(N)
# Limit kSpan to 2 if we only have two cross section
if len(Xsec) == 2:
kSpan = 2
if bluntTe:
if teHeight is None and teHeightScaled is None:
raise Error("teHeight OR teHeightScaled must be supplied for bluntTe option")
if teHeight:
teHeight = np.atleast_1d(teHeight)
if len(teHeight) == 1:
teHeight = np.ones(N) * teHeight
teHeight /= scale
if teHeightScaled:
teHeight = np.atleast_1d(teHeightScaled)
if len(teHeight) == 1:
teHeight = np.ones(N) * teHeight
else:
teHeight = [None for i in range(N)]
# Load in and fit them all
curves = []
knots = []
for i in range(len(xsections)):
if xsections[i] is not None:
x, y = geo_utils.readAirfoilFile(
xsections[i], bluntTe, bluntThickness=teHeight[i], bluntTaperRange=bluntTaperRange
)
weights = np.ones(len(x))
weights[0] = -1
weights[-1] = -1
if nCtl is not None:
c = Curve(x=x, y=y, nCtl=nCtl, k=4, weights=weights)
else:
c = Curve(x=x, y=y, localInterp=True)
curves.append(c)
knots.append(c.t)
else:
curves.append(None)
# If we are fitting curves, blend knot vectors and recompute
if nCtl is not None:
newKnots = geo_utils.blendKnotVectors(knots, True)
for i in range(len(xsections)):
if curves[i] is not None:
curves[i].t = newKnots.copy()
curves[i].recompute(100, computeKnots=False)
# If we want a pinched tip will will zero everything here.
if tip == "pinched":
# Just zero out the last section in y
if curves[-1] is not None:
print("zeroing tip")
curves[-1].coef[:, 1] = 0
else:
# Otherwise do knot inserions
origKnots = [None for i in range(N)]
for i in range(N):
if curves[i] is not None:
origKnots[i] = curves[i].t.copy()
# First take all the knots of the first curve:
baseKnots = []
baseKnots.extend(origKnots[0])
# For the rest of the curves
for i in range(1, N):
if curves[i] is not None:
knots = origKnots[i]
# Search for all indices
indices = np.searchsorted(baseKnots, knots, side="left")
toInsert = []
# Now go over the indices and see if we need to add
for j in range(len(indices)):
if abs(baseKnots[indices[j]] - knots[j]) > 1e-12:
toInsert.append(knots[j])
# Finally add the new indices and resort
baseKnots.extend(toInsert)
baseKnots.sort()
# We have to know determine more information about the set
# of baseKnots: We want just a list of just the knot
# values and their multiplicity.
newKnots = []
mult = []
i = 0
Nmax = len(baseKnots)
while i < len(baseKnots):
curKnot = baseKnots[i]
j = 1
while i + j < Nmax and abs(baseKnots[i + j] - curKnot) < 1e-12:
j += 1
i += j
newKnots.append(curKnot)
mult.append(j)
# Now we have a knot vector that *ALL* curve *MUST* have
# to form our surface. So we loop back over the curves and
# insert the knots as necessary.
for i in range(N):
if curves[i] is not None:
for j in range(len(newKnots)):
if not newKnots[j] in curves[i].t:
curves[i].insertKnot(newKnots[j], mult[j])
# If we want a pinched tip will will zero everything here.
if tip == "pinched":
# Just zero out the last section in y
if curves[-1] is not None:
curves[-1].coef[:, 1] = 0
# Finally force ALL curve to have PRECISELY identical knots
for i in range(len(xsections)):
if curves[i] is not None:
curves[i].t = curves[0].t.copy()
newKnots = curves[0].t.copy()
# end if (nCtl is not none)
# Generate a curve from X just for the parametrization
Xcurve = Curve(X=Xsec, k=kSpan)
# Now blend the missing sections
print("Interpolating missing sections ...")
for i in range(len(xsections)):
if xsections[i] is None:
# Fist two curves bounding this unknown one:
for j in range(i, -1, -1):
if xsections[j] is not None:
istart = j
break
for j in range(i, len(xsections), 1):
if xsections[j] is not None:
iend = j
break
# Now generate blending parameter alpha
sStart = Xcurve.s[istart]
sEnd = Xcurve.s[iend]
s = Xcurve.s[i]
alpha = (s - sStart) / (sEnd - sStart)
coef = curves[istart].coef * (1 - alpha) + curves[iend].coef * (alpha)
curves[i] = Curve(coef=coef, k=4, t=newKnots.copy())
# end for (xsections)
# Before we continue the user may want to artificially scale
# the thickness of the sections. This is useful when a
# different airfoil thickness is desired than the actual
# airfoil coordinates.
if thickness is not None:
thickness = np.atleast_1d(thickness)
if len(thickness) == 1:
thickness = np.ones(len(thickness)) * thickness
for i in range(N):
# Only scale the interior control points; not the first and last
curves[i].coef[1:-1, 1] *= thickness[i]
# Now split each curve at uSplit which roughly corresponds to LE
topCurves = []
botCurves = []
uSplit = curves[0].t[(curves[0].nCtl + 4 - 1) // 2]
for i in range(len(xsections)):
c1, c2 = curves[i].splitCurve(uSplit)
topCurves.append(c1)
c2.reverse()
botCurves.append(c2)
# Note that the number of control points on the upper and
# lower surface MAY not be the same. We can fix this by doing
# more knot insertions.
knotsTop = topCurves[0].t.copy()
knotsBot = botCurves[0].t.copy()
print("Symmetrizing Knot Vectors ...")
eps = 1e-12
for i in range(len(knotsTop)):
# Check if knotsTop[i] is not in knots_bot to within eps
found = False
for j in range(len(knotsBot)):
if abs(knotsTop[i] - knotsBot[j]) < eps:
found = True
if not found:
# Add to all sections
for ii in range(len(xsections)):
botCurves[ii].insertKnot(knotsTop[i], 1)
for i in range(len(knotsBot)):
# Check if knotsBot[i] is not in knotsTop to within eps
found = False
for j in range(len(knotsTop)):
if abs(knotsBot[i] - knotsTop[j]) < eps:
found = True
if not found:
# Add to all sections
for ii in range(len(xsections)):
topCurves[ii].insertKnot(knotsBot[i], 1)
# We now have symmetrized knot vectors for the upper and lower
# surfaces. We will copy the vectors to make sure they are
# precisely the same:
for i in range(len(xsections)):
topCurves[i].t = topCurves[0].t.copy()
botCurves[i].t = topCurves[0].t.copy()
# Now we can set the surfaces
ncoef = topCurves[0].nCtl
coefTop = np.zeros((ncoef, len(xsections), 3))
coefBot = np.zeros((ncoef, len(xsections), 3))
for i in range(len(xsections)):
# Scale, rotate and translate the coefficients
coefTop[:, i, 0] = scale[i] * (topCurves[i].coef[:, 0] - offset[i, 0])
coefTop[:, i, 1] = scale[i] * (topCurves[i].coef[:, 1] - offset[i, 1])
coefTop[:, i, 2] = 0
coefBot[:, i, 0] = scale[i] * (botCurves[i].coef[:, 0] - offset[i, 0])
coefBot[:, i, 1] = scale[i] * (botCurves[i].coef[:, 1] - offset[i, 1])
coefBot[:, i, 2] = 0
for j in range(ncoef):
coefTop[j, i, :] = geo_utils.rotzV(coefTop[j, i, :], rot[i, 2] * np.pi / 180)
coefTop[j, i, :] = geo_utils.rotxV(coefTop[j, i, :], rot[i, 0] * np.pi / 180)
coefTop[j, i, :] = geo_utils.rotyV(coefTop[j, i, :], rot[i, 1] * np.pi / 180)
coefBot[j, i, :] = geo_utils.rotzV(coefBot[j, i, :], rot[i, 2] * np.pi / 180)
coefBot[j, i, :] = geo_utils.rotxV(coefBot[j, i, :], rot[i, 0] * np.pi / 180)
coefBot[j, i, :] = geo_utils.rotyV(coefBot[j, i, :], rot[i, 1] * np.pi / 180)
# Finally translate according to positions specified
coefTop[:, i, :] += Xsec[i, :]
coefBot[:, i, :] += Xsec[i, :]
# Set the two main surfaces
self.surfs.append(Surface(coef=coefTop, ku=4, kv=kSpan, tu=topCurves[0].t, tv=Xcurve.t))
self.surfs.append(Surface(coef=coefBot, ku=4, kv=kSpan, tu=botCurves[0].t, tv=Xcurve.t))
print("Computing TE surfaces ...")
if bluntTe:
if not roundedTe:
coef = np.zeros((len(xsections), 2, 3), "d")
coef[:, 0, :] = coefTop[0, :, :]
coef[:, 1, :] = coefBot[0, :, :]
self.surfs.append(Surface(coef=coef, ku=kSpan, kv=2, tu=Xcurve.t, tv=[0, 0, 1, 1]))
else:
coef = np.zeros((len(xsections), 4, 3), "d")
coef[:, 0, :] = coefTop[0, :, :]
coef[:, 3, :] = coefBot[0, :, :]
# We need to get the tangent for the top and bottom
# surface, multiply by a scaling factor and this gives
# us the two inner rows of control points
for j in range(len(xsections)):
projTop = coefTop[0, j] - coefTop[1, j]
projBot = coefBot[0, j] - coefBot[1, j]
projTop /= np.linalg.norm(projTop)
projBot /= np.linalg.norm(projBot)
curTeThick = np.linalg.norm(coefTop[0, j] - coefBot[0, j])
coef[j, 1] = coef[j, 0] + projTop * 0.5 * curTeThick * teScale
coef[j, 2] = coef[j, 3] + projBot * 0.5 * curTeThick * teScale
self.surfs.append(Surface(coef=coef, ku=kSpan, kv=4, tu=Xcurve.t, tv=[0, 0, 0, 0, 1, 1, 1, 1]))
self.nSurf = len(self.surfs)
print("Computing Tip surfaces ...")
# # Add on additional surfaces if required for a rounded pinch tip
if tip == "rounded":
# Generate the midpoint of the coefficients
midPts = np.zeros([ncoef, 3])
upVec = np.zeros([ncoef, 3])
dsNorm = np.zeros([ncoef, 3])
for j in range(ncoef):
midPts[j] = 0.5 * (coefTop[j, -1] + coefBot[j, -1])
upVec[j] = coefTop[j, -1] - coefBot[j, -1]
ds = 0.5 * ((coefTop[j, -1] - coefTop[j, -2]) + (coefBot[j, -1] - coefBot[j, -2]))
dsNorm[j] = ds / np.linalg.norm(ds)
# Generate "average" projection Vector
projVec = np.zeros((ncoef, 3), "d")
for j in range(ncoef):
offset = teOffset + (float(j) / (ncoef - 1)) * (leOffset - teOffset)
projVec[j] = dsNorm[j] * (np.linalg.norm(upVec[j] * tipScale + offset))
# Generate the tip "line"
tipLine = np.zeros([ncoef, 3])
for j in range(ncoef):
tipLine[j] = midPts[j] + projVec[j]
# Generate a k=4 (cubic) surface
coefTopTip = np.zeros([ncoef, 4, 3])
coefBotTip = np.zeros([ncoef, 4, 3])
for j in range(ncoef):
coefTopTip[j, 0] = coefTop[j, -1]
coefTopTip[j, 1] = coefTop[j, -1] + projVec[j] * spanTang
coefTopTip[j, 2] = tipLine[j] + upTang * upVec[j]
coefTopTip[j, 3] = tipLine[j]
coefBotTip[j, 0] = coefBot[j, -1]
coefBotTip[j, 1] = coefBot[j, -1] + projVec[j] * spanTang
coefBotTip[j, 2] = tipLine[j] - upTang * upVec[j]
coefBotTip[j, 3] = tipLine[j]
# Modify for square_te_tip... taper over last 20%
if squareTeTip and not roundedTe:
tipDist = geo_utils.eDist(tipLine[0], tipLine[-1])
for j in range(ncoef): # Going from back to front:
fraction = geo_utils.eDist(tipLine[j], tipLine[0]) / tipDist
if fraction < 0.10:
fact = (1 - fraction / 0.10) ** 2
omfact = 1.0 - fact
coefTopTip[j, 1] = (
fact * ((5.0 / 6.0) * coefTopTip[j, 0] + (1.0 / 6.0) * coefBotTip[j, 0])
+ omfact * coefTopTip[j, 1]
)
coefTopTip[j, 2] = (
fact * ((4.0 / 6.0) * coefTopTip[j, 0] + (2.0 / 6.0) * coefBotTip[j, 0])
+ omfact * coefTopTip[j, 2]
)
coefTopTip[j, 3] = (
fact * ((1.0 / 2.0) * coefTopTip[j, 0] + (1.0 / 2.0) * coefBotTip[j, 0])
+ omfact * coefTopTip[j, 3]
)
coefBotTip[j, 1] = (
fact * ((1.0 / 6.0) * coefTopTip[j, 0] + (5.0 / 6.0) * coefBotTip[j, 0])
+ omfact * coefBotTip[j, 1]
)
coefBotTip[j, 2] = (
fact * ((2.0 / 6.0) * coefTopTip[j, 0] + (4.0 / 6.0) * coefBotTip[j, 0])
+ omfact * coefBotTip[j, 2]
)
coefBotTip[j, 3] = (
fact * ((1.0 / 2.0) * coefTopTip[j, 0] + (1.0 / 2.0) * coefBotTip[j, 0])
+ omfact * coefBotTip[j, 3]
)
surfTopTip = Surface(coef=coefTopTip, ku=4, kv=4, tu=topCurves[0].t, tv=[0, 0, 0, 0, 1, 1, 1, 1])
surfBotTip = Surface(coef=coefBotTip, ku=4, kv=4, tu=botCurves[0].t, tv=[0, 0, 0, 0, 1, 1, 1, 1])
self.surfs.append(surfTopTip)
self.surfs.append(surfBotTip)
self.nSurf += 2
if bluntTe:
# This is the small surface at the trailing edge
# tip. There are a couple of different things that can
# happen: If we rounded TE we MUST have a
# rounded-spherical-like surface (second piece of code
# below). Otherwise we only have a surface if
# square_te_tip is false in which case a flat curved
# surface results.
if not roundedTe and not squareTeTip:
coef = np.zeros((4, 2, 3), "d")
coef[:, 0] = coefTopTip[0, :]
coef[:, 1] = coefBotTip[0, :]
self.surfs.append(Surface(coef=coef, ku=4, kv=2, tu=[0, 0, 0, 0, 1, 1, 1, 1], tv=[0, 0, 1, 1]))
self.nSurf += 1
elif roundedTe:
coef = np.zeros((4, 4, 3), "d")
coef[:, 0] = coefTopTip[0, :]
coef[:, 3] = coefBotTip[0, :]
# We will actually recompute the coefficients
# on the last sections since we need to do a
# couple of more for this surface
for i in range(4):
projTop = coefTopTip[0, i] - coefTopTip[1, i]
projBot = coefBotTip[0, i] - coefBotTip[1, i]
projTop /= np.linalg.norm(projTop)
projBot /= np.linalg.norm(projBot)
curTeThick = np.linalg.norm(coefTopTip[0, i] - coefBotTip[0, i])
coef[i, 1] = coef[i, 0] + projTop * 0.5 * curTeThick * teScale
coef[i, 2] = coef[i, 3] + projBot * 0.5 * curTeThick * teScale
self.surfs.append(
Surface(coef=coef, ku=4, kv=4, tu=[0, 0, 0, 0, 1, 1, 1, 1], tv=[0, 0, 0, 0, 1, 1, 1, 1])
)
self.nSurf += 1
# end if bluntTe
elif tip == "pinched":
pass
else:
print("No tip specified")
# Cheat and make "original data" so that the edge connectivity works
u = np.linspace(0, 1, 3)
v = np.linspace(0, 1, 3)
[V, U] = np.meshgrid(u, v)
for i in range(self.nSurf):
self.surfs[i].origData = True
self.surfs[i].X = self.surfs[i](U, V)
self.surfs[i].Nu = 3
self.surfs[i].Nv = 3
self._calcConnectivity(1e-6, 1e-6)
sizes = []
for isurf in range(self.nSurf):
sizes.append([self.surfs[isurf].nCtlu, self.surfs[isurf].nCtlv])
self.topo.calcGlobalNumbering(sizes)
self.setSurfaceCoef()
[docs]
def fitGlobal(self):
"""
Perform a global B-spline surface fit to determine the
coefficients of each patch. This is only used with an plot3D
init type
"""
print("Global Fitting")
nCtl = self.topo.nGlobal
print(" -> Copying Topology")
origTopo = copy.deepcopy(self.topo)
print(" -> Creating global numbering")
sizes = []
for isurf in range(self.nSurf):
sizes.append([self.surfs[isurf].Nu, self.surfs[isurf].Nv])
# Get the Global number of the original data
origTopo.calcGlobalNumbering(sizes)
N = origTopo.nGlobal
print(" -> Creating global point list")
pts = np.zeros((N, 3))
for ii in range(N):
pts[ii] = self.surfs[origTopo.gIndex[ii][0][0]].X[origTopo.gIndex[ii][0][1], origTopo.gIndex[ii][0][2]]
# Get the maximum k (ku, kv for each surf)
kmax = 2
for isurf in range(self.nSurf):
if self.surfs[isurf].ku > kmax:
kmax = self.surfs[isurf].ku
if self.surfs[isurf].kv > kmax:
kmax = self.surfs[isurf].kv
nnz = N * kmax * kmax
vals = np.zeros(nnz)
rowPtr = [0]
colInd = np.zeros(nnz, "intc")
for ii in range(N):
isurf = origTopo.gIndex[ii][0][0]
i = origTopo.gIndex[ii][0][1]
j = origTopo.gIndex[ii][0][2]
u = self.surfs[isurf].U[i, j]
v = self.surfs[isurf].V[i, j]
vals, colInd = self.surfs[isurf].getBasisPt(u, v, vals, rowPtr[ii], colInd, self.topo.lIndex[isurf])
kinc = self.surfs[isurf].ku * self.surfs[isurf].kv
rowPtr.append(rowPtr[-1] + kinc)
# Now we can crop out any additional values in col_ptr and vals
vals = vals[: rowPtr[-1]]
colInd = colInd[: rowPtr[-1]]
# Now make a sparse matrix
NN = sparse.csr_matrix((vals, colInd, rowPtr))
print(" -> Multiplying N^T * N")
NNT = NN.T
NTN = NNT * NN
print(" -> Factorizing...")
solve = factorized(NTN)
print(" -> Back Solving...")
self.coef = np.zeros((nCtl, 3))
for idim in range(3):
self.coef[:, idim] = solve(NNT * pts[:, idim])
print(" -> Setting Surface Coefficients...")
self._updateSurfaceCoef()
# ----------------------------------------------------------------------
# Topology Information Functions
# ----------------------------------------------------------------------
[docs]
def doConnectivity(self, fileName=None, nodeTol=1e-4, edgeTol=1e-4):
"""
This is the only public edge connectivity function.
If fileName exists it loads the file OR it calculates the connectivity
and saves to that file.
Parameters
----------
fileName : str
Filename for con file
nodeTol : float
The tolerance for identical nodes
edgeTol : float
The tolerance for midpoint of edges being identical
"""
if fileName is not None and os.path.isfile(fileName):
print("Reading Connectivity File: %s" % (fileName))
self.topo = SurfaceTopology(fileName=fileName)
if self.initType != "iges":
self._propagateKnotVectors()
sizes = []
for isurf in range(self.nSurf):
sizes.append([self.surfs[isurf].nCtlu, self.surfs[isurf].nCtlv])
self.surfs[isurf].recompute()
self.topo.calcGlobalNumbering(sizes)
else:
self._calcConnectivity(nodeTol, edgeTol)
sizes = []
for isurf in range(self.nSurf):
sizes.append([self.surfs[isurf].nCtlu, self.surfs[isurf].nCtlv])
self.topo.calcGlobalNumbering(sizes)
if self.initType != "iges":
self._propagateKnotVectors()
if fileName is not None:
print("Writing Connectivity File: %s" % (fileName))
self.topo.writeConnectivity(fileName)
if self.initType == "iges":
self.setSurfaceCoef()
def _calcConnectivity(self, nodeTol, edgeTol):
"""This function attempts to automatically determine the connectivity
between the patches"""
# Calculate the 4 corners and 4 midpoints for each surface
coords = np.zeros((self.nSurf, 8, 3))
for isurf in range(self.nSurf):
beg, mid, end = self.surfs[isurf].getOrigValuesEdge(0)
coords[isurf][0] = beg
coords[isurf][1] = end
coords[isurf][4] = mid
beg, mid, end = self.surfs[isurf].getOrigValuesEdge(1)
coords[isurf][2] = beg
coords[isurf][3] = end
coords[isurf][5] = mid
beg, mid, end = self.surfs[isurf].getOrigValuesEdge(2)
coords[isurf][6] = mid
beg, mid, end = self.surfs[isurf].getOrigValuesEdge(3)
coords[isurf][7] = mid
self.topo = SurfaceTopology(coords=coords, nodeTol=nodeTol, edgeTol=edgeTol)
[docs]
def printConnectivity(self):
"""
Print the Edge connectivity to the screen
"""
self.topo.printConnectivity()
def _propagateKnotVectors(self):
"""Propagate the knot vectors to make consistent"""
# First get the number of design groups
nDG = -1
ncoef = []
for i in range(self.topo.nEdge):
if self.topo.edges[i].dg > nDG:
nDG = self.topo.edges[i].dg
ncoef.append(self.topo.edges[i].N)
nDG += 1
for isurf in range(self.nSurf):
dgU = self.topo.edges[self.topo.edgeLink[isurf][0]].dg
dgV = self.topo.edges[self.topo.edgeLink[isurf][2]].dg
self.surfs[isurf].nCtlu = ncoef[dgU]
self.surfs[isurf].nCtlv = ncoef[dgV]
if self.surfs[isurf].ku < self.surfs[isurf].nCtlu:
if self.surfs[isurf].nCtlu > 4:
self.surfs[isurf].ku = 4
else:
self.surfs[isurf].ku = self.surfs[isurf].nCtlu
if self.surfs[isurf].kv < self.surfs[isurf].nCtlv:
if self.surfs[isurf].nCtlv > 4:
self.surfs[isurf].kv = 4
else:
self.surfs[isurf].kv = self.surfs[isurf].nCtlv
self.surfs[isurf].calcKnots()
# Now loop over the number of design groups, accumulate all
# the knot vectors that corresponds to this dg, then merge them all
for idg in range(nDG):
knotVectors = []
flip = []
for isurf in range(self.nSurf):
for iedge in range(4):
if self.topo.edges[self.topo.edgeLink[isurf][iedge]].dg == idg:
if self.topo.edgeDir[isurf][iedge] == -1:
flip.append(True)
else:
flip.append(False)
if iedge in [0, 1]:
knotVec = self.surfs[isurf].tu
elif iedge in [2, 3]:
knotVec = self.surfs[isurf].tv
if flip[-1]:
knotVectors.append((1 - knotVec)[::-1].copy())
else:
knotVectors.append(knotVec)
# end for (isurf)
# Now blend all the knot vectors
newKnotVec = geo_utils.blendKnotVectors(knotVectors, False)
newKnotVecFlip = (1 - newKnotVec)[::-1]
counter = 0
for isurf in range(self.nSurf):
for iedge in range(4):
if self.topo.edges[self.topo.edgeLink[isurf][iedge]].dg == idg:
if iedge in [0, 1]:
if flip[counter]:
self.surfs[isurf].tu = newKnotVecFlip.copy()
else:
self.surfs[isurf].tu = newKnotVec.copy()
elif iedge in [2, 3]:
if flip[counter]:
self.surfs[isurf].tv = newKnotVecFlip.copy()
else:
self.surfs[isurf].tv = newKnotVec.copy()
counter += 1
# ----------------------------------------------------------------------
# Surface Writing Output Functions
# ----------------------------------------------------------------------
[docs]
def writeTecplot(
self,
fileName,
orig=False,
surfs=True,
coef=True,
directions=False,
surfLabels=False,
edgeLabels=False,
nodeLabels=False,
):
"""Write the pyGeo Object to Tecplot dat file
Parameters
----------
fileName : str
File name for tecplot file. Should have .dat extension
surfs : bool
Flag to write discrete approximation of the actual surface
coef: bool
Flag to write b-spline coefficients
directions : bool
Flag to write surface direction visualization
surfLabels : bool
Flag to write file with surface labels
edgeLabels : bool
Flag to write file with edge labels
nodeLabels : bool
Falg to write file with node labels
"""
f = openTecplot(fileName, 3)
# Write out the Interpolated Surfaces
if surfs:
for isurf in range(self.nSurf):
self.surfs[isurf].computeData()
writeTecplot2D(f, "interpolated", self.surfs[isurf].data)
# Write out the Control Points
if coef:
for isurf in range(self.nSurf):
writeTecplot2D(f, "control_pts", self.surfs[isurf].coef)
# Write out the Original Data
if orig:
for isurf in range(self.nSurf):
writeTecplot2D(f, "control_pts", self.surfs[isurf].X)
# Write out The Surface Directions
if directions:
for isurf in range(self.nSurf):
self.surfs[isurf].writeDirections(f, isurf)
# Write out The Surface, Edge and Node Labels
dirName, fileName = os.path.split(fileName)
fileBaseName, _ = os.path.splitext(fileName)
if surfLabels:
# Split the filename off
labelFilename = dirName + "./" + fileBaseName + ".surf_labels.dat"
f2 = open(labelFilename, "w")
for isurf in range(self.nSurf):
midu = np.floor(self.surfs[isurf].nCtlu / 2)
midv = np.floor(self.surfs[isurf].nCtlv / 2)
textString = 'TEXT CS=GRID3D, X=%f, Y=%f, Z=%f, ZN=%d, T="S%d"\n' % (
self.surfs[isurf].coef[midu, midv, 0],
self.surfs[isurf].coef[midu, midv, 1],
self.surfs[isurf].coef[midu, midv, 2],
isurf + 1,
isurf,
)
f2.write("%s" % (textString))
f2.close()
if edgeLabels:
# Split the filename off
labelFilename = dirName + "./" + fileBaseName + "edge_labels.dat"
f2 = open(labelFilename, "w")
for iedge in range(self.topo.nEdge):
surfaces = self.topo.getSurfaceFromEdge(iedge)
pt = self.surfs[surfaces[0][0]].edgeCurves[surfaces[0][1]](0.5)
textString = 'TEXT CS=GRID3D X=%f, Y=%f, Z=%f, T="E%d"\n' % (pt[0], pt[1], pt[2], iedge)
f2.write("%s" % (textString))
f2.close()
if nodeLabels:
# First we need to figure out where the corners actually *are*
nNodes = len(geo_utils.unique(self.topo.nodeLink.flatten()))
nodeCoord = np.zeros((nNodes, 3))
for i in range(nNodes):
# Try to find node i
for isurf in range(self.nSurf):
if self.topo.nodeLink[isurf][0] == i:
coordinate = self.surfs[isurf].getValueCorner(0)
break
elif self.topo.nodeLink[isurf][1] == i:
coordinate = self.surfs[isurf].getValueCorner(1)
break
elif self.topo.nodeLink[isurf][2] == i:
coordinate = self.surfs[isurf].getValueCorner(2)
break
elif self.topo.nodeLink[isurf][3] == i:
coordinate = self.surfs[isurf].getValueCorner(3)
break
nodeCoord[i] = coordinate
# Split the filename off
labelFilename = dirName + "./" + fileBaseName + ".node_labels.dat"
f2 = open(labelFilename, "w")
for i in range(nNodes):
textString = 'TEXT CS=GRID3D, X=%f, Y=%f, Z=%f, T="n%d"\n' % (
nodeCoord[i][0],
nodeCoord[i][1],
nodeCoord[i][2],
i,
)
f2.write("%s" % (textString))
f2.close()
# Close out the file
closeTecplot(f)
[docs]
def writeIGES(self, fileName):
"""
Write the surface to IGES format
Parameters
----------
fileName : str
File name of iges file. Should have .igs extension.
"""
f = open(fileName, "w")
# Note: Eventually we may want to put the CORRECT Data here
f.write(" S 1\n")
f.write("1H,,1H;,7H128-000,11H128-000.IGS,9H{unknown},9H{unknown},16,6,15,13,15, G 1\n")
f.write("7H128-000,1.,6,1HM,8,0.016,15H19970830.165254, 0.0001,0., G 2\n")
f.write("[email protected],23HLegacy PDD AP Committee,11,3, G 3\n")
f.write("13H920717.080000,23HMIL-PRF-28000B0,CLASS 1; G 4\n")
Dcount = 1
Pcount = 1
for isurf in range(self.nSurf):
Pcount, Dcount = self.surfs[isurf].writeIGES_directory(f, Dcount, Pcount)
Pcount = 1
counter = 1
for isurf in range(self.nSurf):
Pcount, counter = self.surfs[isurf].writeIGES_parameters(f, Pcount, counter)
# Write the terminate statement
f.write("S%7dG%7dD%7dP%7d%40sT%6s1\n" % (1, 4, Dcount - 1, counter - 1, " ", " "))
f.close()
[docs]
def writeTin(self, fileName):
"""
Write the surfaces to ICEMCFD .tin format
Parameters
----------
fileName : str
File name of tin file. Should have .tin extension.
"""
f = open(fileName, "w")
# Standard Python modules
import datetime
# Write the required header info here:
f.write("// tetin file version 1.1\n")
f.write("// written by pyLayoutGeo - on %s\n" % (datetime.datetime.now()))
f.write("set_triangulation_tolerance 0.001\n")
# Now loop over the componets and each will write the info it
# has to the .tin file:
for i in range(self.nSurf):
if self.surfs[i].name is None:
name = "surface_%d" % i
else:
name = self.surfs[i].name
s = "define_surface name surf.%d family %s tetra_size %f\n" % (i, name, 1.0)
f.write(s)
self.surfs[i].writeTin(f)
# Write the closing info:
f.write("affix 0\n")
f.write("define_model 1e+10 reference_size 1\n")
f.write("return\n")
f.close()
# ----------------------------------------------------------------------
# Update and Derivative Functions
# ----------------------------------------------------------------------
def _updateSurfaceCoef(self):
"""Copy the pyGeo list of control points back to the surfaces"""
for ii in range(len(self.coef)):
for jj in range(len(self.topo.gIndex[ii])):
isurf = self.topo.gIndex[ii][jj][0]
i = self.topo.gIndex[ii][jj][1]
j = self.topo.gIndex[ii][jj][2]
self.surfs[isurf].coef[i, j] = self.coef[ii].astype("d")
for isurf in range(self.nSurf):
self.surfs[isurf].setEdgeCurves()
[docs]
def setSurfaceCoef(self):
"""Set the surface coef list from the pyspline surfaces"""
self.coef = np.zeros((self.topo.nGlobal, 3))
for isurf in range(self.nSurf):
surf = self.surfs[isurf]
for i in range(surf.nCtlu):
for j in range(surf.nCtlv):
self.coef[self.topo.lIndex[isurf][i, j]] = surf.coef[i, j]
[docs]
def getBounds(self, surfs=None):
"""Determine the extents of the collection of surfaces
Parameters
----------
surfs : list or array
Indices of surfaces defining subset for which to get the bounding
box
Returns
-------
xMin : array of length 3
Lower corner of the bounding box
xMax : array of length 3
Upper corner of the bounding box
"""
if surfs is None:
surfs = np.arange(self.nSurf)
Xmin0, Xmax0 = self.surfs[surfs[0]].getBounds()
for i in range(1, len(surfs)):
isurf = surfs[i]
Xmin, Xmax = self.surfs[isurf].getBounds()
# Now check them
if Xmin[0] < Xmin0[0]:
Xmin0[0] = Xmin[0]
if Xmin[1] < Xmin0[1]:
Xmin0[1] = Xmin[1]
if Xmin[2] < Xmin0[2]:
Xmin0[2] = Xmin[2]
if Xmax[0] > Xmax0[0]:
Xmax0[0] = Xmax[0]
if Xmax[1] > Xmax0[1]:
Xmax0[1] = Xmax[1]
if Xmax[2] > Xmax0[2]:
Xmax0[2] = Xmax[2]
return Xmin0, Xmax0
[docs]
def projectCurve(self, curve, *args, surfs=None, **kwargs):
"""
Project a pySpline curve onto the pyGeo object
Parameters
----------
curve : pySpline.Curve object
Curve to use for the intersection
surfs : list or array
Indices of surface defining subset for which to get the bounding
box
Returns
-------
Try it and see!
Notes
-----
This algorithm is not efficient at all. We basically do the
curve-surface projection algorithm for each surface the loop
back over them to see which is the best in terms of closest
distance. This intent is that the curve *actually* intersects
one of the surfaces.
"""
if surfs is None:
surfs = np.arange(self.nSurf)
temp = np.zeros((len(surfs), 4))
result = np.zeros((len(surfs), 4))
patchID = np.zeros(len(surfs), "intc")
for i in range(len(surfs)):
isurf = surfs[i]
u, v, s, d = self.surfs[isurf].projectCurve(curve, *args, **kwargs)
temp[i, :] = [u, v, s, np.linalg.norm(d)]
# Sort the results by distance
index = np.argsort(temp[:, 3])
for i in range(len(surfs)):
result[i] = temp[index[i]]
patchID[i] = surfs[index[i]]
return result, patchID
[docs]
def projectPoints(self, points, *args, surfs=None, **kwargs):
"""Project on or more points onto the nearest surface.
Parameters
----------
points : list or array
Singe point (size 3) or list of points size (N,3) points
to project onto the surfaces
surfs : list or array
Indices of surface defining subset for which to use for
projection
Returns
-------
u : float or array
u parameter values of closest point
v : float or array
v parameter values of closest point
PID : int or int array
Patch index corresponding to the u,v parameter values
"""
if surfs is None:
surfs = np.arange(self.nSurf)
N = len(points)
U = np.zeros((N, len(surfs)))
V = np.zeros((N, len(surfs)))
D = np.zeros((N, len(surfs), 3))
for i in range(len(surfs)):
isurf = surfs[i]
U[:, i], V[:, i], D[:, i, :] = self.surfs[isurf].projectPoint(points, *args, **kwargs)
u = np.zeros(N)
v = np.zeros(N)
patchID = np.zeros(N, "intc")
# Now post-process to get the lowest one
for i in range(N):
d0 = np.linalg.norm(D[i, 0])
u[i] = U[i, 0]
v[i] = V[i, 0]
patchID[i] = surfs[0]
for j in range(len(surfs)):
if np.linalg.norm(D[i, j]) < d0:
d0 = np.linalg.norm(D[i, j])
u[i] = U[i, j]
v[i] = V[i, j]
patchID[i] = surfs[j]
return u, v, patchID