Changeset 2642 for trunk

Timestamp:

07/08/2008 04:44:52 PM (6 months ago)

Author:

wd15

Message:

enabled interpolation of a CellVariable? to multiple points, added special methods for locating nearest cell for UniformGrids?, added first order interpolation

Location:

trunk/fipy

Files:

• meshes/common/mesh.py (modified) (1 diff)
• meshes/numMesh/uniformGrid1D.py (modified) (1 diff)
• meshes/numMesh/uniformGrid2D.py (modified) (1 diff)
• meshes/numMesh/uniformGrid3D.py (modified) (1 diff)
• variables/cellVariable.py (modified) (2 diffs)

trunk/fipy/meshes/common/mesh.py

@@ -611 +611 @@
         return self._getCellsByID([self._getNearestCellID(point)])[0]
 
-    def _getNearestCellID(self, point):
+    def _getNearestCellID(self, points):
+        """
+        Test cases
+
+           >>> from fipy import *
+           >>> m0 = Grid2D(dx=(.1, 1., 10.), dy=(.1, 1., 10.))
+           >>> m1 = Grid2D(nx=2, ny=2, dx=5., dy=5.)
+           >>> print m0._getNearestCellID(m1.getCellCenters())
+           [4 5 7 8]
+           
+        """
+        points = numerix.resize(points, (self.getNumberOfCells(), len(points), len(points[0]))).swapaxes(0,1)
+
         try:
-            tmp = self.getCellCenters() - point
+            tmp = self.getCellCenters()[...,numerix.newaxis] - points
         except TypeError:
-            tmp = self.getCellCenters() - PhysicalField(point)
-        i = numerix.argmin(numerix.add.reduce((tmp * tmp), axis = 0))
-        return i    
+            tmp = self.getCellCenters()[...,numerix.newaxis] - PhysicalField(points)
+
+        return numerix.argmin(numerix.dot(tmp, tmp, axis = 0), axis=0)
 
 ## pickling

trunk/fipy/meshes/numMesh/uniformGrid1D.py

@@ -259 +259 @@
         pass
 
+    def _getNearestCellID(self, points):
+        """
+        Test cases
+
+           >>> from fipy import *
+           >>> m = Grid1D(nx=3)
+           >>> print m._getNearestCellID(([0., .9, 3.],))
+           [0 0 2]
+           >>> print m._getNearestCellID(([1.1],))
+           [1]
+           >>> m0 = Grid1D(nx=2, dx=1.)
+           >>> m1 = Grid1D(nx=4, dx=.5)
+           >>> print m0._getNearestCellID(m1.getCellCenters())
+           [0 0 1 1]
+           
+        """
+        x0, = self.getCellCenters()[...,0]        
+        xi, = points
+        nx, = self.getShape()
+        dx = self.dx
+        
+        i = numerix.array(numerix.rint(((xi - x0) / dx)), 'l')
+        i[i < 0] = 0
+        i[i > nx - 1] = nx - 1
+
+        return i
+
     def _test(self):
         """

trunk/fipy/meshes/numMesh/uniformGrid2D.py

@@ -531 +531 @@
         pass
 
+    def _getNearestCellID(self, points):
+        """
+        Test cases
+
+           >>> from fipy import *
+           >>> m = Grid2D(nx=3, ny=2)
+           >>> eps = numerix.array([[1e-5, 1e-5]])
+           >>> print m._getNearestCellID(((0., .9, 3.), (0., 2., 2.)))
+           [0 3 5]
+           >>> print m._getNearestCellID(([1.1], [1.5]))
+           [4]
+           >>> m0 = Grid2D(nx=2, ny=2, dx=1., dy=1.)
+           >>> m1 = Grid2D(nx=4, ny=4, dx=.5, dy=.5)
+           >>> print m0._getNearestCellID(m1.getCellCenters())
+           [0 0 1 1 0 0 1 1 2 2 3 3 2 2 3 3]
+           
+        """
+        x0, y0 = self.getCellCenters()[...,0]        
+        xi, yi = points
+        nx, ny = self.getShape()
+        dx, dy = self.dx, self.dy
+        
+        i = numerix.array(numerix.rint(((xi - x0) / dx)), 'l')
+        i[i < 0] = 0
+        i[i > nx - 1] = nx - 1
+
+        j = numerix.array(numerix.rint(((yi - y0) / dy)), 'l')
+        j[j < 0] = 0
+        j[j > ny - 1]  = ny - 1
+
+        return j * nx + i
+
+        
     def _test(self):
         """

trunk/fipy/meshes/numMesh/uniformGrid3D.py

@@ -469 +469 @@
     def _calcScaledGeometry(self):
         pass
+    
+    def _getNearestCellID(self, points):
+        x0, y0, z0 = self.getCellCenters()[...,0]        
+        xi, yi, zi = points
+        nx, ny, nz = self.getShape()
+        dx, dy, dz = self.dx, self.dy, self.dz
+        
+        i = numerix.array(numerix.rint(((xi - x0) / dx)), 'l')
+        i[i < 0] = 0
+        i[i > nx - 1] = nx - 1
+
+        j = numerix.array(numerix.rint(((yi - y0) / dy)), 'l')
+        j[j < 0] = 0
+        j[j > ny - 1]  = ny - 1
+
+        k = numerix.array(numerix.rint(((zi - z0) / dz)), 'l')
+        k[k < 0] = 0
+        k[k > nz - 1]  = nz - 1
+        
+        return k * ny * nx + j * nx + i
         
     def _test(self):

trunk/fipy/variables/cellVariable.py

@@ -139 +139 @@
             hasOld = 0)
             
-    def __call__(self, point=None):
-        if point is not None:
-            return self[...,self.getMesh()._getNearestCellID(point)]
+    def __call__(self, points=None, order=0):
+        r"""
+        Interpolates the CellVariable to a set of points using a
+        method that has a memory requirement on the order of Ncells by
+        Npoints in general, but uses only Ncells when the
+        CellVariable's mesh is a UniformGrid object.
+
+        :Parameters:
+
+           - `points`: A point or set of points in the format (X, Y, Z)
+           - `order`: The order of interpolation, 0 or 1, default is 0
+
+        Tests
+
+            >>> from fipy import *
+            >>> m = Grid2D(nx=3, ny=2)
+            >>> v = CellVariable(mesh=m, value=m.getCellCenters()[0])
+            >>> print v(((0., 1.1, 1.2), (0., 1., 1.)))
+            [ 0.5  1.5  1.5]
+            >>> print v(((0., 1.1, 1.2), (0., 1., 1.)), order=1)
+            [ 0.25  1.1   1.2 ]
+            >>> m0 = Grid2D(nx=2, ny=2, dx=1., dy=1.)
+            >>> m1 = Grid2D(nx=4, ny=4, dx=.5, dy=.5)
+            >>> x, y = m0.getCellCenters()
+            >>> v0 = CellVariable(mesh=m0, value=x * y)
+            >>> print v0(m1.getCellCenters())
+            [ 0.25  0.25  0.75  0.75  0.25  0.25  0.75  0.75  0.75  0.75  2.25  2.25
+              0.75  0.75  2.25  2.25]
+            >>> print v0(m1.getCellCenters(), order=1)
+            [ 0.125  0.25   0.5    0.625  0.25   0.375  0.875  1.     0.5    0.875
+              1.875  2.25   0.625  1.     2.25   2.625]
+
+        """           
+        if points is not None:
+
+            if order == 0:
+                return self[...,self.getMesh()._getNearestCellID(points)]
+
+            elif order == 1:
+                cellID = self.getMesh()._getNearestCellID(points)
+                return self[...,self.getMesh()._getNearestCellID(points)] + numerix.dot(points - self.getMesh().getCellCenters()[...,cellID], self.getGrad()[...,cellID])
+
+            else:
+                raise ValueError, 'order should be either 0 or 1'
+
         else:
             return _MeshVariable.__call__(self)
-
+        
     def getCellVolumeAverage(self):
         r"""
@@ -476 +518 @@
             self.old.setValue(dict['old'].getValue())
 
-
 class _ReMeshedCellVariable(CellVariable):
     def __init__(self, oldVar, newMesh):