Index: /codes/mesh_generation/build.sh
===================================================================
--- /codes/mesh_generation/build.sh	(revision 974)
+++ /codes/mesh_generation/build.sh	(revision 974)
@@ -0,0 +1,6 @@
+cd scvt
+python3 setup.py build_ext --inplace
+mv scvt/*.so .
+nm *.so | grep ' T '
+rmdir scvt
+rm -rf build
Index: /codes/mesh_generation/scvt/__init__.py
===================================================================
--- /codes/mesh_generation/scvt/__init__.py	(revision 974)
+++ /codes/mesh_generation/scvt/__init__.py	(revision 974)
@@ -0,0 +1,1 @@
+
Index: /codes/mesh_generation/scvt/delaunay.py
===================================================================
--- /codes/mesh_generation/scvt/delaunay.py	(revision 974)
+++ /codes/mesh_generation/scvt/delaunay.py	(revision 974)
@@ -0,0 +1,60 @@
+import numpy as np
+from scipy.spatial.qhull import _QhullUser
+from scipy.spatial.qhull import _Qhull
+
+from scvt.timer import Timer
+from scvt.sphere import unit_vector
+from scvt.topology import simplex
+
+########################################################
+
+# custom Delaunay class, slightly more efficient
+
+class Delaunay(_QhullUser):
+    def __init__(self, points):
+        qhull_options = b"Qbb Qc Qz Q12"
+        # Run qhull
+        qhull = _Qhull(b"d", points, qhull_options, required_options=b"Qt")
+        _QhullUser.__init__(self, qhull)
+    def _update(self, qhull):
+        qhull.triangulate()
+        self.paraboloid_scale, self.paraboloid_shift = \
+                               qhull.get_paraboloid_shift_scale()
+        self.simplices, self.neighbors, self.equations, self.coplanar = qhull.get_simplex_facet_array()
+        _QhullUser._update(self, qhull)
+
+########################################################        
+
+def testn(z): return z>0.
+def tests(z): return z<=0.
+
+def north_south(points, verbose):
+    with Timer('in stereographic projection', verbose):
+        north,south = [],[]
+        for i in range(points.shape[0]):
+            xyz=unit_vector(points[i,:])
+            points[i,:]=xyz
+            x,y,z=xyz
+            zn, zs = 1.+z, 1.-z
+            if(zn>0.): north.append([i,[x/zn,y/zn]])
+            if(zs>0.):  south.append([i,[x/zs,y/zs]])
+        indn, xyn = zip(*north)
+        inds, xys = zip(*south)
+    with Timer('in planar Delaunay', verbose):
+        trin, tris = Delaunay(np.asarray(xyn)), Delaunay(np.asarray(xys))
+
+    return points[:,2], (indn, trin, testn), (inds, tris, tests)
+
+def enum_simplices(z,domains):
+    # domains is a list of tuples (ind, tri, test)
+    for ind, tri, test in domains:
+        for i,j,k in tri.simplices:
+            ii, jj, kk = ind[i], ind[j], ind[k] # global indices
+            zi, zj, zk = z[ii], z[jj], z[kk]
+            if test(zi) | test(zj) | test(zk): yield simplex(ii,jj,kk)
+
+def delaunay(points, verbose=False):
+    z, north, south = north_south(points, verbose)
+    with Timer('enumerating simplices', verbose):
+        simplices = list(set( enum_simplices(z, (north, south) ) ))
+    return simplices
Index: /codes/mesh_generation/scvt/lloyd.py
===================================================================
--- /codes/mesh_generation/scvt/lloyd.py	(revision 974)
+++ /codes/mesh_generation/scvt/lloyd.py	(revision 974)
@@ -0,0 +1,61 @@
+import numpy as np
+
+from scvt import sphere, sphere_fast
+from scvt.timer import Timer
+from scvt.delaunay import delaunay
+from scvt.topology import all_edges_regions, regions_vertices
+from scvt.sphere import areas, circum
+        
+def flatten_regions(regions, edges):
+    degree = [len(region) for region in regions]
+    vert = list(regions_vertices(regions, edges))
+    return np.asarray(degree), np.asarray(vert)
+
+def flatten_simplices(simplices):
+    for i,j,k in simplices:
+        yield i
+        yield j
+        yield k        
+
+def lloyd(points, weight, verbose):
+    with Timer('in one Lloyd iteration (total)', verbose):
+        N = points.shape[0]
+        simplices = delaunay(points,verbose)
+#    dict_simp=dict(enumerate(simplices))
+        with Timer('in np.asarray(simplices)', verbose):
+            simp = np.asarray(simplices)
+        with Timer('enumerating edges', verbose):
+            simplex_index={ simplex:n for n,simplex in enumerate(simplices) }
+            edges, regions = all_edges_regions(simplices)
+
+#    with Timer('computing circumcenters (Python)', verbose):
+#        vertices = [circum(points[a,:],points[b,:],points[c,:]) for a,b,c in simplices]
+#        vertices = np.asarray(vertices)
+#    with Timer('computing circumcenters (jit)', verbose):
+#        vertices = sphere.circum_fast(simp, points)
+        with Timer('computing circumcenters (C)', verbose):
+            vertices=np.zeros((len(simplices),3))
+            sphere_fast.circum_fast(simp, points,vertices)
+        with Timer('computing barycenters (C)', verbose):
+            regions_size, vert = flatten_regions(regions, edges)
+#        degree, counts = np.unique(regions_size, return_counts=True)
+#        print("Edge counts of Voronoi cells :", dict(zip(degree, counts)))
+            w_i,w_v = weight(points), weight(vertices)
+            bary = points.copy()
+            sphere_fast.barycenters_fast(regions_size, vert, points, w_i, vertices, w_v, bary)
+#        bary = sphere.barycenters_fast(regions_size, vert, points, vertices, weight)
+#        bary = sphere.barycenters(regions, points, edges, vertices)
+    
+    if verbose: 
+        area = np.asarray(list(areas(points, simplices)))
+        print('%d cells, %d edges, %d triangles'%(N, len(edges), len(simplices)))
+        print('Total area = 4*pi + ',area.sum()-4*np.pi)
+
+    if N<50: 
+        print(edges)
+        a,b,c=points[simplices[0],:]
+        test(a,b,c)
+        test(a,c,b)
+        test(c,a,b)
+
+    return bary, vertices, edges
Index: /codes/mesh_generation/scvt/setup.py
===================================================================
--- /codes/mesh_generation/scvt/setup.py	(revision 974)
+++ /codes/mesh_generation/scvt/setup.py	(revision 974)
@@ -0,0 +1,5 @@
+from distutils.core import setup
+from Cython.Build import cythonize
+
+setup(name='Spherical geometry primitives, compiled for speed',
+      ext_modules=cythonize("sphere_fast.pyx", annotate=True))
Index: /codes/mesh_generation/scvt/sphere.py
===================================================================
--- /codes/mesh_generation/scvt/sphere.py	(revision 974)
+++ /codes/mesh_generation/scvt/sphere.py	(revision 974)
@@ -0,0 +1,116 @@
+import numpy as np
+from numpy import arcsin, arctan, tan, sqrt, cross
+from matplotlib.collections import LineCollection
+import matplotlib.pyplot as plt
+
+import numba
+jit=numba.jit(nopython=True, parallel=True)
+
+############### spherical primitives #################
+
+def dotprod(a,b): return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]
+@jit
+def norm2(a): return a[0]*a[0]+a[1]*a[1]+a[2]*a[2]
+@jit
+def norm(a): return np.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2])
+def unit_vector(a): return a/np.sqrt((a*a).sum())
+def arc(a,b): return 2.*np.arcsin(.5*norm(a-b))
+
+@jit
+def cross_jit(a,b, c):
+    c[0]=a[1]*b[2]-b[1]*a[2]
+    c[1]=a[2]*b[0]-b[2]*a[0]
+    c[2]=a[0]*b[1]-b[0]*a[1]
+
+def area(a,b,c):
+# http://mathworld.wolfram.com/LHuiliersTheorem.html
+    a,b,c = .25*arc(b,c), .25*arc(c,a), .25*arc(a,b)
+    area = sqrt(tan(a+b+c)*tan(a+b-c)*tan(b+c-a)*tan(a+c-b))
+#    print("area(a,b,c):",a,b,c, a+b-c, b+c-a, a+c-b, area)
+    return 4.*arctan(area)
+
+def areas(points, simplices):
+    for simplex in simplices:
+        a,b,c = points[simplex,:]
+        yield area(a,b,c)
+    
+def circum(A,B,C):
+# see Dubos et al., GMD 2015, appendix
+# the computation is stable, and valid whether a,b,c is clockwise or not
+    A,B,C=map(unit_vector, (A,B,C))
+    b,c = B-A, C-A
+    p1,p2 = cross(b,c), norm2(c)*b-norm2(b)*c
+    p = A+.5*cross(p1,p2)/norm2(p1)
+    p = p/norm(p)
+    a,b,c=norm(p-A),norm(p-B),norm(p-C)
+    if abs(a-b)>1e-15:
+        print("Check circum", a-b, b-c, a-c)
+        raise(ValueError)
+    return p
+
+@jit
+def circum_fast(simp, points):
+    """ simp : int[M,3] describing M simplices"""
+    M=simp.shape[0]
+    circ = np.zeros((M,3))
+    p1, p2, p = np.zeros(3), np.zeros(3), np.zeros(3)
+    for n in range(M):
+        i,j,k = simp[n,:]
+        A,B,C = points[i,:], points[j,:], points[k,:]
+        b,c = B-A, C-A
+        cross_jit(b,c, p1)
+        p2 = norm2(c)*b-norm2(b)*c
+        cross_jit(p1,p2, p)
+        p = A+p*(.5/norm2(p1))
+        circ[n,:]=p*(1./norm(p))
+    return circ
+
+def barycenters_fast(regions_size, vert, points, vertices, weight):
+    # print("== 4 x number of edges =%d"%len(vert))
+    N = points.shape[0]
+    bary = np.zeros((N,3))
+    i=0 # running index in vert1, vert2
+    A=0.
+    for n in range(N):
+        point = points[n,:]
+        b = 0.
+        for j in range(regions_size[n]):
+            v1 = vertices[vert[i],:]
+            i = i+1
+            v2 = vertices[vert[i],:]
+            i = i+1
+            dA=area(point, v1, v2)
+            A = A+dA
+            bb = unit_vector(point+v1+v2)
+            b = b + dA*bb*weight(bb)
+        bary[n,:] = unit_vector(b)
+    print("== Total area = 4*pi+%f"%(A-4*np.pi)) 
+    return bary
+        
+def barycenters(regions, points, edges, vertices):
+    N = points.shape[0]
+    bary = np.zeros((N,3))
+    A=0.
+    for n in range(N):
+        point, region = points[n,:], regions[n] # generator and Voronoi region
+        b = 0.
+        for edge in region:
+            (n1,p1), (n2,p2) = edges[edge]
+            v1, v2 = vertices[n1,:], vertices[n2,:]
+            dA=area(point, v1, v2)
+            if N<10 :
+                print("n, n1, n2", n,n1,n2)
+                print("point,v1,v2,dA=",point,v1,v2,dA)
+            A = A+dA
+            bb = unit_vector(point+v1+v2)
+            b = b + dA*bb
+        bary[n,:] = unit_vector(b)
+    if N<100: print("==== Total area = 4*pi+%f"%(A-4*np.pi))
+    return bary
+
+def plot_mesh(points,edges,orient, color):
+    xy=points[:,0:2]
+    segments=[ xy[[i, j],:] for i,j in edges if dotprod(orient,points[i,:])>0. ]
+    ax = plt.axes(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1))
+    ax.add_collection(LineCollection(segments, color=color))
+
Index: /codes/mesh_generation/scvt/sphere_fast.h
===================================================================
--- /codes/mesh_generation/scvt/sphere_fast.h	(revision 974)
+++ /codes/mesh_generation/scvt/sphere_fast.h	(revision 974)
@@ -0,0 +1,241 @@
+#include <math.h>
+
+enum { batchsize=8 };
+static double Ax[batchsize], Ay[batchsize], Az[batchsize];
+static double Bx[batchsize], By[batchsize], Bz[batchsize]; 
+static double Cx[batchsize], Cy[batchsize], Cz[batchsize]; 
+static double Dx[batchsize], Dy[batchsize], Dz[batchsize]; 
+
+/*
+
+def norm2(a): return a[0]*a[0]+a[1]*a[1]+a[2]*a[2]
+def norm(a): return np.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2])
+
+def cross_jit(a,b, c):
+    c[0]=a[1]*b[2]-b[1]*a[2]
+    c[1]=a[2]*b[0]-b[2]*a[0]
+    c[2]=a[0]*b[1]-b[0]*a[1]
+
+def circum(A,B,C):
+# see Dubos et al., GMD 2015, appendix
+# the computation is stable, and valid whether a,b,c is clockwise or not
+    b,c = B-A, C-A
+    p1,p2 = cross(b,c), norm2(c)*b-norm2(b)*c
+    p = A+.5*cross(p1,p2)/norm2(p1)
+    return p/norm(p)
+*/
+
+static inline void circum_main(int N)
+{
+  // Any decent compiler should vectorize this
+  for(int i=0 ; i<N ; i++)
+    {
+      //    b,c = B-A, C-A
+      double bx,by,bz,cx,cy,cz;
+      bx=Bx[i]-Ax[i];
+      by=By[i]-Ay[i];
+      bz=Bz[i]-Az[i];
+      cx=Cx[i]-Ax[i];
+      cy=Cy[i]-Ay[i];
+      cz=Cz[i]-Az[i];
+      // p = cross(b,c)
+      double px,py,pz;
+      px = by*cz-bz*cy;
+      py = bz*cx-bx*cz;
+      pz = bx*cy-by*cx;
+      // q = norm2(c)*b-norm2(b)*c
+      double b2, c2, qx,qy,qz ;
+      b2 = bx*bx+by*by+bz*bz;
+      c2 = cx*cx+cy*cy+cz*cz;
+      qx = c2*bx-b2*cx;
+      qy = c2*by-b2*cy;
+      qz = c2*bz-b2*cz;
+      // r = A+.5*cross(p,q)/norm2(p)
+      double p2,rx,ry,rz;
+      p2 = .5/(px*px+py*py+pz*pz);
+      rx = Ax[i]+p2*(py*qz-pz*qy);
+      ry = Ay[i]+p2*(pz*qx-px*qz);
+      rz = Az[i]+p2*(px*qy-py*qx);
+      // circ[i]=r/norm(r)
+      double r=1./sqrt(rx*rx+ry*ry+rz*rz);
+      Dx[i] = r*rx;
+      Dy[i] = r*ry;
+      Dz[i] = r*rz;
+    }
+}
+
+static inline void circum_copyin(int N, long *ind, double *points)
+{
+  for(int n=0 ; n<N ; n++)
+    {
+      int i; const double *ptr;
+      i=*ind++; ptr=points+3*i;
+      Ax[n]=*ptr++;
+      Ay[n]=*ptr++;
+      Az[n]=*ptr++;
+      i=*ind++; ptr=points+3*i;
+      Bx[n]=*ptr++;
+      By[n]=*ptr++;
+      Bz[n]=*ptr++;
+      i=*ind++; ptr=points+3*i;
+      Cx[n]=*ptr++;
+      Cy[n]=*ptr++;
+      Cz[n]=*ptr++;
+    }
+}
+
+static inline void circum_copyout(int N, double *circ)
+{
+  for(int i=0 ; i<N ; i++)
+    {
+      *circ++ = Dx[i];
+      *circ++ = Dy[i];
+      *circ++ = Dz[i];
+    }
+}
+
+static inline void circum_batch(double *points, double *circ)
+{
+  for(int i=0 ; i<batchsize ; i++)
+    {
+      Ax[i]=*points++;
+      Ay[i]=*points++;
+      Az[i]=*points++;
+      Bx[i]=*points++;
+      By[i]=*points++;
+      Bz[i]=*points++;
+      Cx[i]=*points++;
+      Cy[i]=*points++;
+      Cz[i]=*points++;
+    }
+  circum_main(batchsize);
+  circum_copyout(batchsize, circ);
+}
+
+static inline void circum_fast_(int N, long *ind, double *points, double *circ)
+{
+  // work by batches of size batchsize
+  while(N>batchsize)
+    {
+      circum_copyin(batchsize, ind, points);
+      circum_main(batchsize);
+      circum_copyout(batchsize, circ);
+      ind  += 3*batchsize;
+      circ += 3*batchsize;
+      N-=batchsize;
+    }
+  // finish the job
+  circum_copyin(N, ind, points);
+  circum_main(N);
+  circum_copyout(N, circ);
+}
+
+static inline void fast_scale_vec(double scale, double * vec)
+{
+  vec[0]*=scale;
+  vec[1]*=scale;
+  vec[2]*=scale;
+}
+ 
+/*
+
+def arc(a,b): return 2.*np.arcsin(.5*norm(a-b))
+
+def area(a,b,c):
+# http://mathworld.wolfram.com/LHuiliersTheorem.html
+    a,b,c = .25*arc(b,c), .25*arc(c,a), .25*arc(a,b)
+    return 4.*arctan(sqrt(tan(a+b+c)*tan(a+b-c)*tan(b+c-a)*tan(a+c-b)))
+
+def barycenters_fast(regions_size, vert, points, vertices, weight):
+    N = points.shape[0]
+    bary = np.zeros((N,3))
+    i=0 # running index in vert1, vert2
+    for n in range(N):
+        point = points[n,:]
+        b = 0.
+        for j in range(regions_size[n]):
+            v1 = vertices[vert[i],:]
+            i = i+1
+            v2 = vertices[vert[i],:]
+            i = i+1
+            dA=area(point, v1, v2)
+            bb = unit_vector(point+v1+v2)
+            b = b + dA*bb*weight(bb)
+        bary[n,:] = unit_vector(b)
+    return bary
+        
+*/
+
+static inline double arc4(double dx, double dy, double dz) // 1/4 geodesic length, as used in L'Huillier's formula
+{
+  double dist;
+  //  printf("dx dy dz = %f %f %f\n", dx,dy,dz);
+  dist=dx*dx+dy*dy+dz*dz;
+  //  printf("dist = %f\n", dist);
+  dist=sqrt(dist);
+  //  printf("dist = %f\n", dist);
+  dist = .5*asin(.5*dist);
+  //  printf("dist = %f\n", dist);
+  return dist;
+}
+
+
+static inline double barycenters_fast_(int N, const long *degree, const long *ind, 
+				       const double *point, const double *w_i,
+				       const double *vertex, const double *w_v,
+				       double *bary)
+{
+  double A=0; // total area
+  // loop over voronoi cells
+  for(int cell=0; cell<N ; cell++)
+    {
+      double aw,ax,ay,az; // weight and xyz coordinates of generator of voronoi cell
+      int deg;
+      double gx=0., gy=0., gz=0.; // barycenter of current voronoi cell
+      aw=*w_i++;
+      ax=*point++;
+      ay=*point++;
+      az=*point++;
+      deg=*degree++;
+      //      if(cell<10) printf("deg[%d]=%d\n", cell, deg);
+
+      for(int edge=0; edge<deg ; edge++) // loop over edges of voronoi cell
+	{
+	  int i,j; const double *ptr;
+	  double bw,bx,by,bz, cw,cx,cy,cz;
+	  i=*ind++; ptr=vertex+3*i; // first vertex
+	  bw=w_v[i];
+	  bx=*ptr++;
+	  by=*ptr++;
+	  bz=*ptr++;
+	  j=*ind++; ptr=vertex+3*j; // second vertex
+	  cw=w_v[j];
+	  cx=*ptr++;
+	  cy=*ptr++;
+	  cz=*ptr++;
+	  //	  printf("cell, vertex1, vertex2 : %d %d %d\n", cell, i, j);
+	  //	  printf("cell    x y z : %f %f %f\n", ax,ay,az);
+	  //	  printf("vertex1 x y z : %f %f %f\n", bx,by,bz);
+	  //	  printf("vertex2 x y z : %f %f %f\n", cx,cy,cz);
+
+	  double a,b,c,dA; // 1/4 geodesic distances, as needed in L'Huillier's formula
+	  // http://mathworld.wolfram.com/LHuiliersTheorem.html
+	  a=arc4(bx-cx,by-cy,bz-cz);
+	  b=arc4(ax-cx,ay-cy,az-cz);
+	  c=arc4(ax-bx,ay-by,az-bz);
+	  dA=4.*atan(sqrt(tan(a+b+c)*tan(a+b-c)*tan(b+c-a)*tan(a+c-b)));
+	  A+=dA;
+	  gx+=dA*(aw*ax+bw*bx+cw*cx);
+	  gy+=dA*(aw*ay+bw*by+cw*cy);
+	  gz+=dA*(aw*az+bw*bz+cw*cz);
+	  //	  printf("Triangle sides : %f %f %f\n", a,b,c);
+	  //	  printf("Local area : dA=%f\n", dA);
+	}
+      double g=1./sqrt(gx*gx+gy*gy+gz*gz);
+      *bary++ = gx*g;
+      *bary++ = gy*g;
+      *bary++ = gz*g;
+    }
+  // printf("Total area : A=%f\n", A);
+  return A; // total area should be 4*pi
+}
Index: /codes/mesh_generation/scvt/sphere_fast.pyx
===================================================================
--- /codes/mesh_generation/scvt/sphere_fast.pyx	(revision 974)
+++ /codes/mesh_generation/scvt/sphere_fast.pyx	(revision 974)
@@ -0,0 +1,39 @@
+cdef extern from "sphere_fast.h":
+    cdef void fast_scale_vec(double, double[])
+    cdef void circum_batch(double*, double*)
+    cdef void circum_fast_(int, long*, double*, double*)
+    cdef void barycenters_fast_(int, long*, long*, double*, double*, double*, double*, double*)
+
+def scale_vec(double scale, double[::1] vec):
+    fast_scale_vec(scale, &vec[0])
+
+def circum_test(double[:,::1] pts, double[:,::1] circ):
+    circum_batch(&pts[0,0], &circ[0,0])
+
+def circum_fast(long[:,::1] ind, double[:,::1] pts, double[:,::1] circ):
+    cdef int N = ind.shape[0] # number of circumcenters to compute
+#    print('Computing %d circumcenters'%N)
+    circum_fast_(N, &ind[0,0], &pts[0,0], &circ[0,0])
+
+def barycenters_fast(long[::1] degree, long[::1] index, 
+                     double[:,::1] point, double[::1] w_i,
+		     double[:,::1] vertex, double[::1] w_v,
+		     double[:,::1] bary):
+    cdef int N=degree.shape[0] # number of voronoi cells
+
+    cdef int sum=0
+    for i in range(N):
+      sum += degree[i]
+
+#    print(degree[0:10])
+#    print('Computing %d barycenters'%N)
+#    print('degree.sum() = %d'%sum )
+#    print('index.size = %d'%index.size )
+#    print('point.shape = ', point.shape )
+#    print('vertex.shape = ', vertex.shape )
+#    print('bary.shape = ', bary.shape )
+    
+    return barycenters_fast_(N, &degree[0], &index[0], 
+                             &point[0,0], &w_i[0],
+			     &vertex[0,0], &w_v[0],
+			     &bary[0,0])
Index: /codes/mesh_generation/scvt/timer.py
===================================================================
--- /codes/mesh_generation/scvt/timer.py	(revision 974)
+++ /codes/mesh_generation/scvt/timer.py	(revision 974)
@@ -0,0 +1,12 @@
+import time
+
+############### profiling ##################
+
+class Timer(object):
+    def __init__(self, msg, verbose=True):
+        self.msg, self.verbose = msg, verbose
+    def __enter__(self):
+        self.tstart = time.time()
+    def __exit__(self, type, value, traceback):
+        if self.verbose:
+            print('Time spent %s : %s' % (self.msg, time.time() - self.tstart))
Index: /codes/mesh_generation/scvt/topology.py
===================================================================
--- /codes/mesh_generation/scvt/topology.py	(revision 974)
+++ /codes/mesh_generation/scvt/topology.py	(revision 974)
@@ -0,0 +1,60 @@
+import numpy as np
+from scvt.timer import Timer
+
+########################################################        
+
+def simplex(ii,jj,kk):
+    if ii>jj: ii,jj = jj,ii
+    if jj>kk: jj,kk = kk,jj
+    if ii>jj: ii,jj = jj,ii
+    # now ii<jj<kk
+    return (ii,jj,kk)
+
+def enum_edges(simplices):
+    for ii,jj,kk in simplices:
+        # ii<jj<kk
+        yield (ii,jj)
+        yield (jj,kk)
+        yield (ii,kk)
+
+def insert_edge(edges, regions, edge,kk,n):
+    if edge in edges:
+        edges[edge]=(edges[edge], (n,kk))
+    else:
+        edges[edge]=(n,kk)
+    ii,jj = edge
+    if ii in regions:
+        regions[ii].add(edge)
+    else:
+        regions[ii]={edge}
+    if jj in regions:
+        regions[jj].add(edge)
+    else:
+        regions[jj]={edge}
+
+def all_edges_regions(simplices):
+    # returns edges, regions
+    # edge = ordered pair of generators
+    # edges is a dictionary edge -> ((n1,point1, (n2,point2))
+    # where n is the index of a simplex containing edge
+    # and point is the index of the third point of that simplex
+    # regions is a dictionary : point -> dict of edges 
+    edges,regions = {},{}
+    for n,(ii,jj,kk) in enumerate(simplices):
+        # note that ii<jj<kk
+        # edges (a,b) must have a<b
+        insert_edge(edges,regions, (ii,jj),kk,n)
+        insert_edge(edges,regions, (jj,kk),ii,n)
+        insert_edge(edges,regions, (ii,kk),jj,n)
+    return edges, [ regions[i] for i in range(len(regions)) ]
+
+def dual_edges(edges):
+    for (ii,jj),((n1,kk1),(n2,kk2)) in edges.items():
+        yield (n1,n2)        
+
+def regions_vertices(regions, edges):
+    for region in regions:
+        for edge in region:
+            (n1,p1), (n2,p2) = edges[edge]
+            yield n1
+            yield n2
Index: /codes/mesh_generation/test_scvt2.py
===================================================================
--- /codes/mesh_generation/test_scvt2.py	(revision 974)
+++ /codes/mesh_generation/test_scvt2.py	(revision 974)
@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from numba import jit
+
+from scvt import topology
+from scvt import sphere
+from scvt import delaunay
+from scvt import lloyd
+
+def test(a,b,c):
+    p=circum(a,b,c)
+    print(norm(a-p),norm(b-p),norm(c-p)) # should be equal
+
+@jit
+def weight(xyz):
+    y = xyz[:,1];
+    return 1.+15.*np.exp(-16.*y*y)
+
+def main(nbp,inner,outer,filename=None):
+    N=10*nbp*nbp+2
+    points=np.random.randn(N,3)
+    for i in range(outer):
+        verbose = True
+        print('Outer Lloyd iteration %d/%d'%(i,outer))
+        for j in range(inner):
+            points, vertices, edges = lloyd.lloyd(points, weight, verbose)
+            verbose = False
+    plt.figure(figsize=(10,10))
+    sphere.plot_mesh(vertices, topology.dual_edges(edges), orient, "k")
+    if N<1000:
+        sphere.plot_mesh(points,edges,orient,"r")
+    if filename is None:
+        plt.show()
+    else:
+        plt.savefig(filename)
+
+
+
+orient=np.asarray([0,0,1.])
+
+################## Manual tests with octahedron ####################
+
+points=np.asarray([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.],
+                   [-1.,0.,0.],[0.,-1.,0.],[0.,0.,-1.]])
+simplices = delaunay.delaunay(points)
+print("== Simplices ==", simplices)
+edges, regions = topology.all_edges_regions(simplices)
+print('%d cells, %d edges, %d triangles'%(points.shape[0], len(edges), len(simplices)))
+print(edges)
+for n,region in enumerate(regions):
+    print("Voronoi cell %d :"%n, region)
+for n,region in enumerate(regions):
+    print("Voronoi cell %d :"%n, [edges[edge] for edge in region])
+ 
+vertices = [sphere.circum(points[a,:],points[b,:],points[c,:]) for a,b,c in simplices]
+vertices = np.asarray(vertices)
+print("== Vertices ==", vertices)
+bary = sphere.barycenters(regions, points, edges, vertices)
+
+regions_size, vert = lloyd.flatten_regions(regions, edges)
+degree, counts = np.unique(regions_size, return_counts=True)
+print("Edge counts of Voronoi cells :", dict(zip(degree, counts)))
+
+####################### Larger grid #######################
+#main(8, 1,1)
+#main(8, 100,5)
+main(16, 100,5, 'nbp16.png')
+main(40, 100,10, 'nbp40.png')
