diff --git a/pykat/utilities/hermite.py b/pykat/utilities/hermite.py index c9a5ea836f6f6b78abdeb53d8aca1eaa416c7acc..aa4152fa65f3b3351bdaff5ecaeac98870256294 100644 --- a/pykat/utilities/hermite.py +++ b/pykat/utilities/hermite.py @@ -14,7 +14,7 @@ def hermite(n, x): return (4.0 * x * x - 2.0) elif n == 3: - + return (8.0 * x * x * x - 12.0 * x) diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py index 8debfd706fa7379b980bbdcd0fc72e6fbece4522..22641cde471fbfd845110c58f22724a6803253b3 100644 --- a/pykat/utilities/maps.py +++ b/pykat/utilities/maps.py @@ -50,7 +50,7 @@ class surfacemap(object): @property def offset(self): - return np.array(self.step_size)*(self.center - np.array(self.size)/2) + return np.array(self.step_size)*(np.array(self.center) - 1/2. - np.array(self.size)/2.0) @property def ROMWeights(self): @@ -65,14 +65,53 @@ class surfacemap(object): raise BasePyKatException("Map type needs handling") - def generateROMWeights(self, isModeMatched=True, verbose=False): + def generateROMWeights(self, isModeMatched=True, verbose=False, interpolate=False, tolerance = 1e-12, sigma = 1, sort=False): + + if interpolate: + from scipy.interpolate import interp2d + import numpy as np + + D = interp2d(self.x, self.y, self.data, fill_value=0) + + # only want even number of data points spread equally + # about the axes + if self.size[0] % 2 == 0: + Nx = self.size[0] + else: + Nx = self.size[0]-1 + + if self.size[1] % 2 == 0: + Ny = self.size[1] + else: + Ny = self.size[1]-1 + + nx = np.linspace(min(self.x), max(self.x), Nx) + ny = np.linspace(min(self.y), max(self.y), Ny) + + data = D(nx-self.offset[0], ny-self.offset[0]) + + self.name += " [ROMHOM interpolated]" + + self.center = (np.array(data.shape)+1)/2.0 + + self.data = data + xm = self.x[self.x<0] ym = self.y[self.y<0] + symm = False + + if min(xm) == min(ym) and max(xm) == max(ym) and len(xm) == len(ym): + symm = True + EI = {} - if len(xm) > 0: EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched)) - if len(ym) > 0: EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched)) + EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma), sort=sort) + + if symm: + EI["ym"] = EI["xm"] + else: + EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma), sort=sort) EI["limits"] = EI["xm"].limits diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py index ae77c020c611dc66de1b6865c49823e1fa6d3447..350fb6a6e061e4d364a8f6f0a24295e0388a1f49 100644 --- a/pykat/utilities/romhom.py +++ b/pykat/utilities/romhom.py @@ -209,7 +209,7 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1): return ReducedBasis(RB=np.array(RB_matrix), limits=romlimits, x=x) -def makeEmpiricalInterpolant(RB): +def makeEmpiricalInterpolant(RB, sort=False): e_x = RB.RB x = RB.x @@ -244,9 +244,19 @@ def makeEmpiricalInterpolant(RB): V[k][j] = e_x[j][node_indices[k]] invV = inv(V[0:len(node_indices), 0:len(node_indices)]) - B = B_matrix(invV, e_x) + B = np.array(B_matrix(invV, e_x)) - return EmpiricalInterpolant(B=np.array(B), nodes=np.array(x_nodes, dtype=np.float64), node_indices=np.array(node_indices, dtype=np.int32), limits=RB.limits, x=RB.x) + node_indices = np.array(node_indices, dtype=np.int32) + nodes = np.array(x_nodes, dtype=np.float64) + + if sort: + print sort + ix = np.argsort(nodes) + nodes = nodes[ix] + node_indices = node_indices[ix] + B = B[ix, :] + + return EmpiricalInterpolant(B=B, nodes=nodes, node_indices=node_indices, limits=RB.limits, x=RB.x) def makeWeights(smap, EI, verbose=True):