[yt-svn] commit/yt: Andrew Myers: porting streamlines over after the volume rendering refactor and fixing a bug in the star particle rendering

Bitbucket commits-noreply at bitbucket.org
Tue Jul 3 04:42:19 PDT 2012


1 new commit in yt:


https://bitbucket.org/yt_analysis/yt/changeset/dc78a9b9bb47/
changeset:   dc78a9b9bb47
branch:      yt
user:        Andrew Myers
date:        2012-07-01 01:26:41
summary:     porting streamlines over after the volume rendering refactor and fixing a bug in the star particle rendering
affected #:  2 files

diff -r a7255e4514626efb26996347df3347caddbe36b0 -r dc78a9b9bb475f740b4417ac1440786697129444 yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -78,8 +78,12 @@
     cdef public object my_data
     cdef public object LeftEdge
     cdef public object RightEdge
-    cdef int parent_grid_id
+    cdef public int parent_grid_id
     cdef VolumeContainer *container
+    cdef kdtree_utils.kdtree *star_list
+    cdef np.float64_t star_er
+    cdef np.float64_t star_sigma_num
+    cdef np.float64_t star_coeff
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -88,7 +92,8 @@
                   int parent_grid_id, data,
                   np.ndarray[np.float64_t, ndim=1] left_edge,
                   np.ndarray[np.float64_t, ndim=1] right_edge,
-                  np.ndarray[np.int64_t, ndim=1] dims):
+                  np.ndarray[np.int64_t, ndim=1] dims,
+		  star_kdtree_container star_tree = None):
         # The data is likely brought in via a slice, so we copy it
         cdef np.ndarray[np.float64_t, ndim=3] tdata
         self.container = NULL
@@ -111,6 +116,16 @@
         for i in range(n_fields):
             tdata = data[i]
             c.data[i] = <np.float64_t *> tdata.data
+        if star_tree is None:
+            self.star_list = NULL
+        else:
+            self.set_star_tree(star_tree)
+
+    def set_star_tree(self, star_kdtree_container star_tree):
+        self.star_list = star_tree.tree
+        self.star_sigma_num = 2.0*star_tree.sigma**2.0
+        self.star_er = 2.326 * star_tree.sigma
+        self.star_coeff = star_tree.coeff
 
     def __dealloc__(self):
         # The data fields are not owned by the container, they are owned by us!
@@ -119,6 +134,90 @@
         if self.container.data != NULL: free(self.container.data)
         free(self.container)
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    def integrate_streamline(self, pos, np.float64_t h, mag):
+        cdef np.float64_t cmag[1]
+        cdef np.float64_t k1[3], k2[3], k3[3], k4[3]
+        cdef np.float64_t newpos[3], oldpos[3]
+        for i in range(3):
+            newpos[i] = oldpos[i] = pos[i]
+        self.get_vector_field(newpos, k1, cmag)
+        for i in range(3):
+            newpos[i] = oldpos[i] + 0.5*k1[i]*h
+
+        if not (self.LeftEdge[0] < newpos[0] and newpos[0] < self.RightEdge[0] and \
+                self.LeftEdge[1] < newpos[1] and newpos[1] < self.RightEdge[1] and \
+                self.LeftEdge[2] < newpos[2] and newpos[2] < self.RightEdge[2]):
+            if mag is not None:
+                mag[0] = cmag[0]
+            for i in range(3):
+                pos[i] = newpos[i]
+            return
+
+        self.get_vector_field(newpos, k2, cmag)
+        for i in range(3):
+            newpos[i] = oldpos[i] + 0.5*k2[i]*h
+
+        if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
+                self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
+                self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
+            if mag is not None:
+                mag[0] = cmag[0]
+            for i in range(3):
+                pos[i] = newpos[i]
+            return
+
+        self.get_vector_field(newpos, k3, cmag)
+        for i in range(3):
+            newpos[i] = oldpos[i] + k3[i]*h
+
+        if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
+                self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
+                self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
+            if mag is not None:
+                mag[0] = cmag[0]
+            for i in range(3):
+                pos[i] = newpos[i]
+            return
+
+        self.get_vector_field(newpos, k4, cmag)
+
+        for i in range(3):
+            pos[i] = oldpos[i] + h*(k1[i]/6.0 + k2[i]/3.0 + k3[i]/3.0 + k4[i]/6.0)
+
+        if mag is not None:
+            for i in range(3):
+                newpos[i] = pos[i]
+            self.get_vector_field(newpos, k4, cmag)
+            mag[0] = cmag[0]
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef void get_vector_field(self, np.float64_t pos[3],
+                               np.float64_t *vel, np.float64_t *vel_mag):
+        cdef np.float64_t dp[3]
+        cdef int ci[3]
+        cdef VolumeContainer *c = self.container # convenience
+
+        for i in range(3):
+            ci[i] = (int)((pos[i]-self.LeftEdge[i])/c.dds[i])
+            dp[i] = (pos[i] - self.LeftEdge[i])%(c.dds[i])
+
+        cdef int offset = ci[0] * (c.dims[1] + 1) * (c.dims[2] + 1) \
+                          + ci[1] * (c.dims[2] + 1) + ci[2]
+
+        vel_mag[0] = 0.0
+        for i in range(3):
+            vel[i] = offset_interpolate(c.dims, dp, c.data[i] + offset)
+            vel_mag[0] += vel[i]*vel[i]
+        vel_mag[0] = np.sqrt(vel_mag[0])
+        if vel_mag[0] != 0.0:
+            for i in range(3):
+                vel[i] /= vel_mag[0]
+
 cdef struct ImageContainer:
     np.float64_t *vp_pos, *vp_dir, *center, *image,
     np.float64_t pdx, pdy, bounds[4]


diff -r a7255e4514626efb26996347df3347caddbe36b0 -r dc78a9b9bb475f740b4417ac1440786697129444 yt/visualization/volume_rendering/grid_partitioner.py
--- a/yt/visualization/volume_rendering/grid_partitioner.py
+++ b/yt/visualization/volume_rendering/grid_partitioner.py
@@ -97,7 +97,7 @@
             dd = [d[sl[0][0]:sl[0][1]+1,
                     sl[1][0]:sl[1][1]+1,
                     sl[2][0]:sl[2][1]+1].copy() for d in vcds]
-            pgs.append(PartitionedGrid(grid.id, len(self.fields), dd,
+            pgs.append(PartitionedGrid(grid.id, dd,
                         P.LeftEdge, P.RightEdge, sl[-1]))
         return pgs

Repository URL: https://bitbucket.org/yt_analysis/yt/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.



More information about the yt-svn mailing list