Changeset 674 for codes/icosagcm/devel
- Timestamp:
- 01/31/18 15:14:35 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 2 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r672 r674 5 5 import matplotlib.tri as tri 6 6 import matplotlib.pyplot as plt 7 from unstructured import init_mesh, list_stencil 7 from unstructured import init_mesh, list_stencil, ker 8 8 9 9 radian=180/math.pi … … 174 174 return ne 175 175 176 def plot(tri,data ):176 def plot(tri,data, **kwargs): 177 177 plt.figure(figsize=(12,4)) 178 178 plt.gca().set_aspect('equal') 179 plt.tricontourf(tri, data, 20 )179 plt.tricontourf(tri, data, 20, **kwargs) 180 180 plt.colorbar() 181 181 plt.ylim((-90,90)) … … 210 210 211 211 # get indices for stencils 212 # primal -> vertices (unused)213 primal_vertex, dual_vertex = getvars('verticesOnCell','cellsOnVertex')214 212 # primal <-> edges 215 213 primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge') … … 221 219 dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down) 222 220 # primal <-> dual, edges <-> edges (TRiSK) 223 dual_vertex, trisk = getvars('cellsOnVertex','edgesOnEdge')221 primal_vertex, dual_vertex, trisk = getvars('verticesOnCell','cellsOnVertex','edgesOnEdge') 224 222 # get positions, lengths, surfaces and weights 225 223 le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle') … … 296 294 ucov[edge,:] = self.de[edge]*(ulon[edge,:]*math.cos(angle)+ulat[edge,:]*math.sin(angle)) 297 295 return ucov 298 def plot_i(self,data): 299 plot(self.primal,data) 300 def plot_v(self,data): 301 plot(self.dual,data) 302 def plot_e(self,data): 303 plot(self.triedge,data) 296 def plot_i(self,data, **kwargs): 297 plot(self.primal,data, **kwargs) 298 def plot_v(self,data, **kwargs): 299 plot(self.dual,data, **kwargs) 300 def plot_e(self,data, **kwargs): 301 plot(self.triedge,data, **kwargs) 302 303 #-------------------------------------- Mesh reordering ------------------------------------------# 304 305 # NB : Indices start at 0 306 # permutations map an old index to a new index : perm[old_index]=new_index 307 # enumerate(perm) is the list [ (0,perm[0]) , (1,perm[1]), ...] = [(old, new) ... ] 308 309 def reorder_values(perm, *datas): 310 # data in datas : ndarray containing indices ; perm : permutation to apply to values 311 for data in datas: 312 for x in np.nditer(data, op_flags=['readwrite']): 313 x[...]=perm[x] 314 315 def reorder_indices(perm, *datas): 316 # data in datas : ndarray containing values ; perm : permutation to apply to indices 317 for data in datas: 318 data_out, ndim = data.copy(), data.ndim 319 if ndim == 1: 320 for old,new in enumerate(perm): 321 data_out[new]=data[old] 322 if ndim == 2: 323 for old,new in enumerate(perm): 324 data_out[new,:]=data[old,:] 325 data[:]=data_out[:] 326 327 def key_to_perm(keys): 328 # keys : maps an old index to a key 329 # returns : a list mapping old indices to new indices obtained by sorting keys in ascending order 330 ordered = [old for key,old in sorted(zip(keys,range(len(keys))))] 331 # ordered maps new indices to old indices 332 # but we want to map old indices to new indices 333 perm = list(range(len(keys))) 334 for new,old in enumerate(ordered): 335 perm[old]=new 336 return perm # maps old indices to new indices 337 338 def toint(x): return np.int32((x+1.)*65536) 339 340 def morton_index(x,y,z): 341 nn, xx, yy, zz = x.size, toint(x), toint(y), toint(z) 342 idx=np.zeros(nn, dtype=np.int64) 343 ker.dynamico_morton_encode(nn, xx,yy,zz, idx) 344 return idx 304 345 305 346 #-------------------------------------- Mesh partitioning ------------------------------------------# 347 348 # NB : Indices start at 0 306 349 307 350 def partition_from_stencil(owner2, degree, stencil): … … 331 374 332 375 #---------------------------------- Stencil management ---------------------------------------# 376 377 # NB : Indices start at 0 333 378 334 379 # Class Stencil represents an adjacency relationship (e.g. cell=>edges) -
codes/icosagcm/devel/Python/src/partition.c
r620 r674 29 29 tpwgts, ubvec, options, &edgecut, part, &comm); 30 30 } 31 32 33 // http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ 34 // method to seperate bits from a given integer 3 positions apart 35 static inline uint64_t splitBy3(unsigned int a){ 36 uint64_t x = a & 0x1fffff; // we only look at the first 21 bits 37 x = (x | x << 32) & 0x1f00000000ffff; // shift left 32 bits, OR with self, and 00011111000000000000000000000000000000001111111111111111 38 x = (x | x << 16) & 0x1f0000ff0000ff; // shift left 32 bits, OR with self, and 00011111000000000000000011111111000000000000000011111111 39 x = (x | x << 8) & 0x100f00f00f00f00f; // shift left 32 bits, OR with self, and 0001000000001111000000001111000000001111000000001111000000000000 40 x = (x | x << 4) & 0x10c30c30c30c30c3; // shift left 32 bits, OR with self, and 0001000011000011000011000011000011000011000011000011000100000000 41 x = (x | x << 2) & 0x1249249249249249; 42 return x; 43 } 44 45 static inline uint64_t mortonEncode_magicbits(unsigned int x, unsigned int y, unsigned int z){ 46 uint64_t answer = 0; 47 answer |= splitBy3(x) | splitBy3(y) << 1 | splitBy3(z) << 2; 48 return answer; 49 } 50 51 void dynamico_morton_encode(int n, const int *x, const int *y, const int *z, uint64_t *m) 52 { 53 for(int i=0; i<n; i++) 54 m[i]=mortonEncode_magicbits(x[i],y[i],z[i]); 55 } -
codes/icosagcm/devel/Python/src/unstructured.pyx
r663 r674 44 44 except OSError: 45 45 print """ 46 Unable to load shared library 'lib kernels.so' !46 Unable to load shared library 'libicosa.so' ! 47 47 """ 48 48 raise … … 59 59 ['dynamico_caldyn_unstructured', c_double, c_void_p,20], 60 60 ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p], 61 ['dynamico_morton_encode', c_int,c_void_p,4] 61 62 ]) 62 63
Note: See TracChangeset
for help on using the changeset viewer.