From 5b40dea605ebf3f845944ffe22485abe2882b98d Mon Sep 17 00:00:00 2001
From: Michael Krayer <michael.krayer@posteo.de>
Date: Wed, 4 Aug 2021 00:58:01 +0200
Subject: [PATCH] implementing particle masking, unfinished

---
 particle.py | 47 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 46 insertions(+), 1 deletion(-)

diff --git a/particle.py b/particle.py
index 7c842b3..a7d9484 100644
--- a/particle.py
+++ b/particle.py
@@ -123,6 +123,52 @@ 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
+        posd = [pos.copy()]
+        for axis in range(3):
+            if self.period[axis] is not None:
+                nposd = len(posd)
+                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])
+                        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])
+                        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):
+        # 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
         position = np.vstack([self.attr[key] for key in ('x','y','z')]).transpose()
@@ -130,7 +176,6 @@ class Particles:
         for key in self.attr:
             mesh[key] = self.attr[key]
         return mesh
-
     def glyph(self,theta_resolution=30,phi_resolution=30,deep=False):
         import pyvista as pv
         assert self.has_attribute('r'), "Attribute 'r' required."