From 9836fdcbfe4a720fab31bdc8954a7ae347e1f01d Mon Sep 17 00:00:00 2001
From: Stephen Mardson McQuay <stephen@mardson.com>
Date: Thu, 18 Mar 2010 23:32:24 -0600
Subject: [PATCH] added files. still not consistently getting better numbers

---
 bin/test.py                 | 38 +++++++++++++++++++
 lib/grid/DD.py              | 69 +++++++++++++++++++++++++++++++++
 lib/grid/DDD.py             | 76 +++++++++++++++++++++++++++++++++++++
 lib/grid/simplex.py         | 38 +++++++++++++++++++
 test/baker3D.test.py        | 38 +++++++++++++++++++
 tools/blender/maketet.py    | 42 ++++++++++++++++++++
 tools/blender/plot_cloud.py | 12 ++++++
 7 files changed, 313 insertions(+)
 create mode 100755 bin/test.py
 create mode 100644 lib/grid/DD.py
 create mode 100644 lib/grid/DDD.py
 create mode 100644 lib/grid/simplex.py
 create mode 100755 test/baker3D.test.py
 create mode 100644 tools/blender/maketet.py
 create mode 100644 tools/blender/plot_cloud.py

diff --git a/bin/test.py b/bin/test.py
new file mode 100755
index 0000000..2b09a10
--- /dev/null
+++ b/bin/test.py
@@ -0,0 +1,38 @@
+#!/usr/bin/python2.6
+
+import sys
+from grid.DDD import simple_random_grid
+from baker import get_phis_3D, run_baker_3D
+from baker.tools import exact_func_3D
+
+try:
+  total_points = int(sys.argv[1])
+except:
+  total_points = 100
+
+g = simple_random_grid(total_points)
+
+open('/tmp/for_qhull.txt', 'w').write(g.for_qhull())
+
+X = [0.5, 0.3, 0.4]
+
+R, S =  g.get_simplex_and_nearest_points(X, extra_points = 3, simplex_size=4)
+
+
+print "R", R
+print "S", S
+
+exact = exact_func_3D(X)
+print "exact solution: %0.6f" % exact
+
+print "phi values: ", get_phis_3D(X, R.points)
+r = run_baker_3D(X, R, S)
+
+print 'qlin' , r['qlin']
+print 'error', r['error']
+print 'final', r['final']
+
+if abs(r['final'] - exact) <= abs(r['qlin'] - exact):
+  print "win"
+else:
+  print "failure"
diff --git a/lib/grid/DD.py b/lib/grid/DD.py
new file mode 100644
index 0000000..23d31c1
--- /dev/null
+++ b/lib/grid/DD.py
@@ -0,0 +1,69 @@
+from grid import grid
+from baker.tools import exact_func
+
+import numpy as np
+
+
+class simple_rect_grid(grid):
+  def __init__(self, xres = 5, yres = 5):
+    xmin = -1.0
+    xmax =  1.0
+    xspan = xmax - xmin
+    xdel = xspan / float(xres - 1)
+
+    ymin = -1.0
+    ymay =  1.0
+    yspan = ymay - ymin
+    ydel = yspan / float(yres - 1)
+
+
+    points = []
+    q      = []
+    for x in xrange(xres):
+      cur_x = xmin + (x * xdel)
+      for y in xrange(yres):
+        cur_y = ymin + (y * ydel)
+        points.append([cur_x, cur_y])
+        q.append(exact_func(cur_x, cur_y))
+    grid.__init__(self, points, q)
+    self.construct_connectivity()
+
+  def for_qhull_generator(self):
+    """
+      this returns a generator that should be fed into qdelaunay
+    """
+
+    yield '2';
+    yield '%d' % len(self.points)
+
+    for p in self.points:
+      yield "%f %f" % (p[0], p[1])
+
+  def for_qhull(self):
+    """
+      this returns a single string that should be fed into qdelaunay
+    """
+    r  = '2\n'
+    r += '%d\n' % len(self.points)
+    for p in self.points:
+      r += "%f %f\n" % (p[0], p[1])
+    return r
+
+
+class simple_random_grid(simple_rect_grid):
+  def __init__(self, num_points = 10):
+    points = []
+    q = []
+
+    r = np.random
+
+    for i in xrange(num_points):
+      cur_x = r.rand()
+      cur_y = r.rand()
+
+      points.append([cur_x, cur_y])
+      q.append(exact_func(cur_x, cur_y))
+    grid.__init__(self, points, q)
+
+    self.points = np.array(self.points)
+    self.q      = np.array(self.q)
diff --git a/lib/grid/DDD.py b/lib/grid/DDD.py
new file mode 100644
index 0000000..01b3155
--- /dev/null
+++ b/lib/grid/DDD.py
@@ -0,0 +1,76 @@
+from grid import grid
+from baker.tools import exact_func_3D
+
+import numpy as np
+
+
+class simple_rect_grid(grid):
+  def __init__(self, xres = 5, yres = 5, zres = 5):
+    xmin = -1.0
+    xmax =  1.0
+    xspan = xmax - xmin
+    xdel = xspan / float(xres - 1)
+
+    ymin = -1.0
+    ymay =  1.0
+    yspan = ymay - ymin
+    ydel = yspan / float(yres - 1)
+
+    zmin = -1.0
+    zmaz =  1.0
+    zspan = zmaz - zmin
+    zdel = zspan / float(zres - 1)
+
+
+    points = []
+    q      = []
+    for x in xrange(xres):
+      cur_x = xmin + (x * xdel)
+      for y in xrange(yres):
+        cur_y = ymin + (y * ydel)
+        for z in xrange(zres):
+          cur_z = zmin + (z * zdel)
+          points.append([cur_x, cur_y, cur_z])
+          q.append(exact_func_3D((cur_x, cur_y, cur_z)))
+    grid.__init__(self, points, q)
+  # self.construct_connectivity()
+
+  def for_qhull_generator(self):
+    """
+      this returns a generator that should be fed into qdelaunay
+    """
+
+    yield '3';
+    yield '%d' % len(self.points)
+
+    for p in self.points:
+      yield "%f %f %f" % tuple(p)
+
+  def for_qhull(self):
+    """
+      this returns a single string that should be fed into qdelaunay
+    """
+    r  = '3\n'
+    r += '%d\n' % len(self.points)
+    for p in self.points:
+      r += "%f %f %f\n" % tuple(p)
+    return r
+
+class simple_random_grid(simple_rect_grid):
+  def __init__(self, num_points = 10):
+    points = []
+    q = []
+
+    r = np.random
+
+    for i in xrange(num_points):
+      cur_x = r.rand()
+      cur_y = r.rand()
+      cur_z = r.rand()
+
+      points.append([cur_x, cur_y, cur_z])
+      q.append(exact_func_3D((cur_x, cur_y, cur_z)))
+    grid.__init__(self, points, q)
+
+    self.points = np.array(self.points)
+    self.q      = np.array(self.q)
diff --git a/lib/grid/simplex.py b/lib/grid/simplex.py
new file mode 100644
index 0000000..b35d8a9
--- /dev/null
+++ b/lib/grid/simplex.py
@@ -0,0 +1,38 @@
+from baker import get_phis
+
+class face(object):
+  def __init__(self, name):
+    self.name      = name
+    self.verts     = []
+    self.neighbors = []
+
+  def add_vert(self, v):
+    """
+       v should be an index into grid.points
+    """
+    self.verts.append(v)
+
+  def add_neighbor(self, n):
+    """
+       reference to another face object
+    """
+    self.neighbors.append(n)
+
+  def contains(self, X, grid):
+    R    = [grid.points[i] for i in self.verts]
+
+    phis = get_phis(X, R)
+
+    r = True
+    if [i for i in phis if i < 0.0]:
+      r = False
+    return r
+
+  def __str__(self):
+    neighbors = [i.name for i in self.neighbors]
+    return '%s: verts: %s neighbors: [%s]' %\
+           (
+             self.name,
+             self.verts,
+             ", ".join(neighbors)
+           )
diff --git a/test/baker3D.test.py b/test/baker3D.test.py
new file mode 100755
index 0000000..af172c9
--- /dev/null
+++ b/test/baker3D.test.py
@@ -0,0 +1,38 @@
+#!/usr/bin/python
+
+import unittest
+import baker
+import grid
+
+import numpy as np
+import scipy.spatial
+
+class TestSequenceFunctions(unittest.TestCase):
+  def setUp(self):
+    self.X = [0.0, 0.0, 0.0]
+    self.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],
+             ]
+    self.q = [0.0, 0.0, 0.0, 4]
+
+
+  def testGetPhis3D(self):
+    result       = [round(i, 5) for i in baker.get_phis_3D(self.X, self.r)]
+    right_answer = [round(i, 5) for i in [0.25, 0.25, 0.25, 0.25]]
+
+    self.assertEqual(result, right_answer)
+
+
+  def testQlinear3D(self):
+    phi, result  = baker.qlinear_3D(self.X, grid.grid(self.r, self.q))
+    result       = round(result, 5)
+    right_answer = round(1.0, 5)
+    self.assertEqual(result, right_answer)
+
+
+if __name__ == '__main__':
+  suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
+  unittest.TextTestRunner(verbosity=2).run(suite)
diff --git a/tools/blender/maketet.py b/tools/blender/maketet.py
new file mode 100644
index 0000000..27a4913
--- /dev/null
+++ b/tools/blender/maketet.py
@@ -0,0 +1,42 @@
+from Blender import *
+import bpy
+
+import math
+
+phiaa    = -19.471220333
+
+r        = 1.0
+phia     = math.pi * phiaa / 180.0
+the120   = math.pi * 120.0 / 180.0
+
+v = [[0,0,0], [0,0,0], [0,0,0], [0,0,0]]
+
+v[0][0] = 0.0
+v[0][1] = 0.0
+v[0][2] = r
+the     = 0.0
+
+for i in range(1,4):
+  v[i][0] = r * math.cos(the) * math.cos(phia)
+  v[i][1] = r * math.sin(the) * math.cos(phia)
+  v[i][2] = r * math.sin(phia)
+  the     = the + the120
+
+
+print v
+
+# map vertices to 4 faces
+f = []
+f.append([0, 1, 2])
+f.append([0, 2, 3])
+f.append([0, 3, 1])
+f.append([1, 2, 3])
+
+
+me = bpy.data.meshes.new('points')
+
+me.verts.extend([i for i in v])
+me.faces.extend([i for i in f])
+
+scn = bpy.data.scenes.active
+ob  = scn.objects.new(me, 'points_obj')
diff --git a/tools/blender/plot_cloud.py b/tools/blender/plot_cloud.py
new file mode 100644
index 0000000..99ba17a
--- /dev/null
+++ b/tools/blender/plot_cloud.py
@@ -0,0 +1,12 @@
+from Blender import *
+import bpy
+
+from grid.test3d import simple_rect_grid
+
+g = simple_rect_grid(10, 10, 20)
+
+me = bpy.data.meshes.new('points')
+me.verts.extend(g.points)
+# me.faces.extend([i.verts for i in p.faces.itervalues()])
+scn = bpy.data.scenes.active
+ob  = scn.objects.new(me, 'points_obj')