Source code for pygeo.pyGeo

# ======================================================================
#         Imports
# ======================================================================
import os
import copy
import numpy as np
from scipy import sparse
from scipy.sparse.linalg.dsolve import factorized
from pyspline import Curve, Surface
from pyspline.utils import openTecplot, writeTecplot2D, closeTecplot
from . import geo_utils
from .topology import SurfaceTopology
from baseclasses.utils import Error


[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") 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