Source code for pyspline.pyVolume

# External modules
import numpy as np
from scipy.sparse import linalg

# Local modules
from . import libspline
from .pyCurve import Curve
from .pySurface import Surface
from .utils import Error, _assembleMatrix, checkInput, closeTecplot, openTecplot, writeTecplot3D


[docs]class Volume(object): """ Create an instance of a b-spline surface. There are two ways to initialize the class * **Creation**: Create an instance of the Volume class directly by supplying the required information: kwargs MUST contain the following information: ``ku, kv, kw, tu, tv, tw, coef``. * **LMS/Interpolation**: Create an instance of the Volume class by using an interpolating spline to given data points or a LMS spline. The following keyword argument information is required: 1. ``ku`` and ``kv`` and ``kw`` Spline Orders 2. ``X`` real array size (Nu, Nv, Nw, nDim) of data to fit. **OR** 1. ``x`` (3D) and ``y`` (3D) ``z`` (3D) 3D volume interpolation/fitting 3. ``u``, ``v``, ``w`` real array of size (Nu, Nv, Nw). Optional 4. ``nCtlu``, ``nCtlv``, ``nCtlw``. Specify number of control points. Only for LMS fitting. Parameters ---------- ku : int Spline order in u kv : int Spline order in v kw : int Spline order in w nCtlu : int Number of control points in u nCtlv : int Number of control points in v nCtlw : int Number of control points in w coef : array, size (nCtlu, nCtl, nDim) b-spline coefficient array. tu : array, size(nCtlu + ku) knot array in u tv : array, size(nCtlv + kv) knot array in v tw : array, size(nCtlw + kw) knot array in w X : array, size (Nu, Nv, Nw, ndim) Full data array to fit x : array, size (Nu, Nv) Just x data to fit/interpolate y : array, size (Nu, Nv, Nw) Just y data to fit/interpolate z : array, size (Nu, Nv, Nw) Just z data to fit/interpolate u : array, size (Nu, Nv, Nw) Explict u parameters to use. Optional. v : array, size (Nu, Nv, Nw) Explict v parameters to use. Optional. w : array, size (Nu, Nv, Nw) Explict w parameters to use. Optional. nIter : int Number of Hoscheks parater corrections to run recompute : bool Specifies whether the actual fitting is completed. Notes ----- The orientation of the nodes, edges and faces for the volumes is given below:: NODES | EDGES | FACES 6 7| 5 | #-----------#| #------------# | #-----------# / / | /| /| | /| /| / / | / | / | | / | / | / / | 6/ | 7/ | | / | 1 / | / / | / |10 / |11| / | ---------- 5 / / | / | 4 / | | / | / |(back) #-----------# | #------------# | | #-----------# | 4 5 | | | | | | | | | | | | | | | | | | | | <-3 2 3| | | 1 | | |2-> | | | | #-----------#| | #------|-----# | | #-----|-----# / / | |8 / |9 / |4 ---------- | / / / | | / | / | | / | / / / | | /2 | /3 | | / 0 | / / / | | / | / | | / | / / / | |/ |/ | |/ |/ #-----------# | #------------# | #-----------# 0 1 | 0 | """ def __init__(self, recompute=True, **kwargs): self.faceSurfaces = [None, None, None, None, None, None] self.edgeCurves = [None, None, None, None, None, None, None, None, None, None, None, None] self.data = None self.udata = None self.vdata = None self.wdata = None if ( "ku" in kwargs and "kv" in kwargs and "kw" in kwargs and "tu" in kwargs and "tv" in kwargs and "tw" in kwargs and "coef" in kwargs ): self.X = None self.u = None self.v = None self.w = None self.U = None self.V = None self.W = None self.interp = False self.ku = checkInput(kwargs["ku"], "ku", int, 0) self.kv = checkInput(kwargs["kv"], "kv", int, 0) self.kw = checkInput(kwargs["kw"], "kw", int, 0) self.coef = checkInput(kwargs["coef"], "coef", float, 4) self.nCtlu = self.coef.shape[0] self.nCtlv = self.coef.shape[1] self.nCtlw = self.coef.shape[2] self.nDim = self.coef.shape[3] self.tu = checkInput(kwargs["tu"], "tu", float, 1, (self.nCtlu + self.ku,)) self.tv = checkInput(kwargs["tv"], "tv", float, 1, (self.nCtlv + self.kv,)) self.tw = checkInput(kwargs["tw"], "tw", float, 1, (self.nCtlw + self.kw,)) self.umin = self.tu[0] self.umax = self.tu[-1] self.vmin = self.tv[0] self.vmax = self.tv[-1] self.wmin = self.tw[0] self.wmax = self.tw[-1] self.origData = False self.setFaceSurfaces() self.setEdgeCurves() else: # We have LMS/Interpolate # Do some checking on the number of control points if not ( "ku" in kwargs and "kv" in kwargs and "kw" in kwargs and ( "X" in kwargs or "x" in kwargs or ("x" in kwargs and "y" in kwargs) or ("x" in kwargs and "y" in kwargs and "z" in kwargs) ) ): raise ValueError( "ku, kv, and X (or x, or x and y, or x and y and z MUST be defined for task lms or interpolate" ) if "X" in kwargs: self.X = np.array(kwargs["X"]) if len(self.X.shape) == 1: self.nDim = 1 else: self.nDim = self.X.shape[3] elif "x" in kwargs and "y" in kwargs and "z" in kwargs: x = checkInput(kwargs["x"], "x", float, 3) y = checkInput(kwargs["y"], "y", float, 3) z = checkInput(kwargs["z"], "z", float, 3) self.X = np.zeros((x.shape[0], x.shape[1], x.shape[2], 3)) self.X[:, :, :, 0] = x self.X[:, :, :, 1] = y self.X[:, :, :, 2] = z self.nDim = 3 elif "x" in kwargs and "y" in kwargs: x = checkInput(kwargs["x"], "x", float, 3) y = checkInput(kwargs["y"], "y", float, 3) self.X = np.zeros((x.shape[0], x.shape[1], x.shape[3], 3)) self.X[:, :, :, 0] = x self.X[:, :, :, 1] = y self.nDim = 2 elif "x" in kwargs: x = checkInput(kwargs["x"], "x", float, 3) self.X = np.zeros((x.shape[0], x.shape[1], x.shape[3], 3)) self.X[:, :, :, 0] = kwargs["x"] self.nDim = 1 # enf if self.Nu = self.X.shape[0] self.Nv = self.X.shape[1] self.Nw = self.X.shape[2] self.ku = checkInput(kwargs["ku"], "ku", int, 0) self.kv = checkInput(kwargs["kv"], "kv", int, 0) self.kw = checkInput(kwargs["kw"], "kw", int, 0) if "nCtlu" in kwargs and "nCtlv" in kwargs and "nCtlw" in kwargs: self.nCtlu = checkInput(kwargs["nCtlu"], "nCtlu", int, 0) self.nCtlv = checkInput(kwargs["nCtlv"], "nCtlv", int, 0) self.nCtlw = checkInput(kwargs["nCtlw"], "nCtlw", int, 0) self.interp = False else: self.nCtlu = self.Nu self.nCtlv = self.Nv self.nCtlw = self.Nw self.interp = True self.origData = True # Sanity Check on Inputs if self.nCtlu >= self.Nu: self.nCtlu = self.Nu if self.nCtlv >= self.Nv: self.nCtlv = self.Nv if self.nCtlw >= self.Nw: self.nCtlw = self.Nw # Sanity check to make sure k is less than N if self.Nu < self.ku: self.ku = self.Nu if self.Nv < self.kv: self.kv = self.Nv if self.Nw < self.kw: self.kw = self.Nw if self.nCtlu < self.ku: self.ku = self.nCtlu if self.nCtlv < self.kv: self.kv = self.nCtlv if self.nCtlw < self.kw: self.kw = self.nCtlw if "nIter" in kwargs: self.nIter = kwargs["nIter"] else: self.nIter = 1 if "u" in kwargs and "v" in kwargs and "w" in kwargs: self.u = checkInput(kwargs["u"], "u", float, 1) self.v = checkInput(kwargs["v"], "v", float, 1) self.w = checkInput(kwargs["w"], "w", float, 1) else: if self.nDim == 3: self.calcParameterization() else: Error( "Automatic parameterization of ONLY available for spatial data in 3 dimensions. " + "Please supply u and v key word arguments otherwise." ) self.umin = 0 self.umax = 1 self.vmin = 0 self.vmax = 1 self.wmin = 0 self.wmax = 1 self.calcKnots() self.setCoefSize() if recompute: self.recompute() # end if (Interpolation type)
[docs] def recompute(self): """Recompute the volume if any driving data has been modified""" self.setCoefSize() vals, rowPtr, colInd = libspline.volume_jacobian_wrap( self.U, self.V, self.W, self.tu, self.tv, self.tw, self.ku, self.kv, self.kw, self.nCtlu, self.nCtlv, self.nCtlw, ) N = _assembleMatrix(vals, colInd, rowPtr, (self.Nu * self.Nv * self.Nw, self.nCtlu * self.nCtlv * self.nCtlw)) if self.interp: # Factorize once for efficiency solve = linalg.dsolve.factorized(N) for idim in range(self.nDim): self.coef[:, :, :, idim] = solve(self.X[:, :, :, idim].flatten()).reshape( [self.nCtlu, self.nCtlv, self.nCtlw] ) else: solve = linalg.dsolve.factorized(N.transpose() * N) for idim in range(self.nDim): rhs = N.transpose() * self.X[:, :, :, idim].flatten() self.coef[:, :, :, idim] = solve(rhs).reshape([self.nCtlu, self.nCtlv, self.nCtlw]) self.setFaceSurfaces() self.setEdgeCurves()
def setCoefSize(self): self.coef = np.zeros((self.nCtlu, self.nCtlv, self.nCtlw, self.nDim))
[docs] def calcParameterization(self): """Compute distance based parametrization. Use the fortran function for this""" S, u, v, w = libspline.para3d(self.X.T) S = S.T self.u = u self.v = v self.w = w self.U = np.asarray(S[:, :, :, 0], order="c") self.V = np.asarray(S[:, :, :, 1], order="c") self.W = np.asarray(S[:, :, :, 2], order="c") return
[docs] def calcKnots(self): """Determine the knots depending on if it is inerpolated or an LMS fit""" if self.interp: self.tu = libspline.knots_interp(self.u, np.array([], "d"), self.ku) self.tv = libspline.knots_interp(self.v, np.array([], "d"), self.kv) self.tw = libspline.knots_interp(self.w, np.array([], "d"), self.kw) else: self.tu = libspline.knots_lms(self.u, self.nCtlu, self.ku) self.tv = libspline.knots_lms(self.v, self.nCtlv, self.kv) self.tw = libspline.knots_lms(self.w, self.nCtlw, self.kw)
[docs] def getValueCorner(self, corner): """Get the value of the volume spline on the corner. Parameters ---------- corner : int Index of corner, 0<=corner<=7 Returns ------- value : float Volume spline evaluation at corner. """ if corner not in range(0, 8): raise Error("Corner must be in range 0..7 inclusive8") if corner == 0: val = self.getValue(self.umin, self.vmin, self.wmin) elif corner == 1: val = self.getValue(self.umax, self.vmin, self.wmin) elif corner == 2: val = self.getValue(self.umin, self.vmax, self.wmin) elif corner == 3: val = self.getValue(self.umax, self.vmax, self.wmin) elif corner == 4: val = self.getValue(self.umin, self.vmin, self.wmax) elif corner == 5: val = self.getValue(self.umax, self.vmin, self.wmax) elif corner == 6: val = self.getValue(self.umin, self.vmax, self.wmax) elif corner == 7: val = self.getValue(self.umax, self.vmax, self.wmax) return val
[docs] def getOrigValueCorner(self, corner): """Get the value of the original spline data on the corner if it exists Parameters ---------- corner : int Index of corner, 0<=corner<=7 Returns ------- value : float Original data on corner. """ if corner not in range(0, 8): raise Error("Corner must be in range 0..7 inclusive") if corner == 0: val = self.X[0, 0, 0] elif corner == 1: val = self.X[-1, 0, 0] elif corner == 2: val = self.X[0, -1, 0] elif corner == 3: val = self.X[-1, -1, 0] elif corner == 4: val = self.X[0, 0, -1] elif corner == 5: val = self.X[-1, 0, -1] elif corner == 6: val = self.X[0, -1, -1] elif corner == 7: val = self.X[-1, -1, -1] return val
[docs] def getOrigValuesFace(self, face): """For a given face index, face, return the 4 corners and the values of the midpoints of the 4 edges on that face. Parameters ---------- face : int Index of face, 0<=face<=5 Returns ------- coords : array of size (8, ndim) The first 4 entries are the corner, and the last 4 are the midpoints. """ if face not in range(0, 6): raise Error("Face must be in range 0..5 inclusive") if np.mod(self.Nu, 2) == 1: midu = [(self.Nu - 1) // 2, (self.Nu - 1) // 2] else: midu = [self.Nu // 2, self.Nu // 2 - 1] if np.mod(self.Nv, 2) == 1: midv = [(self.Nv - 1) // 2, (self.Nv - 1) // 2] else: midv = [self.Nv // 2, self.Nv // 2 - 1] if np.mod(self.Nw, 2) == 1: midw = [(self.Nw - 1) // 2, (self.Nw - 1) // 2] else: midw = [self.Nw // 2, self.Nw // 2 - 1] if face == 0: values = [ self.X[0, 0, 0], self.X[-1, 0, 0], self.X[0, -1, 0], self.X[-1, -1, 0], 0.5 * (self.X[midu[0], 0, 0] + self.X[midu[1], 0, 0]), 0.5 * (self.X[midu[0], -1, 0] + self.X[midu[1], -1, 0]), 0.5 * (self.X[0, midv[0], 0] + self.X[0, midv[1], 0]), 0.5 * (self.X[-1, midv[0], 0] + self.X[-1, midv[1], 0]), ] elif face == 1: values = [ self.X[0, 0, -1], self.X[-1, 0, -1], self.X[0, -1, -1], self.X[-1, -1, -1], 0.5 * (self.X[midu[0], 0, -1] + self.X[midu[1], 0, -1]), 0.5 * (self.X[midu[0], -1, -1] + self.X[midu[1], -1, -1]), 0.5 * (self.X[0, midv[0], -1] + self.X[0, midv[1], -1]), 0.5 * (self.X[-1, midv[0], -1] + self.X[-1, midv[1], -1]), ] elif face == 2: values = [ self.X[0, 0, 0], self.X[0, -1, 0], self.X[0, 0, -1], self.X[0, -1, -1], 0.5 * (self.X[0, midv[0], 0] + self.X[0, midv[1], 0]), 0.5 * (self.X[0, midv[0], -1] + self.X[0, midv[1], -1]), 0.5 * (self.X[0, 0, midw[0]] + self.X[0, 0, midw[1]]), 0.5 * (self.X[0, -1, midw[0]] + self.X[0, -1, midw[1]]), ] elif face == 3: values = [ self.X[-1, 0, 0], self.X[-1, -1, 0], self.X[-1, 0, -1], self.X[-1, -1, -1], 0.5 * (self.X[-1, midv[0], 0] + self.X[-1, midv[1], 0]), 0.5 * (self.X[-1, midv[0], -1] + self.X[-1, midv[1], -1]), 0.5 * (self.X[-1, 0, midw[0]] + self.X[-1, 0, midw[1]]), 0.5 * (self.X[-1, -1, midw[0]] + self.X[-1, -1, midw[1]]), ] elif face == 4: values = [ self.X[0, 0, 0], self.X[-1, 0, 0], self.X[0, 0, -1], self.X[-1, 0, -1], 0.5 * (self.X[midu[0], 0, 0] + self.X[midu[1], 0, 0]), 0.5 * (self.X[midu[0], 0, -1] + self.X[midu[1], 0, -1]), 0.5 * (self.X[0, 0, midw[0]] + self.X[0, 0, midw[1]]), 0.5 * (self.X[-1, 0, midw[0]] + self.X[-1, 0, midw[1]]), ] elif face == 5: values = [ self.X[0, -1, 0], self.X[-1, -1, 0], self.X[0, -1, -1], self.X[-1, -1, -1], 0.5 * (self.X[midu[0], -1, 0] + self.X[midu[1], -1, 0]), 0.5 * (self.X[midu[0], -1, -1] + self.X[midu[1], -1, -1]), 0.5 * (self.X[0, -1, midw[0]] + self.X[0, -1, midw[1]]), 0.5 * (self.X[-1, -1, midw[0]] + self.X[-1, -1, midw[1]]), ] return np.array(values)
[docs] def getMidPointEdge(self, edge): """Get the midpoint of the edge using the original data. Parameters ---------- edge : int Edge index. Must be 0<edge<11. Returns ------- midpoint : array of length nDim Mid point of edge """ if np.mod(self.Nu, 2) == 1: midu = [(self.Nu - 1) // 2, (self.Nu - 1) // 2] else: midu = [self.Nu // 2, self.Nu // 2 - 1] if np.mod(self.Nv, 2) == 1: midv = [(self.Nv - 1) // 2, (self.Nv - 1) // 2] else: midv = [self.Nv // 2, self.Nv // 2 - 1] if np.mod(self.Nw, 2) == 1: midw = [(self.Nw - 1) // 2, (self.Nw - 1) // 2] else: midw = [self.Nw // 2, self.Nw // 2 - 1] if edge == 0: val = self.X[midu[0], 0, 0] + self.X[midu[1], 0, 0] elif edge == 1: val = self.X[midu[0], -1, 0] + self.X[midu[1], -1, 0] elif edge == 2: val = self.X[0, midv[0], 0] + self.X[0, midv[1], 0] elif edge == 3: val = self.X[-1, midv[0], 0] + self.X[-1, midv[1], 0] elif edge == 4: val = self.X[midu[0], 0, -1] + self.X[midu[1], 0, -1] elif edge == 5: val = self.X[midu[0], -1, -1] + self.X[midu[1], -1, -1] elif edge == 6: val = self.X[0, midv[0], -1] + self.X[0, midv[1], -1] elif edge == 7: val = self.X[-1, midv[0], -1] + self.X[-1, midv[1], -1] elif edge == 8: val = self.X[0, 0, midw[0]] + self.X[0, 0, midw[1]] elif edge == 9: val = self.X[-1, 0, midw[0]] + self.X[-1, 0, midw[1]] elif edge == 10: val = self.X[0, -1, midw[0]] + self.X[0, -1, midw[1]] elif edge == 11: val = self.X[-1, -1, midw[0]] + self.X[-1, -1, midw[1]] return val
[docs] def getMidPointFace(self, face): """Get the midpoint of the face using the original data. Parameters ---------- face : int Face index. Must be 0, 1, 2, 3, 4 or 5 Returns ------- midpoint : array of length nDim Mid point of face """ if face not in range(0, 6): raise Error("Face must be in range 0..5 inclusive") if not self.origData: raise Error("No original data for this surface") if np.mod(self.Nu, 2) == 1: midu = [(self.Nu - 1) // 2, (self.Nu - 1) // 2] else: midu = [self.Nu // 2, self.Nu // 2 - 1] if np.mod(self.Nv, 2) == 1: midv = [(self.Nv - 1) // 2, (self.Nv - 1) // 2] else: midv = [self.Nv // 2, self.Nv // 2 - 1] if np.mod(self.Nw, 2) == 1: midw = [(self.Nw - 1) // 2, (self.Nw - 1) // 2] else: midw = [self.Nw // 2, self.Nw // 2 - 1] if face == 0: val = 0.25 * ( self.X[midu[0], midv[0], 0] + self.X[midu[1], midv[0], 0] + self.X[midu[0], midv[1], 0] + self.X[midu[1], midv[1], 0] ) elif face == 1: val = 0.25 * ( self.X[midu[0], midv[0], -1] + self.X[midu[1], midv[0], -1] + self.X[midu[0], midv[1], -1] + self.X[midu[1], midv[1], -1] ) elif face == 2: val = 0.25 * ( self.X[0, midv[0], midw[0]] + self.X[0, midv[1], midw[0]] + self.X[0, midv[0], midw[1]] + self.X[0, midv[1], midw[1]] ) elif face == 3: val = 0.25 * ( self.X[-1, midv[0], midw[0]] + self.X[-1, midv[1], midw[0]] + self.X[-1, midv[0], midw[1]] + self.X[-1, midv[1], midw[1]] ) elif face == 4: val = 0.25 * ( self.X[midu[0], 0, midw[0]] + self.X[midu[1], 0, midw[0]] + self.X[midu[0], 0, midw[1]] + self.X[midu[1], 0, midw[1]] ) elif face == 5: val = 0.25 * ( self.X[midu[0], -1, midw[0]] + self.X[midu[1], -1, midw[0]] + self.X[midu[0], -1, midw[1]] + self.X[midu[1], -1, midw[1]] ) return val
[docs] def setFaceSurfaces(self): """Create face spline objects for each of the faces""" self.faceSurfaces[0] = Surface(ku=self.ku, kv=self.kv, tu=self.tu, tv=self.tv, coef=self.coef[:, :, 0, :]) self.faceSurfaces[1] = Surface(ku=self.ku, kv=self.kv, tu=self.tu, tv=self.tv, coef=self.coef[:, :, -1, :]) self.faceSurfaces[2] = Surface(ku=self.ku, kv=self.kw, tu=self.tu, tv=self.tw, coef=self.coef[:, 0, :, :]) self.faceSurfaces[3] = Surface(ku=self.ku, kv=self.kw, tu=self.tu, tv=self.tw, coef=self.coef[:, -1, :, :]) self.faceSurfaces[4] = Surface(ku=self.kv, kv=self.kw, tu=self.tv, tv=self.tw, coef=self.coef[0, :, :, :]) self.faceSurfaces[5] = Surface(ku=self.kv, kv=self.kw, tu=self.tv, tv=self.tw, coef=self.coef[-1, :, :, :])
[docs] def setEdgeCurves(self): """Create edge spline objects for each edge""" self.edgeCurves[0] = Curve(k=self.ku, t=self.tu, coef=self.coef[:, 0, 0]) self.edgeCurves[1] = Curve(k=self.ku, t=self.tu, coef=self.coef[:, -1, 0]) self.edgeCurves[2] = Curve(k=self.kv, t=self.tv, coef=self.coef[0, :, 0]) self.edgeCurves[3] = Curve(k=self.kv, t=self.tv, coef=self.coef[-1, :, 0]) self.edgeCurves[4] = Curve(k=self.ku, t=self.tu, coef=self.coef[:, 0, -1]) self.edgeCurves[5] = Curve(k=self.ku, t=self.tu, coef=self.coef[:, -1, -1]) self.edgeCurves[6] = Curve(k=self.kv, t=self.tv, coef=self.coef[0, :, -1]) self.edgeCurves[7] = Curve(k=self.kv, t=self.tv, coef=self.coef[-1, :, -1]) self.edgeCurves[8] = Curve(k=self.kw, t=self.tw, coef=self.coef[0, 0, :]) self.edgeCurves[9] = Curve(k=self.kw, t=self.tw, coef=self.coef[-1, 0, :]) self.edgeCurves[10] = Curve(k=self.kw, t=self.tw, coef=self.coef[0, -1, :]) self.edgeCurves[11] = Curve(k=self.kw, t=self.tw, coef=self.coef[-1, -1, :])
[docs] def getBasisPt(self, u, v, w, vals, istart, colInd, lIndex): """This function should only be called from pyBlock The purpose is to compute the basis function for a u, v, w point and add it to pyBlocks's global dPt/dCoef matrix. vals, rowPtr, colInd is the CSR data and lIndex in the local -> global mapping for this volume""" return libspline.getbasisptvolume( u, v, w, self.tu, self.tv, self.tw, self.ku, self.kv, self.kw, vals, colInd, istart, lIndex.T )
def __call__(self, u, v, w): """ Equivalant to getValue() """ return self.getValue(u, v, w)
[docs] def getValue(self, u, v, w): """Get the value at the volume points(s) u, v, w. This is the main evaluation routine for the volume object. Parameters ---------- u : scalar, vector or matrix or tensor of values u position v : scalar, vector or matrix or tensor of values v position w : scalar, vector or matrix or tensor of values w position Returns ------- values : scalar, vector, matrix or tensor of values The spline evaluation at (u, v, w) """ u = np.atleast_3d(u).T v = np.atleast_3d(v).T w = np.atleast_3d(w).T if not u.shape == v.shape == w.shape: raise Error("u and v must have the same shape") vals = libspline.eval_volume(u, v, w, self.tu, self.tv, self.tw, self.ku, self.kv, self.kw, self.coef.T) return vals.squeeze().T
[docs] def getValueEdge(self, edge, s): """Get the value at the volume points(s) u, v, w Parameters ---------- edge : int Index of edge. Must be between 0 and 11. s : float or array Parameter position(s) along edge to evaluate. Returns ------- values : array Array of values evaluated along edge. """ if edge == 0: u = s v = self.vmin w = self.wmin elif edge == 1: u = s v = self.vmax w = self.wmin elif edge == 2: u = self.umin v = s w = self.wmax elif edge == 3: u = self.umax v = s w = self.umin elif edge == 4: u = s v = self.vmin w = self.wmax elif edge == 5: u = s v = self.vmax w = self.wmax elif edge == 6: u = self.umin v = s w = self.wmax elif edge == 7: u = self.umax v = s w = self.wmax elif edge == 8: u = self.umin v = self.vmin w = s elif edge == 9: u = self.umax v = self.vmin w = s elif edge == 10: u = self.umin v = self.vmax w = s elif edge == 11: u = self.umax v = self.vmax w = s u = np.atleast_3d(u).T v = np.atleast_3d(v).T w = np.atleast_3d(w).T if not u.shape == v.shape == w.shape: raise Error("u, v, and w must have the same shape") vals = libspline.eval_volume(u, v, w, self.tu, self.tv, self.tw, self.ku, self.kv, self.kw, self.coef.T) return vals.squeeze().T
[docs] def getBounds(self): """Determine the extents of the volume Returns ------- xMin : array of length 3 Lower corner of the bounding box xMax : array of length 3 Upper corner of the bounding box """ if self.nDim != 3: raise Error("getBounds is only defined for nDim = 3") cx = self.coef[:, :, :, 0].flatten() cy = self.coef[:, :, :, 1].flatten() cz = self.coef[:, :, :, 2].flatten() Xmin = np.zeros(self.nDim) Xmin[0] = min(cx) Xmin[1] = min(cy) Xmin[2] = min(cz) Xmax = np.zeros(self.nDim) Xmax[0] = max(cx) Xmax[1] = max(cy) Xmax[2] = max(cz) return Xmin, Xmax
[docs] def projectPoint(self, x0, nIter=25, eps=1e-10, **kwargs): """ Project a point x0 or points x0 onto the volume and return parametric positions Parameters ---------- x0 : array of length 3 or array of size (N, 3) Points to embed in the volume. If the points do not **actually** lie in the volume, the closest point is returned nIter : int Maximum number of Newton iterations to perform eps : float Tolerance for the Newton iteration u : float or array of len(X0) Optional initial guess for u position. v : float or array of len(X0) Optional initial guess for v position. w : float or array of len(X0) Optional initial guess for w position. Returns ------- u : float or array of length N u parametric position of closest point v : float or array of length N v parametric position of closest point w : float or array of length N w parametric position of closest point D : float or array of length N Distance between projected point and closest point If the points are 'inside' the volume, D should be less than eps. """ x0 = np.atleast_2d(x0) if "u" in kwargs and "v" in kwargs and "w" in kwargs: u = np.atleast_1d(kwargs["u"]) v = np.atleast_1d(kwargs["v"]) w = np.atleast_1d(kwargs["w"]) else: u = -1 * np.ones(len(x0)) v = -1 * np.ones(len(x0)) w = -1 * np.ones(len(x0)) if not len(x0) == len(u) == len(v) == len(w): raise Error("The length of x0 and u, v, w must be the same") # If necessary get brute-force starting point if np.any(u < 0) or np.any(u > 1) or np.any(v < 0) or np.any(v > 1): self.computeData() u, v, w = libspline.point_volume_start(x0.real.T, self.udata, self.vdata, self.wdata, self.data.T) D = np.zeros_like(x0) for i in range(len(x0)): u[i], v[i], w[i], D[i] = libspline.point_volume( x0[i].real, self.tu, self.tv, self.tw, self.ku, self.kv, self.kw, self.coef.T, nIter, eps, u[i], v[i], w[i], ) return u.squeeze(), v.squeeze(), w.squeeze(), D.squeeze()
[docs] def computeData(self, recompute=False): """ Compute discrete data that is used for the Tecplot Visualization as well as the data for doing the brute-force checks Parameters ---------- recompute : bool If True, recompute the data even if it has already been computed. """ # Only recompute if it doesn't exist already if self.data is None or recompute: self.edgeCurves[0].calcInterpolatedGrevillePoints() self.udata = self.edgeCurves[0].sdata self.edgeCurves[2].calcInterpolatedGrevillePoints() self.vdata = self.edgeCurves[2].sdata self.edgeCurves[8].calcInterpolatedGrevillePoints() self.wdata = self.edgeCurves[8].sdata U = np.zeros((len(self.udata), len(self.vdata), len(self.wdata))) V = np.zeros((len(self.udata), len(self.vdata), len(self.wdata))) W = np.zeros((len(self.udata), len(self.vdata), len(self.wdata))) for i in range(len(self.udata)): for j in range(len(self.vdata)): for k in range(len(self.wdata)): U[i, j, k] = self.udata[i] V[i, j, k] = self.vdata[j] W[i, j, k] = self.wdata[k] self.data = self.getValue(U, V, W)
[docs] def insertKnot(self, direction, s, r): """ Insert a knot into the volume along either u, v, w: Parameters ---------- direction : str Parameteric direction to insert. Either 'u', 'v', or 'w' s : float Parametric position along 'direction' to insert r : int Desired number of times to insert. Returns ------- r : int The **actual** number of times the knot was inserted. """ if direction not in ["u", "v", "w"]: raise Error("Direction must be one of 'u' or 'v' or 'w'") s = checkInput(s, "s", float, 0) r = checkInput(r, "r", int, 0) if s <= 0.0: return if s >= 1.0: return # This is relatively inefficient, but we'll do it for # simplicity just call insertknot for each slab if direction == "u": # Insert once to know how many times it was actually inserted # so we know how big to make the new coef: actualR, tNew, coefNew, breakPt = libspline.insertknot(s, r, self.tu, self.ku, self.coef[:, 0, 0].T) newCoef = np.zeros((self.nCtlu + actualR, self.nCtlv, self.nCtlw, self.nDim)) for k in range(self.nCtlvw): for j in range(self.nCtlv): actualR, tNew, coefSlice, breakPt = libspline.insertknot( s, r, self.tu, self.ku, self.coef[:, j, k].T ) newCoef[:, j, k] = coefSlice[:, 0 : self.nCtlu + actualR].T self.tu = tNew[0 : self.nCtlu + self.ku + actualR] self.nCtlu = self.nCtlu + actualR elif direction == "v": actualR, tNew, coefNew, breakPt = libspline.insertknot(s, r, self.tv, self.kv, self.coef[0, :, 0].T) newCoef = np.zeros((self.nCtlu, self.nCtlv + actualR, self.nCtlw, self.nDim)) for k in range(self.nCtlw): for i in range(self.nCtlu): actualR, tNew, coefSlice, breakPt = libspline.insertknot( s, r, self.tv, self.kv, self.coef[i, :, k].T ) newCoef[i, :, k] = coefSlice[:, 0 : self.nCtlv + actualR].T self.tv = tNew[0 : self.nCtlv + self.kv + actualR] self.nCtlv = self.nCtlv + actualR elif direction == "w": actualR, tNew, coefNew, breakPt = libspline.insertknot(s, r, self.tw, self.kw, self.coef[0, 0, :].T) newCoef = np.zeros((self.nCtlu, self.nCtlv, self.nCtlw + actualR, self.nDim)) for j in range(self.nCtlv): for i in range(self.nCtlu): actualR, tNew, coefSlice, breakPt = libspline.insertknot( s, r, self.tw, self.kw, self.coef[i, j, :].T ) newCoef[i, j, :] = coefSlice[:, 0 : self.nCtlw + actualR].T self.tw = tNew[0 : self.nCtlw + self.kw + actualR] self.nCtlw = self.nCtlw + actualR self.coef = newCoef # break_pt is converted to zero based ordering here!!! return actualR, breakPt - 1
[docs] def writeTecplot(self, fileName, vols=True, coef=True, orig=False): """Write the volume to a tecplot data file. Parameters ---------- fileName : str Tecplot filename. Should end in .dat vols : bool Flag specifiying whether the interpolated volume should be used. This is usually True if you want to get an approximation of the entire volume. coef : bool Flag specifiying if the control points are to be plotted orig : bool Flag specifiying if original data (used for fitting) is to be included. If on original data exists, this argument is ignored. """ f = openTecplot(fileName, self.nDim) if vols: self.computeData() writeTecplot3D(f, "interpolated", self.data) if coef: writeTecplot3D(f, "control_pts", self.coef) if orig and self.origData: writeTecplot3D(f, "orig_data", self.X) closeTecplot(f)