field.py

@@ -1,13 +1,9 @@
 import numpy as np
 class Field3d:
-    def __init__(self,data,origin,spacing,deep=False):
+    def __init__(self,data,origin,spacing):
         assert len(origin)==3,  "'origin' must be of length 3"
         assert len(spacing)==3, "'spacing' must be of length 3"
-        assert isinstance(data,np.ndarray), "'data' must be numpy.ndarray."
-        if deep:
-            self.data = data.copy()
-        else:
-            self.data = data
+        self.data      = np.array(data)
         self.origin    = tuple([float(x) for x in origin])
         self.spacing   = tuple([float(x) for x in spacing])
         self.eps_collapse = 1e-8
@@ -191,28 +187,16 @@ 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),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."
-        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))
+    def extract_subfield(self,idx_origin,dim,stride=(1,1,1)):
+        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))
         origin  = self.coordinate(idx_origin)
         spacing = tuple(self.spacing[ii]*stride[ii] for ii in range(3))
-        data = self.data[sl]
-        return Field3d(data,origin,spacing,deep=deep)
+        data = self.data[sl].copy()
+        return Field3d(data,origin,spacing)
 
     def coordinate(self,idx,axis=None):
         if axis is None:

parallel.py

@@ -728,8 +728,7 @@ 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),
-                                            deep=True):
+                                            no_upper_ghost=(False,False,False)):
         '''Returns the field with a stride applied.'''
         stride    = np.array(stride,dtype=int)
         num_ghost = np.array(num_ghost,dtype=int)
@@ -766,7 +765,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,deep=deep)
+        return self.field[key].extract_subfield(idx_origin,num_points,stride=stride)
 
     def _merge_at_root(self,key,stride=(1,1,1)):
         '''Returns the entire field gathered from all processors with a 

particle.py

@@ -123,12 +123,9 @@ class Particles:
                 key = ('x','y','z')[axis]
                 self.attr[key] %= self.period[axis]
         return
-    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
+    def position_duplicates(self,ipart,padding=0.0):
+        pos  = np.array((px[ipart],py[ipart],pz[ipart]))
+        rp   = pr[ipart]+padding
         posd = [pos.copy()]
         for axis in range(3):
             if self.period[axis] is not None:
@@ -136,29 +133,41 @@ class Particles:
                 if pos[axis]-rp<0.0:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = tmp[axis]+self.period[axis]
+                        tmp[axis] = np.mod(tmp[axis]-rp,self.period[axis])
                         posd.append(tmp)
                 if pos[axis]+rp>self.period[axis]:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = tmp[axis]-self.period[axis]
+                        tmp[axis] = np.mod(tmp[axis]+rp,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.'''
-        assert self.has_attribute('r'), "Attribute 'r' required."
-        for ipart in range(self.num):
-            # Slice a box from the field around the particle
-            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
+      '''Fills grid points which lie inside of solid phase with values.'''
+      for ipart in range(0,part.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
         return
     def to_vtk(self,deep=False):
         import pyvista as pv