Changeset 674 for codes/icosagcm

Timestamp:

01/31/18 15:14:35 (6 years ago)

Author:

dubos

Message:

devel/unstructured : mesh reordering (Morton)

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/meshes.py (modified) (6 diffs)
• src/partition.c (modified) (1 diff)
• src/unstructured.pyx (modified) (2 diffs)
• test/fig_morton (added)
• test/py/reorder.py (added)

codes/icosagcm/devel/Python/dynamico/meshes.py

@@ -5 +5 @@
 import matplotlib.tri as tri
 import matplotlib.pyplot as plt
-from unstructured import init_mesh, list_stencil
+from unstructured import init_mesh, list_stencil, ker
 
 radian=180/math.pi
@@ -174 +174 @@
     return ne
 
-def plot(tri,data):
+def plot(tri,data, **kwargs):
     plt.figure(figsize=(12,4))
     plt.gca().set_aspect('equal')
-    plt.tricontourf(tri, data, 20)
+    plt.tricontourf(tri, data, 20, **kwargs)
     plt.colorbar()
     plt.ylim((-90,90))
@@ -210 +210 @@
 
         # get indices for stencils 
-        # primal -> vertices (unused)
-        primal_vertex, dual_vertex = getvars('verticesOnCell','cellsOnVertex')
         # primal <-> edges
         primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge')
@@ -221 +219 @@
         dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
         # primal <-> dual, edges <-> edges (TRiSK)
-        dual_vertex, trisk = getvars('cellsOnVertex','edgesOnEdge')
+        primal_vertex, dual_vertex, trisk = getvars('verticesOnCell','cellsOnVertex','edgesOnEdge')
         # get positions, lengths, surfaces and weights
         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))
         return ucov
-    def plot_i(self,data):
-        plot(self.primal,data)
-    def plot_v(self,data):
-        plot(self.dual,data)
-    def plot_e(self,data):
-        plot(self.triedge,data)
+    def plot_i(self,data, **kwargs):
+        plot(self.primal,data, **kwargs)
+    def plot_v(self,data, **kwargs):
+        plot(self.dual,data, **kwargs)
+    def plot_e(self,data, **kwargs):
+        plot(self.triedge,data, **kwargs)
+
+#--------------------------------------  Mesh reordering  ------------------------------------------#
+
+# NB : Indices start at 0
+# permutations map an old index to a new index : perm[old_index]=new_index
+# enumerate(perm) is the list [ (0,perm[0]) , (1,perm[1]), ...] = [(old, new) ... ]
+
+def reorder_values(perm, *datas): 
+    # data in datas : ndarray containing indices ; perm : permutation to apply to values
+    for data in datas:
+        for x in np.nditer(data, op_flags=['readwrite']):
+            x[...]=perm[x]
+
+def reorder_indices(perm, *datas): 
+    # data in datas : ndarray containing values ; perm : permutation to apply to indices
+    for data in datas:
+        data_out, ndim = data.copy(), data.ndim
+        if ndim == 1:
+            for old,new in enumerate(perm):
+                data_out[new]=data[old]
+        if ndim == 2:
+            for old,new in enumerate(perm):
+                data_out[new,:]=data[old,:]
+        data[:]=data_out[:]
+
+def key_to_perm(keys): 
+    # keys : maps an old index to a key
+    # returns : a list mapping old indices to new indices obtained by sorting keys in ascending order
+    ordered = [old for key,old in sorted(zip(keys,range(len(keys))))] 
+    # ordered maps new indices to old indices
+    # but we want to map old indices to new indices
+    perm = list(range(len(keys)))
+    for new,old in enumerate(ordered):
+        perm[old]=new
+    return perm  # maps old indices to new indices
+
+def toint(x): return np.int32((x+1.)*65536)
+    
+def morton_index(x,y,z):
+    nn, xx, yy, zz = x.size, toint(x), toint(y), toint(z)
+    idx=np.zeros(nn, dtype=np.int64)
+    ker.dynamico_morton_encode(nn, xx,yy,zz, idx)
+    return idx
 
 #-------------------------------------- Mesh partitioning ------------------------------------------#
+
+# NB : Indices start at 0
 
 def partition_from_stencil(owner2, degree, stencil):
@@ -331 +374 @@
 
 #---------------------------------- Stencil management ---------------------------------------#
+
+# NB : Indices start at 0
 
 # Class Stencil represents an adjacency relationship (e.g. cell=>edges)

codes/icosagcm/devel/Python/src/partition.c

@@ -29 +29 @@
                        tpwgts, ubvec, options, &edgecut, part, &comm); 
 }
+
+
+// http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ 
+// method to seperate bits from a given integer 3 positions apart
+static inline uint64_t splitBy3(unsigned int a){
+  uint64_t x = a & 0x1fffff; // we only look at the first 21 bits
+  x = (x | x << 32) & 0x1f00000000ffff;  // shift left 32 bits, OR with self, and 00011111000000000000000000000000000000001111111111111111
+  x = (x | x << 16) & 0x1f0000ff0000ff;  // shift left 32 bits, OR with self, and 00011111000000000000000011111111000000000000000011111111
+  x = (x | x << 8) & 0x100f00f00f00f00f; // shift left 32 bits, OR with self, and 0001000000001111000000001111000000001111000000001111000000000000
+  x = (x | x << 4) & 0x10c30c30c30c30c3; // shift left 32 bits, OR with self, and 0001000011000011000011000011000011000011000011000011000100000000
+  x = (x | x << 2) & 0x1249249249249249;
+  return x;
+}
+ 
+static inline uint64_t mortonEncode_magicbits(unsigned int x, unsigned int y, unsigned int z){
+  uint64_t answer = 0;
+  answer |= splitBy3(x) | splitBy3(y) << 1 | splitBy3(z) << 2;
+  return answer;
+}
+
+void dynamico_morton_encode(int n, const int *x, const int *y, const int *z, uint64_t *m)
+{
+  for(int i=0; i<n; i++)
+    m[i]=mortonEncode_magicbits(x[i],y[i],z[i]);
+}

codes/icosagcm/devel/Python/src/unstructured.pyx

@@ -44 +44 @@
 except OSError:
     print """
-Unable to load shared library 'libkernels.so' !
+Unable to load shared library 'libicosa.so' !
     """
     raise
@@ -59 +59 @@
                      ['dynamico_caldyn_unstructured', c_double, c_void_p,20],
                      ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p],
+                     ['dynamico_morton_encode', c_int,c_void_p,4]
                      ])