mwtkrayer
/
suspendtools
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added field masking with particle data

This commit is contained in:
Michael Krayer 2021-08-05 00:21:22 +02:00
parent 9527941936
commit 5bf6b6e102
1 changed files with 22 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
particle.py
View File

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