made the get_phis and qlinear functions work independant of dimension

This commit is contained in:
Stephen Mardson McQuay 2011-02-02 10:56:20 -07:00
parent f1cd3061d7
commit 046e24b67d
4 changed files with 49 additions and 60 deletions

View File

@ -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,10 +33,22 @@ 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
# 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]
])
elif len(X) == 3:
A = np.array([
[ 1, 1, 1, 1 ],
[R[0][0], R[1][0], R[2][0], R[3][0]],
@ -66,15 +60,19 @@ def get_phis_3D(X, R):
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))

View File

@ -337,6 +337,7 @@ 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()
@ -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')

View File

@ -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)
r = True
if [i for i in phis if i < 0.0 - TOL]:

View File

@ -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)