added field masking with particle data

field.py

@@ -1,9 +1,13 @@
 import numpy as np
 class Field3d:
-    def __init__(self,data,origin,spacing):
+    def __init__(self,data,origin,spacing,deep=False):
         assert len(origin)==3,  "'origin' must be of length 3"
         assert len(spacing)==3, "'spacing' must be of length 3"
-        self.data      = np.array(data)
+        assert isinstance(data,np.ndarray), "'data' must be numpy.ndarray."
+        if deep:
+            self.data = data.copy()
+        else:
+            self.data = data
         self.origin    = tuple([float(x) for x in origin])
         self.spacing   = tuple([float(x) for x in spacing])
         self.eps_collapse = 1e-8
@@ -187,16 +191,28 @@ class Field3d:
         self.data[ib:ib+nx,jb:jb+ny,kb:kb+nz] = subfield.data[:,:,:]
         return
 
-    def extract_subfield(self,idx_origin,dim,stride=(1,1,1)):
+    def extract_subfield(self,idx_origin,dim,stride=(1,1,1),deep=False,strict_bounds=True):
+        if strict_bounds:
             assert all(idx_origin[ii]>=0 and idx_origin[ii]<self.dim(axis=ii) for ii in range(3)),\
                     "'origin' is out-of-bounds."
             assert all(idx_origin[ii]+stride[ii]*(dim[ii]-1)<self.dim(axis=ii) for ii in range(3)),\
                     "endpoint is out-of-bounds."
-        sl = tuple(slice(idx_origin[ii],idx_origin[ii]+stride[ii]*dim[ii],stride[ii]) for ii in range(3))
+        else:
+            for ii in range(3):
+                if idx_origin[ii]<0:
+                    dim[ii] += idx_origin[ii]
+                    idx_origin[ii] = 0
+                if idx_origin[ii]+stride[ii]*(dim[ii]-1)>=self.dim(axis=ii):
+                    dim[ii] = (self.dim(axis=ii)-idx_origin[ii])//stride[ii]
+                if dim[ii]<=0:
+                    return None
+        sl = tuple(slice(idx_origin[ii],
+                         idx_origin[ii]+stride[ii]*dim[ii],
+                         stride[ii]) for ii in range(3))
         origin  = self.coordinate(idx_origin)
         spacing = tuple(self.spacing[ii]*stride[ii] for ii in range(3))
-        data = self.data[sl].copy()
-        return Field3d(data,origin,spacing)
+        data = self.data[sl]
+        return Field3d(data,origin,spacing,deep=deep)
 
     def coordinate(self,idx,axis=None):
         if axis is None:

parallel.py

@@ -728,7 +728,8 @@ class PPP:
         return self._subfield(key,stride,(0,0,0))
 
     def _subfield(self,key,stride,num_ghost,no_lower_ghost=(False,False,False),
-                                            no_upper_ghost=(False,False,False)):
+                                            no_upper_ghost=(False,False,False),
+                                            deep=True):
         '''Returns the field with a stride applied.'''
         stride    = np.array(stride,dtype=int)
         num_ghost = np.array(num_ghost,dtype=int)
@@ -765,7 +766,7 @@ class PPP:
         assert all(idx_origin>=0)
         assert all(idx_endpoint<np.array(self.field[key].dim()))
         # Return the subfield
-        return self.field[key].extract_subfield(idx_origin,num_points,stride=stride)
+        return self.field[key].extract_subfield(idx_origin,num_points,stride=stride,deep=deep)
 
     def _merge_at_root(self,key,stride=(1,1,1)):
         '''Returns the entire field gathered from all processors with a 

particle.py

@@ -123,9 +123,12 @@ class Particles:
                 key = ('x','y','z')[axis]
                 self.attr[key] %= self.period[axis]
         return
-    def position_duplicates(self,ipart,padding=0.0):
-        pos  = np.array((px[ipart],py[ipart],pz[ipart]))
-        rp   = pr[ipart]+padding
+    def position_with_duplicates(self,ipart,padding=0.0):
+        pos  = np.array(
+            (self.attr['x'][ipart],
+             self.attr['y'][ipart],
+             self.attr['z'][ipart]))
+        rp = self.attr['r'][ipart]+padding
         posd = [pos.copy()]
         for axis in range(3):
             if self.period[axis] is not None:
@@ -133,41 +136,29 @@ class Particles:
                 if pos[axis]-rp<0.0:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]-rp,self.period[axis])
+                        tmp[axis] = tmp[axis]+self.period[axis]
                         posd.append(tmp)
                 if pos[axis]+rp>self.period[axis]:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]+rp,self.period[axis])
+                        tmp[axis] = tmp[axis]-self.period[axis]
                         posd.append(tmp)
         return posd
     def mask_field(self,fld,cval=np.nan,padding=0.0):
         '''Fills grid points which lie inside of solid phase with values.'''
-      for ipart in range(0,part.num):
+        assert self.has_attribute('r'), "Attribute 'r' required."
+        for ipart in range(self.num):
             # Slice a box from the field around the particle
-        xp,yp,zp,rp = (px[ipart],py[ipart],pz[ipart],pr[ipart])
-        rp += padding
-        idxlo = np.array(fld.nearest_gridpoint(np.array(xp,yp,zp)-rp,lower=True))
-        idxhi = idxlo+2*rp/fld.spacing
-        
-
-        # Get bounding box of particle
-        idx_x = np.nonzero((xg>=xp-rp) & (xg<=xp+rp))[0]
-        idx_y = np.nonzero((yg>=yp-rp) & (yg<=yp+rp))[0]
-        idx_z = np.nonzero((zg>=zp-rp) & (zg<=zp+rp))[0]
-        # Triple for loop
-        for ii in range(idx_x[0],idx_x[-1]+1):
-          Dx = xg[ii]-xp
-          for jj in range(idx_y[0],idx_y[-1]+1):
-            Dy = yg[jj]-yp
-            for kk in range(idx_z[0],idx_z[-1]+1):
-              Dz = zg[kk]-zp
-              isInside = Dx*Dx+Dy*Dy+Dz*Dz <= rp*rp
-              if isInside:
-                if reconstruct:
-                  self.field[key][ii,jj,kk] = coeff_lin + coeff_rotx*Dx + coeff_roty*Dy + coeff_rotz*Dz
-                else:
-                  self.field[key][ii,jj,kk] = cval
+            rp = self.attr['r'][ipart]+padding
+            for pos in self.position_with_duplicates(ipart,padding=padding):
+                idxlo = np.array(fld.nearest_gridpoint(pos-rp,lower=True))
+                sfdim = np.ceil(2*rp/fld.spacing+1).astype('int')
+                sf = fld.extract_subfield(idxlo,sfdim,deep=False,strict_bounds=False)
+                xsf = sf.x().reshape((-1,1,1))
+                ysf = sf.y().reshape((1,-1,1))
+                zsf = sf.z().reshape((1,1,-1))
+                dist = (xsf-pos[0])**2 + (ysf-pos[1])**2 + (zsf-pos[2])**2
+                sf.data[dist<=rp*rp] = cval
         return
     def to_vtk(self,deep=False):
         import pyvista as pv