working on getting baker to work in 3D. the code runs, but the numbers are odd. I suspect the rectangular grid, and am going to try a random cloud of points.

lib/baker/baker.py

@@ -41,9 +41,14 @@ def get_phis_3D(X, r):
     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]]
+         r = [
+               [0.0, 0.0, 1.0],
+               [0.94280904333606508, 0.0, -0.3333333283722672],
+               [-0.47140452166803232, 0.81649658244673617, -0.3333333283722672],
+               [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672],
+             ]
 
-    this will return [0.333, 0.333, 0.333]
+    this will return [0.25, 0.25, 0.25, 0.25]
   """
 
   # baker: eq 7
@@ -82,7 +87,7 @@ def qlinear(X, R):
   qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
   return phis, qlin
 
-def qlinear_3D(X, R, q):
+def qlinear_3D(X, R):
   """
     this calculates the linear portion of q from X to r
 
@@ -91,8 +96,8 @@ def qlinear_3D(X, R, q):
     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)])
+  phis = get_phis_3D(X, R.points)
+  qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
   return phis, qlin
 
 def run_baker(X, R, S):
@@ -103,13 +108,11 @@ def run_baker(X, R, S):
          X = [0,0]
 
     R = Simplex
-
     S = extra points
-
   """
 
   # calculate values only for the triangle
-  phi, qlin = qlinear (X, R)
+  phi, qlin = qlinear(X, R)
 
   if len(S.points) == 0:
     answer = {
@@ -129,7 +132,13 @@ def run_baker(X, R, S):
     cur_phi, cur_qlin = qlinear(s, R)
     (phi1, phi2, phi3) = cur_phi
 
-    B.append([phi1 * phi2, phi2 * phi3, phi3 * phi1])
+    B.append(
+              [
+                phi1 * phi2,
+                phi2 * phi3,
+                phi3 * phi1,
+              ]
+            )
     w.append(q - cur_qlin)
 
   B = np.array(B)
@@ -161,3 +170,87 @@ def run_baker(X, R, S):
             }
 
   return answer
+
+def run_baker_3D(X, R, S):
+  """
+    This is the main function to call to get an interpolation to X from the input meshes
+
+    X -- the destination point (3D)
+         X = [0,0,0]
+
+    R = Simplex (4 points, contains X)
+    S = extra points (surrounding, in some manner, R and X, but not in R)
+  """
+
+  # calculate values only for the triangle
+  phi, qlin = qlinear_3D(X, R)
+
+  if len(S.points) == 0:
+    answer = {
+               'a': None,
+               'b': None,
+               'c': None,
+               'd': None,
+               'e': None,
+               'f': 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):
+    cur_phi, cur_qlin = qlinear_3D(s, R)
+    (phi1, phi2, phi3, phi4) = cur_phi
+
+    B.append(
+              [
+                phi1 * phi2,
+                phi1 * phi3,
+                phi1 * phi4,
+                phi2 * phi3,
+                phi2 * phi4,
+                phi3 * phi4,
+              ]
+            )
+
+    w.append(q - cur_qlin)
+
+  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, d, e, f) = np.linalg.solve(A,b)
+  except:
+    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv"
+    (a, b, c, d, e, f) = np.dot(np.linalg.pinv(A), b)
+
+  error_term =  a * phi[0] * phi[1]\
+              + b * phi[0] * phi[2]\
+              + c * phi[0] * phi[3]\
+              + d * phi[1] * phi[2]\
+              + e * phi[1] * phi[3]\
+              + f * phi[2] * phi[3]
+
+  q_final = qlin + error_term
+
+  answer = {
+             'a': a,
+             'b': b,
+             'c': c,
+             'd': d,
+             'e': e,
+             'f': f,
+             'qlin': qlin,
+             'error': error_term,
+             'final': q_final,
+            }
+
+  return answer

lib/baker/tools.py

@@ -24,3 +24,12 @@ def exact_func(x, y):
     the exact function used from baker's article (for testing)
   """
   return np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2)
+
+def exact_func_3D(X):
+  """
+    the exact function (3D) used from baker's article (for testing)
+  """
+  x = X[0]
+  y = X[1]
+  z = X[2]
+  return np.power((np.sin(x * np.pi / 2.0) * np.sin(y * np.pi / 2.0) * np.sin(z * np.pi / 2.0)), 2)