Source code for openmdao.lib.geometry.bspline

import cPickle
import os.path

from numpy import linspace, hstack, dstack, less ,less_equal, logical_and, \
    array, empty, matrix, dot
    
from scipy.optimize import fsolve, newton
from scipy.sparse import csr_matrix

[docs]class Bspline(object): def __init__(self,controls,points,order=3): #controls and points are 2-d arrays of points self.controls = controls self.order = order self.degree = order-1 self.n = len(controls) self.knots = hstack(([0,]*(self.degree), hstack((linspace(0,1,self.n-self.order+2),[1,]*(self.degree))) )) self.__b_cache = {} #uses for memoizing the b_jn function self.max_x = max(points[:,0]) #see if we can h1 = str(hash(tuple(points.flatten()))).replace("-","n") h2 = str(hash(tuple(controls.flatten()))).replace("-","n") pkl_file_name = "%s__%s.bspline_pkl"%(h1,h2) pkl_folder = "pyBspline_pkl" pkl_file_name = os.path.join(pkl_folder,pkl_file_name) if not os.path.exists(pkl_folder): os.mkdir(pkl_folder) if os.path.exists(pkl_file_name): self.B = cPickle.load(open(pkl_file_name)) else: self.B = self._calc_jacobian(points) cPickle.dump(self.B,open(pkl_file_name,'w')) def _calc_jacobian(self,points): #pre-calculate the B matrix n_p = points.shape[0] B = matrix(empty((n_p,self.n))) #1 row per point, one column per control_point r = range(0,self.n) for i,p in enumerate(points): t = self.find(p[0]) for j in r: B[i,j] = self.b_jn_wrapper(j,self.degree,t) #self.B = csr_matrix(B) self.B = B return self.B
[docs] def calc(self,C,points=None): self.controls = C if points: self.B = self._calc_jacobian(points) return array(self.B.dot(C))
[docs] def find(self,X): """returns the parametric coordinate that matches the given x location""" try: return array([fsolve(lambda f: self(f)[:,0][0] - x,[x/self.max_x,],xtol=1e-5) for x in X])[:,0] except TypeError: return fsolve(lambda f: self(f)[:,0] - X,[X/self.max_x,],xtol=1e-5)
[docs] def b_jn(self,j,n,t): t_j = self.knots[j] t_j1 = self.knots[j+1] t_jn = self.knots[j+n] t_jn1 = self.knots[j+n+1] if n==0: return logical_and(less_equal(t_j,t),less(t,t_j1)) if t_jn-t_j: q1 = (t-t_j)/(t_jn-t_j) else: q1 = 0 if t_jn1-t_j1: q2 = (t_jn1-t)/(t_jn1-t_j1) else: q2 = 0 B = q1*self.b_jn(j,n-1,t) + q2*self.b_jn(j+1,n-1,t) return B #dirty hack to fix some weird corner case that shows up
[docs] def b_jn_wrapper(self,j,n,t): tuple_t = tuple(t) try: return self.__b_cache[(j,n,tuple_t)] except KeyError: B = self.b_jn(j,n,t) if j == self.n-1: try: B[t==1]=1 #anywhere t=1, set B = 1 except TypeError: B = 1 self.__b_cache[(j,n,tuple_t)] = B return B
def __call__(self,t): rng = range(0,self.n) b = [self.b_jn_wrapper(i,self.degree,t) for i in rng] X = dot(self.controls[:,0],b) Y = dot(self.controls[:,1],b) return dstack((X,Y))[0]
OpenMDAO Home