Stephen Mardson McQuay 9 years ago
commit
961b3b7b26

+ 1
- 0
.hgignore View File

@@ -0,0 +1 @@
1
+\.pyc$

+ 2
- 0
.setup View File

@@ -0,0 +1,2 @@
1
+export PATH=~/src/baker/bin:~/src/baker/test:$PATH
2
+export PYTHONPATH=~/src/baker/lib

+ 43
- 0
bin/driver-random.py View File

@@ -0,0 +1,43 @@
1
+#!/usr/bin/python
2
+
3
+import sys
4
+from optparse import OptionParser
5
+
6
+import numpy as np
7
+import scipy.spatial
8
+
9
+from grid import simple_random_grid, exact_func
10
+import baker
11
+
12
+def rms(errors):
13
+  r = 0.0
14
+  for i in errors:
15
+    r += np.power(i, 2)
16
+  r = np.sqrt(r / len(errors))
17
+  return r
18
+
19
+if __name__ == '__main__':
20
+  parser = OptionParser()
21
+  parser.add_option("-o", "--output-file",      dest="output",           type='str', default = '/tmp/for_qhull.txt',       help = "qhull output file")
22
+  parser.add_option("-e", "--extra-points",     dest="extra",            type='int', default = 3,       help = "how many extra points")
23
+  parser.add_option("-s", "--source-total",     dest="source_total",     type='int', default = 100,     help = "total number of source points")
24
+  parser.add_option("-d", "--destination-total",dest="destination_total",type='int', default = 100,     help = "total number of destination points")
25
+  parser.add_option("-v", "--verbose", action = 'store_true', default = False, help = "verbosity")
26
+
27
+  (options, args) = parser.parse_args()
28
+  print >> sys.stderr, options
29
+
30
+  mesh_source = simple_random_grid(options.source_total)
31
+  tree = scipy.spatial.KDTree(mesh_source.points)
32
+  mesh_dest   = simple_random_grid(options.destination_total)
33
+
34
+  open(options.output, 'w').write(mesh_source.for_qhull())
35
+  print >> sys.stderr, "wrote output to %s" % options.output
36
+
37
+  errors = []
38
+  for x in mesh_dest.points:
39
+    (final, exact) = baker.run_baker(x, mesh_source, tree, options.extra, options.verbose)
40
+    cur_error = np.abs(final - exact)
41
+    errors.append(cur_error)
42
+
43
+  print rms(errors)

+ 45
- 0
bin/driver-structured.py View File

@@ -0,0 +1,45 @@
1
+#!/usr/bin/python
2
+
3
+import sys
4
+from optparse import OptionParser
5
+
6
+import numpy as np
7
+import scipy.spatial
8
+
9
+from grid import simple_rect_grid, exact_func
10
+import baker
11
+from tools import rms
12
+
13
+
14
+if __name__ == '__main__':
15
+  parser = OptionParser()
16
+  parser.add_option("-e", "--extra-points",   dest="extra",   type='int',  default = 3,      help = "how many extra points")
17
+  parser.add_option("-s", "--source-density", dest="source",  type='int',  default = 11,     help = "resolution of source mesh")
18
+  parser.add_option("-d", "--dest-density",   dest="dest",    type='int',  default = 11,     help = "resolution of dest mesh")
19
+  parser.add_option("-t", "--random-total",   dest="random_total",    type='int',  default = 100,     help = "total number of random points")
20
+  parser.add_option("-r", "--random-dest", action = 'store_true', default = False, help = "verbosity")
21
+  parser.add_option("-v", "--verbose", action = 'store_true', default = False, help = "verbosity")
22
+
23
+  (options, args) = parser.parse_args()
24
+  print >> sys.stderr, options
25
+
26
+  mesh_source = simple_rect_grid(options.source, options.source)
27
+  tree = scipy.spatial.KDTree(mesh_source.points)
28
+
29
+  mesh_dest   = simple_rect_grid(options.dest, options.dest)
30
+
31
+  if options.random_dest:
32
+    x = []
33
+    for i in xrange(options.random_total):
34
+      rx = np.random.rand() * 2 - 1
35
+      ry = np.random.rand() * 2 - 1
36
+      x.append([rx, ry])
37
+    mesh_dest.points = np.array(x)
38
+
39
+  errors = []
40
+  for x in mesh_dest.points:
41
+    (final, exact) = baker.run_baker(x, mesh_source, tree, options.extra, options.verbose)
42
+    cur_error = np.abs(final - exact)
43
+    errors.append(cur_error)
44
+
45
+  print rms(errors)

+ 16
- 0
bin/generate_qhull.pl View File

@@ -0,0 +1,16 @@
1
+#!/usr/bin/perl
2
+
3
+use strict;
4
+use warnings;
5
+
6
+my $output_filename = $ARGV[0] or die "no args";
7
+my $tmp_grid_file   = "$output_filename.grid.txt";
8
+
9
+system("grid.py > $tmp_grid_file");
10
+
11
+
12
+system("cat $tmp_grid_file | qdelaunay  GD2 s > $output_filename.txt\n"       );
13
+system("cat $tmp_grid_file | qdelaunay GD2 s Qt > ${output_filename}_qt.txt\n");
14
+system("cat $tmp_grid_file | qdelaunay GD2 s QJ > ${output_filename}_qj.txt\n");
15
+
16
+unlink($tmp_grid_file);

BIN
data/mesh.shelve View File


+ 11
- 0
data/qhull_input.txt View File

@@ -0,0 +1,11 @@
1
+2
2
+9
3
+-1.0000  -1.0000
4
+-1.0000   0.0000
5
+-1.0000   1.0000
6
+ 0.0000  -1.0000
7
+ 0.0000   0.0000
8
+ 0.0000   1.0000
9
+ 1.0000  -1.0000
10
+ 1.0000   0.0000
11
+ 1.0000   1.0000

+ 106
- 0
lib/baker.py View File

@@ -0,0 +1,106 @@
1
+from grid import exact_func
2
+import numpy as np
3
+import sys
4
+
5
+def get_phis(X, r):
6
+  """
7
+    The get_phis function is used to get barycentric coordonites for a point on a triangle.
8
+
9
+    X -- the destination point (2D)
10
+         X = [0,0]
11
+    r -- the three points that make up the triangle (2D)
12
+         r = [[-1, -1], [0, 2], [1, -1]]
13
+
14
+    this will return [0.333, 0.333, 0.333]
15
+  """
16
+
17
+  # baker: eq 7
18
+  A = np.array([
19
+                [1,   1,   1  ],
20
+                [r[0][0], r[1][0], r[2][0]],
21
+                [r[0][1], r[1][1], r[2][1]],
22
+                ])
23
+  b = np.array([1, X[0], X[1]])
24
+  try:
25
+    phi = np.linalg.solve(A,b)
26
+  except:
27
+    print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system"
28
+    phi = np.dot(np.linalg.pinv(A), b)
29
+
30
+  return phi
31
+
32
+def qlinear(X, r, q):
33
+  """
34
+    this calculates the linear portion of q from X to r
35
+
36
+    X = destination point
37
+    r = simplex points
38
+    q = CFD quantities of interest at the simplex points
39
+  """
40
+
41
+  phis = get_phis(X, r)
42
+  qlin = sum([q_i * phi_i for q_i, phi_i in zip(q[:len(phis)], phis)])
43
+  return qlin
44
+
45
+def run_baker(X, g, tree, extra_points = 3, verbose = False):
46
+  """
47
+    This is the main function to call to get an interpolation to X from the tree
48
+
49
+    X -- the destination point (2D)
50
+         X = [0,0]
51
+
52
+    g -- the grid object
53
+
54
+    tree -- the kdtree search object (built from the g mesh)
55
+  """
56
+
57
+  (dist, indicies) = tree.query(X, 3 + extra_points)
58
+
59
+  nn = [g.points[i] for i in indicies]
60
+  nq = [g.q[i]      for i in indicies]
61
+
62
+  phi = get_phis(X, nn[:3])
63
+  qlin =  nq[0] * phi[0] + nq[1] * phi[1] + nq[2] * phi[2]
64
+
65
+  error_term = 0.0
66
+
67
+  if extra_points != 0:
68
+    B = [] # baker eq 9
69
+    w = [] # baker eq 11
70
+
71
+    for index in indicies[3:]:
72
+      (phi1,phi2,phi3) = get_phis(g.points[index], nn)
73
+      B.append([phi1 * phi2, phi2*phi3, phi3*phi1])
74
+
75
+      w.append(g.q[index] - qlinear(g.points[index], nn, nq))
76
+
77
+    B = np.array(B)
78
+    w = np.array(w)
79
+
80
+    A = np.dot(B.T, B)
81
+    b = np.dot(B.T, w)
82
+
83
+    # baker solve eq 10
84
+    try:
85
+      (a, b, c) = np.linalg.solve(A,b)
86
+    except:
87
+      print >> sys.stderr, "warning: linear calculation went bad, resorting to np.linalg.pinv"
88
+      (a, b, c) = np.dot(np.linalg.pinv(A), b)
89
+
90
+    error_term =  a * phi[0] * phi[1]\
91
+                + b * phi[1] * phi[2]\
92
+                + c * phi[2] * phi[0]
93
+
94
+  exact = exact_func(X[0], X[1])
95
+  q_final = qlin + error_term
96
+
97
+  if verbose:
98
+    print "current point : %s" % X
99
+    print "exact         : %0.4f" % exact
100
+    print "qlin          : %0.4f" % qlin
101
+    print "qlinerr       : %0.4f" % np.abs(exact - qlin)
102
+    print "q_final       : %0.4f" % q_final
103
+    print "q_final_err   : %0.4f" % np.abs(exact - q_final)
104
+    print
105
+
106
+  return (q_final, exact)

+ 77
- 0
lib/grid.py View File

@@ -0,0 +1,77 @@
1
+#!/usr/bin/python
2
+
3
+import sys
4
+import numpy as np
5
+import scipy.spatial
6
+
7
+def exact_func(x, y):
8
+  return np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2)
9
+  return np.sin(x * np.pi) * np.cos(y * np.pi)
10
+
11
+
12
+class grid(object):
13
+  def __init__(self, points, q):
14
+    self.points = np.array(points)
15
+    self.q      = np.array(q)
16
+
17
+  def __str__(self):
18
+    r = ''
19
+    assert( len(self.points) == len(self.q) )
20
+    for i in xrange(len(self.points)):
21
+      r += "%r: %0.4f\n" % ( self.points[i], self.q[i] )
22
+    return r
23
+
24
+class simple_rect_grid(grid):
25
+  def __init__(self, xres = 5, yres = 5):
26
+    xmin = -1.0
27
+    xmax =  1.0
28
+    xspan = xmax - xmin
29
+    xdel = xspan / float(xres - 1)
30
+
31
+    ymin = -1.0
32
+    ymay =  1.0
33
+    yspan = ymay - ymin
34
+    ydel = yspan / float(yres - 1)
35
+
36
+
37
+    self.points = []
38
+    self.q      = []
39
+    for x in xrange(xres):
40
+      cur_x = xmin + (x * xdel)
41
+      for y in xrange(yres):
42
+        cur_y = ymin + (y * ydel)
43
+        self.points.append([cur_x, cur_y])
44
+        self.q.append(exact_func(cur_x, cur_y))
45
+    self.points = np.array(self.points)
46
+    self.q      = np.array(self.q)
47
+
48
+
49
+  def for_qhull(self):
50
+    r = '2\n'
51
+    r += '%d\n' % len(self.points)
52
+    for p in self.points:
53
+      r += "%f %f\n" % (p[0], p[1])
54
+    return r
55
+
56
+
57
+class simple_random_grid(simple_rect_grid):
58
+  def __init__(self, num_points = 10):
59
+    self.points = []
60
+    self.q = []
61
+
62
+    r = np.random
63
+
64
+    for i in xrange(num_points):
65
+      cur_x = r.rand()
66
+      cur_y = r.rand()
67
+
68
+      self.points.append([cur_x, cur_y])
69
+      self.q.append(exact_func(cur_x, cur_y))
70
+
71
+    self.points = np.array(self.points)
72
+    self.q      = np.array(self.q)
73
+
74
+
75
+if __name__ == '__main__':
76
+  g = simple_random_grid(100)
77
+  print g.for_qhull()

+ 4
- 0
lib/smcqhull.py View File

@@ -0,0 +1,4 @@
1
+#!/usr/bin/python
2
+
3
+if __name__ == '__main__':
4
+  print "hello world"

+ 12
- 0
lib/tools.py View File

@@ -0,0 +1,12 @@
1
+import numpy as np
2
+
3
+
4
+def rms(errors):
5
+  """
6
+    root mean square calculation
7
+  """
8
+  r = 0.0
9
+  for i in errors:
10
+    r += np.power(i, 2)
11
+  r = np.sqrt(r / len(errors))
12
+  return r

+ 75
- 0
test/baker.test.py View File

@@ -0,0 +1,75 @@
1
+#!/usr/bin/python
2
+
3
+import unittest
4
+import baker
5
+
6
+import numpy as np
7
+import scipy.spatial
8
+
9
+class TestSequenceFunctions(unittest.TestCase):
10
+  def setUp(self):
11
+    self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]]
12
+    self.approx_fmt = "%0.6f"
13
+
14
+  def testGetPhis(self):
15
+
16
+    X = [0,0]
17
+    r = [[-1, -1], [0, 2], [1, -1]]
18
+
19
+    result = baker.get_phis(X, r)
20
+    result = [self.approx_fmt % i for i in result]
21
+
22
+    right_answer = [self.approx_fmt % i for i in [1/3.0, 1/3.0, 1/3.0]]
23
+
24
+    for a,b in zip(result, right_answer):
25
+      self.assertEqual(a,b)
26
+
27
+  def testGetPhis2(self):
28
+
29
+    X = [0.5,0.25]
30
+    r = [[0, 0], [1, 0], [1, 1]]
31
+
32
+    result = baker.get_phis(X, r)
33
+
34
+    right_answer = [0.5, 0.25, 0.25]
35
+
36
+    for a,b in zip(result, right_answer):
37
+      self.assertEqual(a,b)
38
+
39
+  def testQlinear(self):
40
+    X = [0.5, 0.25]
41
+    r = [[0, 0], [1, 0], [1, 1]]
42
+    q = [1, 0, 0]
43
+
44
+    result = baker.qlinear(X, r, q)
45
+
46
+    right_answer = 0.5
47
+
48
+    self.assertEqual(result, right_answer)
49
+
50
+  def testRunBaker(self):
51
+    X = [0.5, 0.25]
52
+
53
+    all_points = [
54
+                   [ 0, 0], # 0
55
+                   [ 1, 0], # 1
56
+                   [ 1, 1], # 2
57
+                   [ 0, 1], # 3
58
+                   [ 1,-1], # 4
59
+                   [ 0,-1], # 5
60
+                   [-1, 1], # 6
61
+                   [-1, 0], # 7
62
+                   [-1,-1], # 8
63
+                 ]
64
+    q = [1, 0, 0, 0, 0, 0, 0, 0, 0]
65
+
66
+    tree = scipy.spatial.KDTree(all_points)
67
+    baker.run_baker(X, 
68
+    result = 3
69
+    right_answer = 5
70
+
71
+    self.assertEqual(result, right_answer)
72
+
73
+if __name__ == '__main__':
74
+  suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
75
+  unittest.TextTestRunner(verbosity=2).run(suite)

+ 33
- 0
test/qhull.test.py View File

@@ -0,0 +1,33 @@
1
+#!/usr/bin/python
2
+
3
+import delaunay
4
+import subprocess
5
+
6
+QFILE = '/tmp/forQhull.txt'
7
+
8
+l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]]
9
+
10
+qdelauney_file = open(QFILE, 'w')
11
+
12
+qdelauney_file.write("%d\n" % len(l[0]))
13
+qdelauney_file.write("%d\n" % len(l))
14
+for i in l:
15
+  qdelauney_file.write("%0.3f %0.3f\n" % (i[0], i[1]))
16
+
17
+dt = delaunay.Triangulation(l)
18
+print  dt.indices
19
+
20
+answer = [
21
+           [4,1,3],
22
+           [1,5,0],
23
+           [5,1,4],
24
+           [7,3,6],
25
+           [7,4,3],
26
+           [7,5,4],
27
+           [5,7,8],
28
+         ]
29
+
30
+print  answer == dt.indices
31
+print dt.indices
32
+
33
+print "now run /usr/bin/qdelaunay Qt i < %s" % QFILE

+ 29
- 0
test/utest.py View File

@@ -0,0 +1,29 @@
1
+#!/usr/bin/python
2
+
3
+
4
+import random
5
+import unittest
6
+import delaunay
7
+
8
+class TestSequenceFunctions(unittest.TestCase):
9
+  def setUp(self):
10
+    self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]]
11
+
12
+
13
+  def testQhull(self):
14
+    dt = delaunay.Triangulation(self.l)
15
+    answer = [
16
+               [4,1,3],
17
+               [1,5,0],
18
+               [5,1,4],
19
+               [7,3,6],
20
+               [7,4,3],
21
+               [7,5,4],
22
+               [5,7,8],
23
+             ]
24
+
25
+    self.assertEqual(dt.indices, answer)
26
+
27
+if __name__ == '__main__':
28
+  suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
29
+  unittest.TextTestRunner(verbosity=5).run(suite)

+ 26
- 0
tools/mkpoints.py View File

@@ -0,0 +1,26 @@
1
+#!/usr/bin/python
2
+
3
+import sqlite3
4
+import sys, os
5
+import numpy as np
6
+
7
+output_file = 'points.db'
8
+
9
+if os.path.exists(output_file):
10
+  print "found old db, erasing"
11
+  os.unlink(output_file)
12
+
13
+create = 'CREATE TABLE points (point_id integer primary key, x float, y float)'
14
+
15
+sqconn = sqlite3.connect(output_file)
16
+sqc = sqconn.cursor()
17
+sqc.execute(create)
18
+
19
+for i in xrange(int(1e6)):
20
+  rx = np.random.rand() * 2 - 1
21
+  ry = np.random.rand() * 2 - 1
22
+  sqc.execute('INSERT INTO points VALUES (NULL, ?, ?)', (rx, ry))
23
+
24
+
25
+sqconn.commit()
26
+sqconn.close()

+ 22
- 0
tools/multi/forktest.py View File

@@ -0,0 +1,22 @@
1
+#!/usr/bin/python
2
+
3
+import os
4
+import time
5
+
6
+def work():
7
+  j = 0.0
8
+  for i in xrange(10000000):
9
+    j = j + j/2.0
10
+
11
+
12
+
13
+pid = os.fork()
14
+if pid != 0:
15
+  print "i spawned %d and am working" % pid
16
+  work()
17
+else:
18
+  print "i am spawn, and am working"
19
+  work()
20
+  os._exit(0)
21
+
22
+print "parent done"

+ 12
- 0
tools/multi/pool.py View File

@@ -0,0 +1,12 @@
1
+#!/usr/bin/python
2
+
3
+from multiprocessing import Pool
4
+
5
+def f(x):
6
+  return x*x
7
+
8
+if __name__ == '__main__':
9
+  pool = Pool(processes = 4)
10
+  result = pool.apply_async(f, (10,))
11
+  print result.get()
12
+  print pool.map(f, range(10))

+ 24
- 0
tools/multi/smp.py View File

@@ -0,0 +1,24 @@
1
+#!/usr/bin/python
2
+
3
+from multiprocessing import Process, Lock
4
+
5
+def f(l, i, d):
6
+  d['a'] = 'shutup'
7
+  print 'hello world', i, d
8
+  j = 0.0
9
+  for i in xrange(1000000):
10
+    j = j + j/2.0
11
+
12
+if __name__ == '__main__':
13
+  lock = Lock()
14
+
15
+  d = {'a': 1, 'b': 2}
16
+
17
+  ps = []
18
+  for num in range(2):
19
+    ps.append(Process(target=f, args=(lock, num, d)))
20
+
21
+  for p in ps: p.start()
22
+  for p in ps: p.join()
23
+
24
+  print d

+ 22
- 0
tools/multi/threadtest.py View File

@@ -0,0 +1,22 @@
1
+#!/usr/bin/python
2
+
3
+import thread, time
4
+
5
+def work():
6
+  j = 0.0
7
+  for i in xrange(10000000):
8
+    j = j + j/2.0
9
+
10
+def child(tid):
11
+  print "%d start" % tid
12
+  work()
13
+  print "%d  stop" % tid
14
+
15
+def parent():
16
+  i = 0
17
+  while True:
18
+    i += 1
19
+    thread.start_new(child, (i,))
20
+    if raw_input() == 'q': break
21
+
22
+parent()

Loading…
Cancel
Save