Index: /trunk/anuga_work/development/2010-projects/kv/compile_all.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/compile_all.py	(revision 8052)
+++ /trunk/anuga_work/development/2010-projects/kv/compile_all.py	(revision 8052)
@@ -0,0 +1,16 @@
+import os
+
+buildroot = os.getcwd()
+
+#Attempt to compile all extensions
+
+#os.chdir('utilities')
+execfile('compile.py')
+
+
+
+os.chdir(buildroot)    
+#execfile('test_all.py')
+    
+if sys.platform == 'win32':
+    raw_input('Press the RETURN key')
Index: /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity.py	(revision 8051)
+++ /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity.py	(revision 8052)
@@ -1,3 +1,3 @@
-from anuga.abstract_2d_finite_volumes.domain import Domain
+from anuga import Domain
 #from anuga.utilities.sparse import Sparse, Sparse_CSR
 from sparse import Sparse, Sparse_CSR #the new import
Index: /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity_ext.c
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity_ext.c	(revision 8051)
+++ /trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity_ext.c	(revision 8052)
@@ -40,5 +40,12 @@
 }
 
-int build_geo_structure(int n, int tot_len, double *centroids, int *neighbours, double *edgelengths, double *edge_midpoints, int *geo_indices, double *geo_values) {
+int build_geo_structure(int n, 
+                        int tot_len, 
+                        double *centroids, 
+                        long *neighbours, 
+                        double *edgelengths, 
+                        double *edge_midpoints, 
+                        long *geo_indices, 
+                        double *geo_values) {
     int i, edge, edge_counted, j;
 	double dist, this_x, this_y, other_x, other_y, edge_length;
@@ -75,5 +82,12 @@
 }
 
-int build_operator_matrix(int n, int tot_len, int *geo_indices, double *geo_values, double *h, double *boundary_heights, double *data, int *colind) {
+int build_operator_matrix(int n, 
+                          int tot_len, 
+                          long *geo_indices, 
+                          double *geo_values, 
+                          double *h, 
+                          double *boundary_heights, 
+                          double *data, 
+                          long *colind) {
 	int i, k, edge, j[4], *sorted_j, this_index;
 	double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
@@ -151,5 +165,11 @@
     Py_DECREF(mesh);
     
-    err = build_geo_structure(n,tot_len, (double *)centroid_coordinates -> data, (int *)neighbours -> data, (double *)edgelengths->data, (double *)edge_midpoint_coordinates -> data, (int *)geo_indices -> data, (double *)geo_values -> data);
+    err = build_geo_structure(n,tot_len, 
+                             (double *)centroid_coordinates -> data, 
+                             (long *)neighbours -> data, 
+                             (double *)edgelengths->data, 
+                             (double *)edge_midpoint_coordinates -> data, 
+                             (long *)geo_indices -> data, 
+                             (double *)geo_values -> data);
     if (err != 0) {
         PyErr_SetString(PyExc_RuntimeError, "Could not build geo structure");
@@ -190,5 +210,11 @@
 	colind = get_consecutive_array(kv_operator,"operator_colind");
 	
-	err = build_operator_matrix(n,tot_len, (int *)geo_indices -> data, (double *)geo_values -> data, (double *)h -> data, (double *)boundary_heights -> data, (double *)_data -> data, (int *)colind -> data);
+	err = build_operator_matrix(n,tot_len, 
+								(long *)geo_indices -> data, 
+								(double *)geo_values -> data, 
+								(double *)h -> data, 
+								(double *)boundary_heights -> data, 
+								(double *)_data -> data, 
+								(long *)colind -> data);
     if (err != 0) {
         PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
Index: unk/anuga_work/development/2010-projects/kv/kv_solver_new.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/kv_solver_new.py	(revision 8051)
+++ 	(revision )
@@ -1,312 +1,0 @@
-from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
-import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory
-from anuga.utilities.sparse import Sparse, Sparse_CSR
-from anuga.utilities.cg_solve import conjugate_gradient
-import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
-from generic_boundary_conditions import Dirichlet_boundary
-import numpy as num
-import anuga.utilities.log as log
-"""
-Kinematic Viscosity Class for ANUGA
-Lindon Roberts, 2010
-
-Possible extensions:
-(1) Currently assume that the line between centroids
-is parallel to the outward normal. Change formula
-for this.
-(2) Higher order approximation (eg. Least squares
-fitting or quadratic etc.)
-(3) Vectorise the for loops - multiply by diagonal matrices?
-(4) Tidy up direct application calling structure: __mul__ for matrices,
-then ensure correct setting of include_boundaries
-(5) More comprehensive testing (larger problems)
-"""
-
-
-"""
-How to use:
-
-(1) Direct application
-domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...)
-KV_Operator = Kinematic_Viscosity_Operator(domain)
-h = [n*1 numpy column vector]
-KV_Operator.apply_stage_heights(h)
-uh = [n*1 numpy column vector]
-vh = [n*1 numpy column vector]
-KV_Operator.set_qty_considered('u') #or KV_Operator.set_qty_considered(1)
-KV_Operator.build_boundary_vector()
-div_u = KV_Operator.apply_kv_operator(uh)
-KV_Operator.set_qty_considered('v') #or KV_Operator.set_qty_considered(2)
-KV_Operator.build_boundary_vector()
-div_v = KV_Operator.apply_kv_operator(vh)
-
-[Now div_u and div_v are (n,) numpy vectors]
-
-------------------------------------------------------------------------
-(2) Implicit use
-domain = anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain(...)
-KV_Operator = Kinematic_Viscosity_Operator(domain)
-h = [n*1 numpy column vector]
-KV_Operator.apply_stage_heights(h)
-div = [n*2 numpy matrix] #column 0 = div_u, column 1 = div_v
-velocities = KV_Operator.cg_solve(div)
-
-[Now velocities is a (n,2) numpy matrix = [(uh) (vh)]]
-
-If you just want to solve for (uh), say:
-KV_Operator.apply_stage_heights(h)
-div_u = [(n,) numpy vector]
-uh = KV_Operator.cg_solve(div_u, 'u') #or KV_Operator.cg_solve(div_u,1)
-"""
-
-class Kinematic_Viscosity_Operator:
-    ## @brief Initialise a KV operator
-    # @param domain The anuga.abstract_2d_finite_volumes.generic_domain.Generic_Domain for the problem
-    # @param triangle_areas Whether we should scale the divergence by 1/|Triangle_i|
-    # @param verbose
-    def __init__(self, domain, triangle_areas=True, verbose=False):
-        #Expose the domain attributes
-        self.domain = domain
-        self.mesh = domain.mesh
-        self.boundary = domain.boundary
-        self.boundary_enum = self.enumerate_boundary()
-        self.n = len(self.domain)
-        self.boundary_len = len(domain.boundary)
-        self.qty_considered = 1 #1 or 2 (uh or vh respectively)
-        self.verbose = verbose
-        if self.verbose: log.critical('Kinematic Viscosity: Building structure')
-        self.geo_structure = self.build_geo_structure() #list of Sparse matrices
-        self.operator_matrix = Sparse(self.n,self.n+self.boundary_len)
-        self.stage_heights = num.zeros((self.n,1),num.float)
-        self.stage_heights_scaling = Sparse(self.n,self.n)
-        self.apply_triangle_areas = triangle_areas
-        self.boundary_vector = num.zeros((self.n,),num.float)
-        
-        self.triangle_areas = Sparse(self.n,self.n)
-        for i in range(self.n):
-            self.triangle_areas[i,i] = 1.0 / self.mesh.areas[i]
-        if self.verbose: log.critical('Kinematic Viscosity: Initialisation Done')
-    
-    ##
-    # Enumerate the boundary conditions in some way
-    def enumerate_boundary(self):
-        #Just enumerate by the dictionary.keys() list order
-        #NOTE If this is changed, need to fix build_boundary_vector() method
-        enumeration = {}
-        keys = self.boundary.keys()
-        for k,key in enumerate(keys):
-            enumeration[key] = k
-        return enumeration
-    
-    #Create the matrix structure for the geometric component of the operator
-    # (constant across all applications of the operator)
-    def build_geo_structure(self):
-        A = [Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len)]
-        for i in range(self.n):
-            this_centroid_coords = self.mesh.centroid_coordinates[i,:]
-            for edge in range(3):
-                j = self.mesh.neighbours[i,edge]
-                edge_length = self.mesh.edgelengths[i,edge]
-                dist = 0.0
-                index = j
-                if j==-1: #Boundary edge
-                    boundary_coords = self.mesh.get_edge_midpoint_coordinate(i,edge)
-                    dist = num.linalg.norm(this_centroid_coords - boundary_coords)
-                    index = self.n + self.boundary_enum[(i,edge)]
-                else: #Interior edge
-                    other_centroid_coords = self.mesh.centroid_coordinates[j,:]
-                    dist = num.linalg.norm(this_centroid_coords - other_centroid_coords)
-                A[edge][i,i] = -edge_length / dist
-                A[edge][i,index] = edge_length / dist
-        return A
-    
-    def set_qty_considered(self,qty):
-        if qty == 1 or qty == 'u':
-            self.qty_considered = 1
-        elif qty == 2 or qty == 'v':
-            self.qty_considered = 2
-        else: #Raise an exception
-            msg = "Incorrect input qty"
-            assert 0==1, msg
-    
-    def apply_stage_heights(self,h):
-        msg = "(KV_Operator.apply_stage_heights) h vector has incorrect length"
-        assert h.size == self.n, msg
-        
-        #h is a vector of self.n stage heights (1 per triangle)
-        #scaled_geo_structure = num.zeros((3,self.n,self.n+self.boundary_len),num.float)
-        scaled_geo_structure = [Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len), Sparse(self.n,self.n+self.boundary_len)]
-        for i in range(self.n):
-            h_i = h[i]
-            for edge in range(3):
-                j = self.mesh.neighbours[i,edge]
-                h_j = 0.0
-                index = j
-                if j==-1: #Boundary edge
-                    h_j = self.boundary[(i,edge)].evaluate()[0]
-                    index = self.n + self.boundary_enum[(i,edge)]
-                else: #Interior edge
-                    h_j = h[j]
-                scaled_geo_structure[edge][i,i] = 0.5*(h_i+h_j)*self.geo_structure[edge][i,i]
-                scaled_geo_structure[edge][i,index] = 0.5*(h_i+h_j)*self.geo_structure[edge][i,index]
-        self.operator_matrix = Sparse_CSR(scaled_geo_structure[0]+scaled_geo_structure[1]+scaled_geo_structure[2])
-        #self.operator_matrix = scaled_geo_structure[0]+scaled_geo_structure[1]+scaled_geo_structure[2]
-        self.stage_heights = h
-        #Set up the scaling matrix
-        for i in range(self.n):
-            if h[i] == 0.0:
-                self.stage_heights_scaling[i,i] = 0.0
-            else:
-                self.stage_heights_scaling[i,i] = 1.0/h[i]
-        return self.operator_matrix
-    
-    def build_boundary_vector(self):
-        #Require self.qty_considered already set and stage heights applied
-        X = num.zeros((self.n+self.boundary_len,),num.float)
-        #TODO Can vectorise?
-        for i,key in enumerate(self.boundary.keys()):
-            quantities = self.boundary[key].evaluate()
-            h_i = quantities[0]
-            if h_i == 0.0:
-                X[self.n+i] = 0.0
-            else:
-                X[self.n+i] = quantities[self.qty_considered] / h_i
-        self.boundary_vector = self.operator_matrix * X
-        #Tidy up
-        if self.apply_triangle_areas:
-            self.boundary_vector = self.triangle_areas * self.boundary_vector
-
-        return self.boundary_vector
-    
-    def apply_kv_operator(self,V,qty_considered=None,include_boundary=True):
-        msg = "(KV_Operator.apply_vector) V vector has incorrect dimensions"
-        assert V.shape == (self.n,) or V.shape == (self.n,1), msg
-        
-        if not qty_considered==None:
-            print "Quantity Considered changed!"
-            self.set_qty_considered(qty_considered)
-            self.build_boundary_vector()
-        
-        D = num.zeros((self.n,),num.float)
-        #Change (uh) to (u), (vh) to (v)
-        X = self.stage_heights_scaling * V
-        
-        if self.apply_triangle_areas:
-            X = self.triangle_areas * X
-        
-        if include_boundary:
-            #TODO Tidy this up
-            #A = self.operator_matrix.todense()
-            #D = num.mat(A[:,:self.n]) * num.mat(X.reshape(self.n,1)) + self.boundary_vector
-            A = self.operator_matrix
-            D = A * num.concatenate((X, num.zeros(self.boundary_len,num.float))) + self.boundary_vector
-        else:
-            A = self.operator_matrix
-            D = A * num.concatenate((X, num.zeros(self.boundary_len,num.float)))
-            #A = self.operator_matrix.todense()
-            #D = num.mat(A[:,:self.n]) * num.mat(X.reshape(self.n,1))
-
-        return num.array(D).reshape(self.n,)
-    
-    def __mul__(self,other):
-        try:
-            B = num.array(other)
-        except:
-            msg = "Trying to multiply the Kinematic Viscosity Operator against a non-numeric type"
-            raise msg
-        
-        if len(B.shape) == 0:
-            #Scalar
-            R = B*self
-        elif len(B.shape) == 1 or (len(B.shape) == 2 and B.shape[1]==1):
-            #Vector
-            #include_boundary=False is this is *only* used for cg_solve()
-            R = self.apply_kv_operator(other,include_boundary=False)
-        else:
-            raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
-        return R
-    
-    def __rmul__(self, other):
-        #Right multiply with scalar
-        try:
-            other = float(other)
-        except:
-            msg = 'Sparse matrix can only "right-multiply" onto a scalar'
-            raise TypeError, msg
-        else:
-            new = self.operator_matrix * new
-        return new
-    
-    def cg_solve(self,B,qty_to_consider=None):
-        if len(B.shape) == 1:
-            return self.cg_solve_vector(B,qty_to_consider)
-        elif len(B.shape) == 2:
-            return self.cg_solve_matrix(B)
-        else:
-            raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
-    
-    def cg_solve_matrix(self,B):
-        assert B.shape[1] < 3, "Input matrix has too many columns (max 2)"
-        
-        X = num.zeros(B.shape,num.float)
-        for i in range(B.shape[1]):
-            X[:,i] = self.cg_solve_vector(B[:,i],i+1) #assuming B columns are (uh) and (vh)
-        return X
-    
-    def cg_solve_vector(self,b,qty_to_consider=None):
-        if not qty_to_consider==None:
-            self.set_qty_considered(qty_to_consider)
-        self.build_boundary_vector()
-        if self.verbose :
-            iprint = 1
-        else:
-            iprint = 0
-        x = conjugate_gradient(self,b - self.boundary_vector, iprint=iprint) #Call the ANUGA conjugate gradient utility
-        return x
-
-################################################################################
-
-if __name__ == "__main__":
-    #Set up a rectangular grid and check minimum/maximum principle
-
-    import time
-    t0 = time.time()
-    points,elements,boundary = mesh_factory.rectangular_cross(10, 15)
-    mesh = neighbour_mesh.Mesh(points,elements)
-    n = len(mesh)
-    boundary_map = {}
-    for (vol_id,edge) in boundary.keys():
-        x = mesh.get_edge_midpoint_coordinate(vol_id,edge)[0]
-        y = mesh.get_edge_midpoint_coordinate(vol_id,edge)[1]
-        #Define boundary conditions: u=x^2-y^2, v=x-y+2
-        boundary_map[(vol_id,edge)] = Dirichlet_boundary([1,x*x-y*y,x-y+2])
-    #Now define an operator
-    domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
-    KV_Operator = Kinematic_Viscosity_Operator(domain, verbose = True)
-    
-    #Set up the operator and RHS
-    h = num.ones((n,),num.float)
-    KV_Operator.apply_stage_heights(h)
-    rhs = num.zeros((n,2),num.float) #we will solve for u and v
-    
-    #Solve
-    velocities = KV_Operator.cg_solve(rhs)
-    #Alternatively, velocities = [u v] where:
-    u = KV_Operator.cg_solve(num.zeros((n,),num.float),'u')
-    v = KV_Operator.cg_solve(num.zeros((n,),num.float),'v')
-    vel_min = velocities.min(axis=0)
-    vel_max = velocities.max(axis=0)
-    
-    #Check a minimum/maximum principle
-    boundary_velocities = num.zeros((len(boundary_map),2),num.float)
-    for i,key in enumerate(boundary_map.keys()):
-        quantity = boundary_map[key].evaluate()
-        boundary_velocities[i,:] = quantity[1:]
-    boundary_vel_min = boundary_velocities.min(axis=0)
-    boundary_vel_max = boundary_velocities.max(axis=0)
-
-
-    print "Time ",time.time() - t0
-    print "Kinematic Viscosity Min/Max Principle Test"
-    print "Minimum Principle?", (vel_min > boundary_vel_min).all()
-    print "Maximum Principle?", (boundary_vel_max > vel_max).all()
Index: /trunk/anuga_work/development/2010-projects/kv/test_all.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/test_all.py	(revision 8052)
+++ /trunk/anuga_work/development/2010-projects/kv/test_all.py	(revision 8052)
@@ -0,0 +1,83 @@
+"""Regression testing framework
+This module will search for scripts in the same directory named
+test_*.py.  Each such script should be a test suite that tests a
+module through PyUnit. This script will aggregate all
+found test suites into one big test suite and run them all at once.
+"""
+
+# Author: Mark Pilgrim
+# Modified by Ole Nielsen
+
+import unittest
+import os
+
+
+#List files that should be excluded from the testing process.
+#E.g. if they are known to fail and under development
+exclude = [] 
+
+
+def get_test_files(path):
+
+    import sys
+
+    files = os.listdir(path)
+
+    #Check sub directories
+    test_files = []
+    for file in files:
+        if os.path.isdir(file):
+            sys.path.append(file)
+            #print 'Recursing into', file
+            test_files += get_test_files(path + os.sep + file)
+        elif file[:5] == 'test_' and file[-2:] == 'py':
+            #print 'Appending', file
+            test_files.append(file)
+        else:
+            pass
+    return test_files
+
+
+
+def regressionTest():
+    import sys, os, re, unittest
+    path = os.path.split(sys.argv[0])[0] or os.getcwd()
+
+
+    files = get_test_files(path)
+
+
+
+    #test = re.compile('^test_[\w]*.py$', re.IGNORECASE)
+    #files = filter(test.search, files)
+
+
+    try:
+        files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
+    except:
+        files.remove('test_all.py')
+
+    print 'Testing:'
+    for file in files:
+        print '  ' + file
+
+    if globals().has_key('exclude'):
+        for file in exclude:
+            files.remove(file)
+            print 'WARNING: File '+ file + ' excluded from testing'
+
+
+    filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    #print "files",files 
+    moduleNames = map(filenameToModuleName, files)
+    #print "moduleNames",moduleNames
+    modules = map(__import__, moduleNames)
+    load = unittest.defaultTestLoader.loadTestsFromModule
+    return unittest.TestSuite(map(load, modules))
+
+################################################################################
+
+if __name__ == '__main__':   
+    suite = regressionTest()
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
Index: /trunk/anuga_work/development/2010-projects/kv/test_kinematic_viscosity.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/test_kinematic_viscosity.py	(revision 8051)
+++ /trunk/anuga_work/development/2010-projects/kv/test_kinematic_viscosity.py	(revision 8052)
@@ -1,4 +1,4 @@
-from anuga.abstract_2d_finite_volumes.domain import Domain
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
+from anuga import Dirichlet_boundary
 from kinematic_viscosity import Kinematic_Viscosity_Operator
 import numpy as num
@@ -21,5 +21,5 @@
 		boundary_map[(0,1)] = Dirichlet_boundary([1,1,2])
 		boundary_map[(0,2)] = Dirichlet_boundary([1,1,0])
-		domain = Domain(source=points,triangles=elements,boundary=boundary_map)
+		domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
 		return Kinematic_Viscosity_Operator(domain)
 
@@ -33,5 +33,5 @@
 		boundary_map[(1,0)] = Dirichlet_boundary([1,1,0])
 		boundary_map[(1,2)] = Dirichlet_boundary([1,2,1])
-		domain = Domain(source=points,triangles=elements,boundary=boundary_map)
+		domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
 		return Kinematic_Viscosity_Operator(domain)
 
@@ -146,8 +146,9 @@
 		U_mod[0,:] = U_new
 		assert num.allclose(U_new - operator1.dt * 2 * num.mat(A)*num.mat(U_mod), U[0,:])
+
 ################################################################################
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Kinematic_Viscosity, 'test')
+    suite = unittest.makeSuite(Test_Kinematic_Viscosity, 'test', verbose=True)
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: unk/anuga_work/development/2010-projects/kv/test_kv_solver_new.py
===================================================================
--- /trunk/anuga_work/development/2010-projects/kv/test_kv_solver_new.py	(revision 8051)
+++ 	(revision )
@@ -1,153 +1,0 @@
-from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
-from kv_solver_new import Kinematic_Viscosity_Operator
-import numpy as num
-from math import sqrt
-import unittest
-
-class Test_kv_solver_new(unittest.TestCase):
-    def setUp(self):
-        pass
-
-    def tearDown(self):
-        pass
-    
-    #First test operator class (1 triangle)
-    def operator1(self):
-        points = num.array([[0.0,0.0],[1.0,0.0],[0.0,1.0]])
-        elements = num.array([[0,1,2]])
-        boundary_map = {}
-        boundary_map[(0,0)] = Dirichlet_boundary([1,2,1])
-        boundary_map[(0,1)] = Dirichlet_boundary([1,1,2])
-        boundary_map[(0,2)] = Dirichlet_boundary([1,1,0])
-        #We get boundary_may.keys() = [(0, 1), (0, 0), (0, 2)]
-        domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
-        return Kinematic_Viscosity_Operator(domain)
-  
-    #Second test operator class (2 triangles)
-    def operator2(self):
-        points = num.array([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0]])
-        elements = num.array([[0,1,3],[1,2,3]])
-        boundary_map = {}
-        boundary_map[(0,1)] = Dirichlet_boundary([1,1,2])
-        boundary_map[(0,2)] = Dirichlet_boundary([1,2,2])
-        boundary_map[(1,0)] = Dirichlet_boundary([1,1,0])
-        boundary_map[(1,2)] = Dirichlet_boundary([1,2,1])
-        #We get boundary_map.keys() = [(0, 1), (1, 2), (1, 0), (0, 2)]
-        domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
-        return Kinematic_Viscosity_Operator(domain)
-    
-    def test_build_index_boundary_map(self):
-        operator1 = self.operator1()
-        boundary_map_index = operator1.boundary_enum
-        assert boundary_map_index[(0,1)] == 0
-        assert boundary_map_index[(0,0)] == 1
-        assert boundary_map_index[(0,2)] == 2
-
-        operator2 = self.operator2()
-        boundary_map_index = operator2.boundary_enum
-        assert boundary_map_index[(0,1)] == 0
-        assert boundary_map_index[(1,2)] == 1
-        assert boundary_map_index[(1,0)] == 2
-        assert boundary_map_index[(0,2)] == 3
-
-    def test_kv_system_geometry(self):
-        operator1 = self.operator1()
-        operator2 = self.operator2()
-
-        kv_geo_structure = operator1.geo_structure
-        assert num.allclose(kv_geo_structure[0].todense(), num.array([[-6.0,0.0,6.0,0.0]]))
-        assert num.allclose(kv_geo_structure[1].todense(), \
-                            num.array([[-6.0/sqrt(5.0),6.0/sqrt(5.0),0.0,0.0]]))
-        assert num.allclose(kv_geo_structure[2].todense(),\
-                            num.array([[-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0)]]))
-
-        kv_geo_structure = operator2.geo_structure
-        assert num.allclose(kv_geo_structure[0].todense(), num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],\
-                                                                [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]]))
-        assert num.allclose(kv_geo_structure[1].todense(), num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\
-                                                                [3.0,-3.0,0.0,0.0,0.0,0.0]]))
-        assert num.allclose(kv_geo_structure[2].todense(), num.array([[-6.0/sqrt(5.0),0.0,0.0,0.0,0.0,6.0/sqrt(5.0)],\
-                                                                [0.0, -6.0/sqrt(5.0), 0.0, 6.0/sqrt(5.0), 0.0, 0.0]]))
-    
-    def test_apply_heights(self):
-        operator1 = self.operator1()
-        operator2 = self.operator2()
-        
-        h = num.array([1.0])
-        A = operator1.apply_stage_heights(h)
-        assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)]))
-        
-        h = num.array([2.0])
-        A = operator1.apply_stage_heights(h)
-        assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)]))
-        
-        h = num.array([1.0, 1.0])
-        A = operator2.apply_stage_heights(h)
-        #From test_kv_system_geometry
-        A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],
-                        [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]])
-        A1 = num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\
-                        [3.0,-3.0,0.0,0.0,0.0,0.0]])
-        A2 = num.array([[-6.0/sqrt(5.0),0.0,0.0,0.0,0.0,6.0/sqrt(5.0)],\
-                        [0.0, -6.0/sqrt(5.0), 0.0, 6.0/sqrt(5.0), 0.0, 0.0]])
-        assert num.allclose(A.todense(), A0+A1+A2)
-        
-        h = num.array([2.0, 1.0])
-        A = operator2.apply_stage_heights(h)
-        assert num.allclose(A.todense()[0,:], 1.5*A0[0,:]+1.5*A1[0,:]+1.5*A2[0,:])
-        assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:])
-        
-        h = num.array([-2.0, -2.0])
-        A = operator2.apply_stage_heights(h)
-        assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:])
-        assert num.allclose(A.todense()[1,:], -0.5*A0[1,:]-2*A1[1,:]-0.5*A2[1,:])
-    
-    def test_apply_kv_operator(self):
-        operator1 = self.operator1()
-        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
-        operator1.build_boundary_vector()
-        V1 = num.array([2.0]) #(uh)=2
-        V2 = num.array([2.0]) #(vh)=2
-        A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)])
-        U1 = num.array([[2.0],[1.0],[2.0],[1.0]])
-        U2 = num.array([[2.0],[2.0],[1.0],[0.0]])
-        
-        operator1.set_qty_considered('u')
-        operator1.build_boundary_vector()
-        D1 = operator1.apply_kv_operator(V1)
-        assert num.allclose(D1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))
-        
-        operator1.set_qty_considered(2)
-        operator1.build_boundary_vector()
-        D2 = operator1.apply_kv_operator(V2)
-        assert num.allclose(D2, 2*num.array(num.mat(A)*num.mat(U2)).reshape(1,))
-    
-    def test_mul(self):
-        operator1 = self.operator1()
-        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
-        operator1.set_qty_considered(1)
-        operator1.build_boundary_vector()
-        A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)])
-        V1 = num.array([2.0]) #(uh)=2
-        U1 = num.array([[2.0],[0.0],[0.0],[0.0]])
-        assert num.allclose(operator1 * V1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))
-
-    def test_cg_solve(self):
-        #cf self.test_mul()
-        operator1 = self.operator1()
-        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
-        operator1.set_qty_considered('u')
-        V = num.array([2.0]) #h=1, (uh)=2
-        A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)])
-        U = num.array([[2.0,2.0],[1.0,2.0],[2.0,1.0],[1.0,0.0]])
-        test = 2*num.mat(A)*num.mat(U[:,0].reshape(4,1))
-        X = operator1.cg_solve(num.array(test).reshape(1,))
-        assert num.allclose(V, X)
-
-################################################################################
-
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_kv_solver_new, 'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
