a3d72ed200
I'm going to hurry and clean up this method, but for now it is verbose, and uses sets. i'm going to just check it in for future reference
161 lines
3.7 KiB
Python
161 lines
3.7 KiB
Python
import numpy as np
|
|
import sys
|
|
|
|
def get_phis(X, R):
|
|
"""
|
|
The get_phis function is used to get barycentric coordonites for a point on a triangle.
|
|
|
|
X -- the destination point (2D)
|
|
X = [0,0]
|
|
r -- the three points that make up the triangular simplex (2D)
|
|
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:
|
|
print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system"
|
|
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 triangle.
|
|
|
|
X -- the destination point (3D)
|
|
X = [0,0,0]
|
|
r -- the four points that make up the tetrahedron (3D)
|
|
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, 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]
|
|
])
|
|
try:
|
|
phi = np.linalg.solve(A,b)
|
|
except:
|
|
print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system"
|
|
phi = np.dot(np.linalg.pinv(A), b)
|
|
|
|
return phi
|
|
|
|
|
|
def qlinear(X, R, q):
|
|
"""
|
|
this calculates the linear portion of q from X to r
|
|
|
|
also, this is baker eq 3
|
|
|
|
X = destination point
|
|
R = simplex points
|
|
q = CFD quantities of interest at the simplex points
|
|
"""
|
|
|
|
phis = get_phis(X, R)
|
|
qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)])
|
|
return qlin
|
|
|
|
def qlinear_3D(X, R, q):
|
|
"""
|
|
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)
|
|
qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)])
|
|
return qlin
|
|
|
|
def run_baker(X, R, S):
|
|
"""
|
|
This is the main function to call to get an interpolation to X from the input meshes
|
|
|
|
X -- the destination point (2D)
|
|
X = [0,0]
|
|
|
|
R = Simplex
|
|
|
|
S = extra points
|
|
|
|
"""
|
|
|
|
# calculate values only for the triangle
|
|
phi = get_phis(X, R.points)
|
|
qlin = qlinear (X, R.points, R.q)
|
|
|
|
if len(S.points) == 0:
|
|
answer = {
|
|
'a': None,
|
|
'b': None,
|
|
'c': None,
|
|
'qlin': qlin,
|
|
'error': None,
|
|
'final': None,
|
|
}
|
|
return answer
|
|
|
|
B = [] # baker eq 9
|
|
w = [] # baker eq 11
|
|
|
|
for (s, q) in zip(S.points, S.q):
|
|
(phi1, phi2, phi3) = get_phis(s, R.points)
|
|
B.append([phi1 * phi2, phi2*phi3, phi3*phi1])
|
|
|
|
w.append(q - qlinear(s, R.points, R.q))
|
|
|
|
B = np.array(B)
|
|
w = np.array(w)
|
|
|
|
A = np.dot(B.T, B)
|
|
b = np.dot(B.T, w)
|
|
|
|
# baker solve eq 10
|
|
try:
|
|
(a, b, c) = np.linalg.solve(A,b)
|
|
except:
|
|
print >> sys.stderr, "warning: linear calculation went bad, resorting to np.linalg.pinv"
|
|
(a, b, c) = np.dot(np.linalg.pinv(A), b)
|
|
|
|
error_term = a * phi[0] * phi[1]\
|
|
+ b * phi[1] * phi[2]\
|
|
+ c * phi[2] * phi[0]
|
|
|
|
q_final = qlin + error_term
|
|
|
|
answer = {
|
|
'a': a,
|
|
'b': b,
|
|
'c': c,
|
|
'qlin': qlin,
|
|
'error': error_term,
|
|
'final': q_final,
|
|
}
|
|
|
|
return answer
|