Commits

Matthew Turk  committed 40fa770

Adding a bunch of neighbor-searching to the SPH/Particle/Octree code.

  • Participants
  • Parent commits 7d64152

Comments (0)

Files changed (1)

 dle = na.zeros(3)
 dre = na.ones(3) * f["/Header"].attrs["BoxSize"]
 poc = ParticleOctreeContainer([16,16,16], dle, dre)
+poc.allocate_root()
 lc = None
 do_plot = False
 
 for i in range(8):
-    f = h5py.File("snapshot_033/snap_033.%s.hdf5" % i)
+    f = h5py.File("snapshot_033/snap_033.%s.hdf5" % (i % 8))
     pos = f["/PartType0/Coordinates"][:].astype("float64")
     poc.add(pos, i)
     nlc = poc.recursively_count()
     if lc is not None:
         for l, v in sorted(nlc.items()):
             if l not in lc: continue
-            assert(v > lc[l])
+            assert(v >= lc[l])
     lc = nlc
     if do_plot:
         pylab.clf()
 
 positions = [
     (25.0, (12.5, 12.5, 12.5)),
-    (2.5, (12.5, 12.5, 12.5)),
-    (5.0, (20.0, 20.0, 20.0)),
-    (1.0, (5.0, 5.0, 5.0)),
-    (3.0, (5.0, 5.0, 5.0)),
+    #(2.5, (12.5, 12.5, 12.5)),
+    #(5.0, (20.0, 20.0, 20.0)),
+    #(1.0, (5.0, 5.0, 5.0)),
+    #(3.0, (5.0, 5.0, 5.0)),
 ]
 
-for center, radius in positions:
-    sp = SphereSelector(FakeSphere(center, radius))
+for rad, cen in positions:
+    sp = SphereSelector(FakeSphere(rad, cen))
     mask = sp.select_octs(poc)
     domains = poc.domain_identify(mask)
-    print center, radius, mask.sum(), domains
+    print rad, cen, mask.sum(), domains
+    icoords = poc.icoords(mask, mask.sum())
+    ires = poc.ires(mask, mask.sum())
+    mask[:] = True
+    ic = poc.icoords(mask, mask.size)
+    ir = poc.ires(mask, mask.size)
+    rl = ic[ir == 0,:]
+    fc = poc.fcoords(mask, mask.size)
+    max_point = fc[na.argmax(ir), :]
+    print ic.max(), fc.max(), ir.max()
 
+tnp = poc.count_neighbor_particles(max_point)
+print "Got %s particles" % (tnp)
+N = 30000
+ps = na.random.random((N,3)) * (dre - dle) + dle
+t1 = time.time()
+for i in range(N):
+    poc.count_neighbor_particles((ps[i,0], ps[i,1], ps[i,2]))
+t2 = time.time()
+print "Took %0.3e" % ((t2-t1)/N)
+
+N = 20
+ps = na.random.random((N,3)) * (dre - dle) + dle
+for i in range(N):
+    bounds = poc.get_neighbor_boundaries((ps[i,0], ps[i,1], ps[i,2]))
+    pylab.clf()
+    pylab.plot(ps[i,0], ps[i,1], 'x')
+    for ri in range(27):
+        print bounds[ri,:]
+        r = pylab.Rectangle((bounds[ri, 0], bounds[ri, 1]),
+                             bounds[ri, 3], bounds[ri, 4],
+                             fill = False)
+        pylab.gca().add_patch(r)
+    pylab.xlim(dle[0], dre[0])
+    pylab.ylim(dle[1], dre[1])
+    pylab.savefig("neighbors_%03i.png" % i)