DISPATCH
_planes.py
1 # -*- coding: utf-8 -*-
2 """
3 Created on Sun Jul 29 20:35:19 2018
4 
5 @author: Aake
6 """
7 # Pythn 2/3 compatibility
8 from __future__ import print_function
9 
10 import numpy as np
11 
12 def contribution(p,dir=2,point=0.5,iv=0,i4=0,all=False,verbose=0):
13  """ Return the contribution from a patch to a unitgrid_plane """
14  f=0.0
15  v=0.0
16  ds=p.ds[dir]
17  # Bounds of interval along the dir axes
18  lc=p.position[dir]-p.size[dir]/2.0
19  uc=p.position[dir]+p.size[dir]/2.0
20  if all and p.guard_zones:
21  lc=lc-ds*p.ng[dir]
22  uc=uc+ds*p.ng[dir]
23  n=p.gn[dir]
24  else:
25  n=p.n[dir]
26  ds2 = ds*0.5
27  # For patches that contain the given point along the axis
28  if uc+ds2 >= point and lc-ds2 <= point:
29  if p.no_mans_land:
30  lc=lc+ds2
31  # Floating point index => integer index + floating remainder
32  fi=(point-lc)/ds
33  i=int(max(0,min(n-1,np.floor(fi))))
34  fi=fi-i
35  if verbose>2:
36  print(p.id,i,fi)
37  # The 1st no-mans-land half interval
38  if fi < 0:
39  w0=1.0+fi
40  w1=None
41  # The 2nd no-mans-land half interval
42  elif i == n-1:
43  w0=1.0-fi
44  w1=None
45  # Internal interval
46  else:
47  w0=1.0-fi
48  w1=fi
49  if verbose>1:
50  if w1:
51  print ('id:{:3d} i:{:3d} w:{:5.2f}{:5.2f}'.format(p.id,i,w0,w1))
52  else:
53  print ('id:{:3d} i:{:3d} w:{:5.2f}'.format(p.id,i,w0))
54  n+=1
55  dv=np.product(p.ds)
56  data=p.var(iv,i4=i4,all=all)
57  v=v+dv*w0
58  if w1:
59  v=v+dv*w1
60  if dir==0:
61  f=f+dv*data[i,:,:]*w0
62  if w1:
63  f=f+dv*data[i+1,:,:]*w1
64  elif dir==1:
65  f=f+dv*data[:,i,:]*w0
66  if w1:
67  f=f+dv*data[:,i+1,:]*w1
68  else:
69  f=f+dv*data[:,:,i]*w0
70  if w1:
71  f=f+dv*data[:,:,i+1]*w1
72  if all and p.guard_zones:
73  m=list(p.ng)
74  m.pop(dir)
75  f=f[m[0]:-m[0],m[1]:-m[1]]
76  return f/v,i
77  else:
78  return None,0
79 
80 def corner_indices(sn,p,dir=-1):
81  """ Return the corder indices of patch p in snapshot sn in direction dir """
82  i=(p.position-p.size/2-sn.cartesian.origin)/p.ds
83  i=[int(k+0.5) for k in i]
84  n=list(p.n)
85  if dir>=0 and dir<3:
86  i.pop(dir)
87  n.pop(dir)
88  return i[0],i[0]+n[0],i[1],i[1]+n[1]
89  else:
90  return i[0],i[0]+n[0],i[1],i[1]+n[1],i[2],i[2]+n[2]
91 
92 def unigrid_plane(snap,iv=0,x=None,y=None,z=None,point=0.5,dir=2,all=False,i4=0,verbose=0):
93  pp=snap.patches
94  p=pp[0]
95  dim=list(snap.cartesian.dims*p.n)
96  dim.pop(dir)
97  ff=np.zeros(dim)
98  if x:
99  dir=0
100  point=x
101  elif y:
102  dir=1
103  point=y
104  elif z:
105  dir=2
106  point=z
107  # List of patches that contain the plane
108  if all:
109  ds=p.ds[dir]*(2*p.ng[dir]+1)
110  else:
111  ds=p.ds[dir]
112  pp=[p for p in pp if abs(point-p.position[dir]) <= (p.size[dir]+ds)/2]
113  if verbose>0:
114  print('number of patches:',len(pp))
115  n=0
116  for p in pp:
117  n+=1
118  if verbose>1:
119  print(n)
120  elif verbose>0:
121  if n%100==0:
122  print(n)
123  k=corner_indices(snap,p,dir)
124  f,i=contribution(p,dir=dir,point=point,iv=iv,i4=i4,all=all,
125  verbose=verbose)
126  if f is not None:
127  ff[k[0]:k[1],k[2]:k[3]]+=f
128  else:
129  if verbose>1:
130  print('patch:',p.id,'range:',k,'shape: None')
131  return ff
132 
133 def unigrid_volume(snap,iv=0,i4=0,verbose=0):
134  """
135  assemble a volume from a unigrid experiment, where all
136  patches have the same size and resolution
137  """
138  pp=snap.patches
139  p=pp[0]
140  dim=tuple(snap.cartesian.dims*p.n)
141  ff=np.zeros(dim)
142  n=0
143  for p in pp:
144  n+=1
145  if verbose>1:
146  print(n)
147  elif verbose>0:
148  if n%100==0:
149  print(n)
150  k=corner_indices(snap,p)
151  ff[k[0]:k[1],k[2]:k[3],k[4]:k[5]]=p.var(iv)
152  return ff