From 046e24b67d352d7bd22a2b2a8b59772632462da6 Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Wed, 2 Feb 2011 10:56:20 -0700 Subject: [PATCH] made the get_phis and qlinear functions work independant of dimension --- interp/baker/__init__.py | 85 +++++++++++++++++----------------------- interp/grid/__init__.py | 5 ++- interp/grid/simplex.py | 9 ++--- test/baker3D.test.py | 10 ++--- 4 files changed, 49 insertions(+), 60 deletions(-) diff --git a/interp/baker/__init__.py b/interp/baker/__init__.py index a83e55d..ee4109f 100644 --- a/interp/baker/__init__.py +++ b/interp/baker/__init__.py @@ -7,7 +7,11 @@ from interp.tools import log def get_phis(X, R): """ - The get_phis function is used to get barycentric coordonites for a point on a triangle. + The get_phis function is used to get barycentric coordonites for a point on + a triangle or tetrahedron: + + + in 2D: X -- the destination point (2D) X = [0,0] @@ -15,31 +19,9 @@ def get_phis(X, R): r = [[-1, -1], [0, 2], [1, -1]] this will return [0.333, 0.333, 0.333] - """ - # baker: eq 7 - A = np.array([ - [ 1, 1, 1], - [R[0][0], R[1][0], R[2][0]], - [R[0][1], R[1][1], R[2][1]], - ]) - b = np.array([ 1, - X[0], - X[1] - ]) - try: - phi = np.linalg.solve(A,b) - except np.linalg.LinAlgError as e: - msg = "calculation of phis yielded a linearly dependant system (%s)" % e - log.error(msg) - raise Exception(msg) - phi = np.dot(np.linalg.pinv(A), b) - return phi - -def get_phis_3D(X, R): - """ - The get_phis function is used to get barycentric coordonites for a point on a tetrahedron. + in 3D: X -- the destination point (3D) X = [0,0,0] @@ -51,30 +33,46 @@ def get_phis_3D(X, R): [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672], ] - this (should) will return [0.25, 0.25, 0.25, 0.25] + this will return [0.25, 0.25, 0.25, 0.25] """ # baker: eq 7 - A = np.array([ - [ 1, 1, 1, 1 ], - [R[0][0], R[1][0], R[2][0], R[3][0]], - [R[0][1], R[1][1], R[2][1], R[3][1]], - [R[0][2], R[1][2], R[2][2], R[3][2]], + # TODO: perhaps also test R .. ? + if len(X) == 2: + A = np.array([ + [ 1, 1, 1], + [R[0][0], R[1][0], R[2][0]], + [R[0][1], R[1][1], R[2][1]], + ]) + b = np.array([ 1, + X[0], + X[1] ]) - b = np.array([ 1, - X[0], - X[1], - X[2] - ]) + elif len(X) == 3: + A = np.array([ + [ 1, 1, 1, 1 ], + [R[0][0], R[1][0], R[2][0], R[3][0]], + [R[0][1], R[1][1], R[2][1], R[3][1]], + [R[0][2], R[1][2], R[2][2], R[3][2]], + ]) + b = np.array([ 1, + X[0], + X[1], + X[2] + ]) + else: + raise Exception("inapropriate demension on X") + try: phi = np.linalg.solve(A,b) except np.linalg.LinAlgError as e: - log.error("calculation of phis yielded a linearly dependant system: %s" % e) + msg = "calculation of phis yielded a linearly dependant system (%s)" % e + log.error(msg) + raise Exception(msg) phi = np.dot(np.linalg.pinv(A), b) return phi - def qlinear(X, R): """ this calculates the linear portion of q from X to R @@ -90,18 +88,7 @@ def qlinear(X, R): qlin = np.sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) return phis, qlin -def qlinear_3D(X, R): - """ - this calculates the linear portion of q from X to R - - X = destination point - R = simplex points - q = CFD quantities of interest at the simplex points(R) - """ - - phis = get_phis_3D(X, R.verts) - qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) - return phis, qlin +qlinear_3D = qlinear def get_error(phi, R, S, order = 2): log.debug("len(phi): %d"% len(phi)) diff --git a/interp/grid/__init__.py b/interp/grid/__init__.py index e88cc03..d313bfa 100644 --- a/interp/grid/__init__.py +++ b/interp/grid/__init__.py @@ -337,12 +337,13 @@ class delaunay_grid(grid): """ log.debug('start') qdelaunay_string = get_qdelaunay_dump_str(self) + # log.debug(qdelaunay_string) cell_to_cells = [] for matcher in delaunay_grid.cell_re.finditer(qdelaunay_string): d = matcher.groupdict() cell_name = d['cell'] - verticies = d['verts'] + verticies = d['verts'] neighboring_cells = d['neigh'] cur_cell = cell(cell_name) @@ -355,9 +356,11 @@ class delaunay_grid(grid): nghbrs = [(cell_name, i) for i in neighboring_cells.split()] cell_to_cells.extend(nghbrs) + log.debug(cell_to_cells) for rel in cell_to_cells: if rel[1] in self.cells: self.cells[rel[0]].add_neighbor(self.cells[rel[1]]) + log.debug(self.cells) log.debug('end') diff --git a/interp/grid/simplex.py b/interp/grid/simplex.py index b2d857f..1086280 100644 --- a/interp/grid/simplex.py +++ b/interp/grid/simplex.py @@ -1,4 +1,4 @@ -from interp.baker import get_phis, get_phis_3D +from interp.baker import get_phis from interp.tools import log TOL = 1e-8 @@ -9,11 +9,10 @@ def contains(X, R): represented by a list of n-degree coordinates) it now correctly checks for 2/3-D verts + + TODO: write unit test ... """ - if len(R) == 3: - phis = get_phis(X, R) - elif len(R) == 4: - phis = get_phis_3D(X, R) + phis = get_phis(X, R) r = True if [i for i in phis if i < 0.0 - TOL]: diff --git a/test/baker3D.test.py b/test/baker3D.test.py index 5683821..433e3d0 100755 --- a/test/baker3D.test.py +++ b/test/baker3D.test.py @@ -1,7 +1,7 @@ #!/usr/bin/env python import unittest -from interp.baker import get_phis_3D, qlinear_3D +from interp.baker import get_phis, qlinear from interp.grid import grid import numpy as np @@ -19,16 +19,16 @@ class TestSequenceFunctions(unittest.TestCase): self.q = [0.0, 0.0, 0.0, 4] - def testGetPhis3D(self): - result = get_phis_3D(self.X, self.r) + def testGetPhis(self): + result = get_phis(self.X, self.r) right_answer = [0.25, 0.25, 0.25, 0.25] for a,b in zip(result, right_answer): self.assertAlmostEqual(a,b) - def testQlinear3D(self): - phi, result = qlinear_3D(self.X, grid(self.r, self.q)) + def testQlinear(self): + phi, result = qlinear(self.X, grid(self.r, self.q)) result = result right_answer = 1.0 self.assertAlmostEqual(result, right_answer)