Index: trunk/anuga_work/development/2010-projects/kv/Report.tex
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/Report.tex	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/Report.tex	(revision 8051)
@@ -0,0 +1,95 @@
+\documentclass[a4paper]{article}
+\usepackage[margin=1in]{geometry}
+\usepackage{amsfonts}
+\usepackage{amsmath}
+\usepackage{amssymb}
+
+\newcommand{\divergence}{\operatorname{div}}
+
+\begin{document}
+\title{ANUGA Kinematic Viscosity Module}
+\author{Lindon Roberts}
+\date{\today}
+\maketitle
+
+\section{Kinematic Viscosity Equation}
+The equation we are interested in are the 2-dimensional parabolic equation given by
+\begin{align*}
+\frac{\partial u}{\partial t} &= (hu_{x})_{x} + (hu_{y})_{y} \\
+\frac{\partial v}{\partial t} &= (hv_{x})_{x} + (hv_{y})_{y}
+\end{align*}
+where $h$ is the stage height of the water, and $u$ and $v$ are the $x$- and $y$- velocities of the water. 
+
+\subsection{Elliptic Approximation}
+Given some discretisation of the 2d domain into triangles $\{T_{i}\}_{i=1}^{n}$, and letting $\partial T_{i}$ be the set of indices $j$ such that $T_{i}$ and $T_{j}$ share an edge, we approximate as follows
+\begin{align*}
+\int_{T_{i}}\divergence(u_{x}h,u_{y}h) dA &= \int_{\partial T_{i}}(u_{x}h,u_{y}h)\cdot\nu dS \\
+&= \sum_{j\in\partial T_{i}}\ell_{ij}((u_{x}h,u_{y}h)\cdot\nu) \\
+&\approx \sum_{j\in\partial T_{i}}\ell_{ij}\left(\frac{h_{i}+h_{j}}{2}\right)(Du\cdot\nu) \\
+&\approx \sum_{j\in\partial T_{i}}\ell_{ij}\left(\frac{h_{i}+h_{j}}{2}\right)\left(\frac{u_{i}-u_{j}}{d_{ij}}\right)
+\end{align*}
+where $\ell_{ij}$ is the length of the edge between $T_{i}$ and $T_{j}$, and $d_{ij}$ is the (straight-line) distance between the centroids of $T_{i}$ and $T_{j}$, and $\nu$ is the outward normal. 
+
+We now have the option of whether to approximate the divergence at the centroid as the integral of the divergence over the whole triangle, or whether it should be the integral average of the divergence over the triangle, as given by
+\[\frac{1}{|T_{i}|}\int_{T_{i}}\divergence(u_{x}h,u_{y}h) dA \approx \sum_{j\in\partial T_{i}}\frac{\ell_{ij}}{|T_{i}|}\left(\frac{h_{i}+h_{j}}{2}\right)\left(\frac{u_{i}-u_{j}}{d_{ij}}\right)\]
+Which approximation we use is governed by the boolean field \verb|apply_triangle_areas|, which is set in the constructor (default value True). The module has the ability to solve the elliptic equation $\divergence(u_{x}h,u_{y}h)=f$ for arbitrary right-hand side $f$ also. This is implemented in the \verb|cg_solve| method.
+
+\subsection{Parabolic Solver}
+To solve the full parabolic equation, we use Implicit Euler for timestepping. We will see below that the elliptic operator $A[h,u]\approx\divergence(u_{x}h,u_{y}h)$ can be written in the form $A[h,u]=B(h)u + c(h)$ where $B(h)$ is a $n\times n$ matrix and $c(h)$ is a length $n$ vector. Therefore our timestepping is given by
+\[\frac{u^{n+1}-u^{n}}{\Delta t} = A[h^{n},u^{n+1}] = B(h^{n})u^{n+1} + c(h^{n})\]
+which rearranges to give
+\[(I-\Delta t B(h^{n}))u^{n+1} = u^{n}+c(h^{n})\Delta t\]
+This is the equation we solve in the \verb|parabolic_solve| method.
+
+\section{Setting up the matrix problem}
+Since the elliptic interaction is on all the $n$ interior nodes by the interior nodes \emph{and} the $m$ boundary nodes, we actually get a $n\times(n+m)$ sparse matrix of interactions. Along the $i$th row, we have four nonzero entries. We define
+\[g_{ij} = -\frac{\ell_{ij}}{d_{ij}} = g_{i}[edge]\]
+where $edge=0,1,2$ (the edge opposite the vertex numbered 0, 1 or 2), and $T_{j}$ for $j\in\{0,\ldots,n+m-1\}$ is the triangle which neighbours $T_{i}$ along $T_{i}$'s edge number $edge$. If this edge is a boundary, then $j\geq n$ and we are interested in the boundary point at the midpoint of the edge. We enumerate the boundary edges $(i,edge)$ by first sorting by node $i$, then by $edge$. 
+
+We store the values $g_{i}[edge]$ in a $n\times3$ matrix \verb|geo_structure_values| in which the $(i,edge)$ entry is $g_{i}[edge]$. We also store another $n\times3$ matrix \verb|geo_structure_indices| in which the $(i,edge)$ entry is $j$. These matrices are computed during the initialisation process using the C extension module, as they are geometry-dependent only. We use these matrices and the stage heights $h$ to define the elliptic operator $A[h,\cdot]$.
+
+\subsection{Applying Stage Heights}
+We are given a length $n$ vector of stage heights $h$, from which we need to form the elliptic operator. As stated above, the operator is actually a $n\times(n+m)$ matrix, where we have the off-diagonal entries (in the appropriate places, as given by the matrix \verb|geo_structure_values|, which maps $(i,edge)$ to $j$)
+\[A_{ij} = \left(\frac{h_{i}+h_{j}}{2}\right)g_{ij}\]
+and diagonal entries
+\[A_{ii} = -\sum_{j\in\partial T_{i}}A_{ij}\]
+We of course need to evaluate the stage height on the boundary edges, which we do using the evaluate method associated with each boundary class (the collection of boundary classes is stored in a dictionary with key $(i,edge)$). All the height values (the length $n$ vector $h$ and a length $m$ vector of boundary values) is read into the C extension, along with the geometric structure matrices. The output from this is a tuple $(data,colind,rowptr)$, which are the three arrays corresponding to the CSR Sparse format for matrices (implemented in the ANUGA utilities).
+
+Given that our operator will receive the $uh$ and $vh$ values as inputs, rather than $u$ and $v$, we also create a matrix \verb|stage_height_scaling| which is a diagonal matrix with entries $1/h[i]$ ($h_{i}>0$) or 0 if $h_{i}=0$ (in this case, the entire interaction should be zero, so this works). Thus applying this matrix to $uh$ and $vh$ will give us $u$ and $v$ respectively. Thus we construct this matrix during this process.
+
+We now notice that we can only apply this operator matrix to a length $(n+m)$ vector of $u$ or $v$ values. The first $n$ entries of this vector will be our input, but the last $m$ entries are the values of $u$ and $v$ at the boundary, which we can evaluate already. If we write $A=[B|B']$ where $B$ is a $n\times n$ matrix and $B'$ is a $n\times m$ matrix, then the result of multiplyng the full vector of $u$ (or $v$) by $A$ is the same as $Bu+B'b$ where $b$ is the $m$-vector of boundary values of $u$. Therefore we compute the length $n$ vector $c(h)=B'b$, stored as \verb|boundary_vector|, and computed in \verb|build_boundary_vector|. Since we need $c(h)$ for $u$ and $v$, we store this as a $n\times 2$ matrix. In the solving, we will only use one of these columns depending on whether or not we are solving for $u$ or $v$. This is achieved by having the \verb|qty_considered| value equal to 1 ($u$) or 2 ($v$).
+
+\section{Matrix Multiplication and Solving}
+\subsection{Elliptic Multiplication}
+For the elliptic multiplication, we take in a length $n$ vector (either $uh$ or $vh$). We first need to change $uh$ to $u$ (or $vh$ to $v$), which we do by multiplying by the \verb|stage_heights_scaling| matrix (see above). If we are using the integral average formula rather than the integral formula (see the initial equations), we need to scale down each entry by the area of each triangle. We then extend this vector by zeros to make it of length $n+m$ (so we can multiply by $A$, but have the same effect as multiplying by $B$ only), and then multiply by $A$. We have now computed $B(h)u$. In general, if we just want to evaluate $A[h,u]$ (when \verb|include_boundary| is True), then we need to add $c(h)$ to the result, so we do this. We do not want to do this when we are solving (see below), which is why we set it as an argument.
+
+We then scale the result back from $u$ to $uh$, in order to get consistency in the solver.
+
+\subsection{Parabolic Multiplication}
+This is almost identical to the elliptic case, except we calculate $(I-\Delta t B(h))u$ (note we do not include the boundary). We do not include the boundary term $c(h)$ because this method is only ever called during the parabolic solve. In the implementation, we assume that we get $u$ (not $uh$) as the input. The scaling is done in the solver wrapper.
+
+\subsection{Solving}
+The solving is done with the ANUGA conjugate gradient solver. In the elliptic case, we read in a length $n$ or $n\times 2$ right-hand side. If it is a vector, then we need to ensure we are considering the right quantity (see above, \verb|qty_considered|). If it is a matrix, then solve for $uh$ and $vh$ using the two columns respectively. This is split up using the various \verb|cg_solve| methods. We then solve the equation
+\[B(h)u = f - c(h)\]
+for $u$. This equivalently solves $A[h,u]=f$ for $u$ (or $v$) and return either a vector of $uh$ values or a matrix of $uh$ and $vh$ values (depending on the input). The elliptic solver is \verb|cg_solve|.
+
+For the parabolic solver, we assume that we are given both $(uh)^{n}$ and $(vh)^{n}$ as input. We then scale down to $u^{n}$ and $v^{n}$ (given that the most recent call of \verb|apply_stage_heights| would have been for $h^{n}$), and solve the equation
+\[(I-\Delta t B(h^{n}))u^{n+1} = u^{n}+c(h^{n})\Delta t\]
+
+Note that the elliptic solver returns $uh$, but the parabolic solver only returns $u^{n+1}$ and $v^{n+1}$ (the non-scaled values), since we do not have the values of $h^{n+1}$ available.
+
+\section{Next Steps}
+The module is not yet complete. It functions as a stand-alone module, but is not yet fully integrated into the ANUGA code. The primary things that need to be addressed are
+\begin{itemize}
+\item Accessing the value of $\Delta t$ from the domain input.
+\item Accessing the boundary values of $h$, $uh$ and $vh$ from the domain class. This may need to involve some fixing of the enumeration of the boundary sides used in the body of the code.
+\item Moving from accepting as input stage heights to just heights (stage heights + elevation).
+\item Simplify the code to take in a matrix with the $h$, $uh$ and $vh$ values in one go, apply stage heights and then perform a parabolic solve.
+\item Including the module as part of the evolve process for the ANUGA program.
+\item Setting up the compiler for Linux (all testing and development done on Windows).
+\end{itemize}
+
+\section{More Information}
+For some more information, including sample usage, go to the Readme, or the test file.
+
+\end{document}
Index: trunk/anuga_work/development/2010-projects/kv/compile.py
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/compile.py	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/compile.py	(revision 8051)
@@ -0,0 +1,452 @@
+"""compile.py - compile Python C-extension
+
+   Commandline usage: 
+     python compile.py <filename>
+
+   Usage from within Python:
+     import compile
+     compile.compile(<filename>,..)
+
+   Ole Nielsen, Duncan Gray Oct 2001      
+"""     
+
+#NumPy ------------------------------------
+# Something like these lines recommended in "Converting from NUMARRAY to NUMPY"
+import numpy
+I_dirs = '-I"%s" ' % numpy.get_include()
+#NumPy ------------------------------------
+
+# FIXME (Ole): Although this script says it works with a range of compilers,
+# it has only really been used with gcc.
+
+import os, string, sys, types
+
+separation_line = '---------------------------------------'      
+ 
+def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
+  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
+  
+     Compile FN(s) using compiler CC (e.g. mpicc), 
+     Loader LD and shared flag SFLAG.
+     If CC is absent use default compiler dependent on platform
+     if LD is absent CC is used.
+     if SFLAG is absent platform default is used
+     FNs can be either one filename or a list of filenames
+     In the latter case, the first will be used to name so file.
+  """
+  
+  # Input check
+  #
+  assert not FNs is None, 'No filename provided'
+
+  if not type(FNs) == types.ListType:
+    FNs = [FNs]
+
+
+  libext = 'so' #Default extension (Unix)
+  libs = ''
+  version = sys.version[:3]
+  
+  # Determine platform and compiler
+  #
+  if sys.platform == 'sunos5':  #Solaris
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+      
+  elif sys.platform == 'osf1V5':  #Compaq AlphaServer
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'linux2':  #Linux
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'    
+      
+  elif sys.platform == 'darwin':  #Mac OS X:
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'bundle -flat_namespace -undefined suppress'
+
+  elif sys.platform == 'cygwin':  #Cygwin (compilation same as linux)
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc'
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+    libext = 'dll'
+    libs = '/lib/python%s/config/libpython%s.dll.a' %(version,version)
+      
+  elif sys.platform == 'win32':  #Windows
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'gcc.exe' #Some systems require this (a security measure?) 
+    if LD:
+      loader = LD
+    else:  
+      loader = compiler
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'shared'
+     
+    # As of python2.5, .pyd is the extension for python extension
+    # modules.
+    if sys.version_info[0:2] >= (2, 5):
+      libext = 'pyd'
+    else:
+      libext = 'dll'
+
+    libs, is_found = set_python_dll_path()
+      
+  else:
+    if verbose: print "Unrecognised platform %s - revert to default"\
+                %sys.platform
+    
+    if CC:
+      compiler = CC
+    else:  
+      compiler = 'cc'
+    if LD:
+      loader = LD
+    else:  
+      loader = 'ld'
+    if SFLAG:
+      sharedflag = SFLAG
+    else:  
+      sharedflag = 'G'
+
+
+  # Verify that compiler can be executed
+  print 'Compiler: %s, version ' %compiler,
+  sys.stdout.flush()
+  s = '%s -dumpversion' %(compiler)
+  err = os.system(s)
+  print
+  
+  if err != 0:
+      msg = 'Unable to execute compiler: %s. ' %compiler
+      msg += 'Make sure it is available on the system path.\n'
+      msg += 'One way to check this is to run %s on ' %compiler
+      msg += 'the commandline.'
+      raise Exception, msg    
+
+  
+  # Find location of include files
+  #
+  if sys.platform == 'win32':  #Windows
+    python_include = os.path.join(sys.exec_prefix, 'include')    
+  else:  
+    python_include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
+                                  'python' + version)
+
+  # Check existence of Python.h
+  #
+  headerfile = python_include + os.sep + 'Python.h'
+  try:
+    open(headerfile, 'r')
+  except:
+    raise """Did not find Python header file %s.
+    Make sure files for Python C-extensions are installed. 
+    In debian linux, for example, you need to install a
+    package called something like python2.3-dev""" %headerfile
+
+
+
+  # Add Python path + utilities to includelist (see ticket:31)
+  # Assume there is only one 'utilities' dir under path dirs
+  
+  utilities_include_dir = None
+  for pathdir in sys.path:
+
+    utilities_include_dir = pathdir + os.sep + 'utilities'
+    #print pathdir
+    #print utilities_include_dir
+    try:
+      os.stat(utilities_include_dir)
+    except OSError:
+      pass
+    else:
+      #print 'Found %s to be used as include dir' %utilities_include_dir
+      break
+
+  # This is hacky since it
+  # assumes the location of the compile_all that determines buildroot
+  try:
+    utilities_include_dir = buildroot + os.sep + "source" + os.sep + "anuga" \
+                            + os.sep + 'utilities'
+  except:
+    # This will make compile work locally
+    utilities_include_dir = '.'
+    utilities_include_dir = '../utilities'
+
+
+    
+  try:
+    os.stat(utilities_include_dir)
+  except OSError: 
+    utilities_include_dir = buildroot + os.sep + 'utilities'
+  
+  
+  
+  # Check filename(s)
+  #
+  object_files = ''
+  for FN in FNs:        
+    root, ext = os.path.splitext(FN)
+    if ext == '':
+      FN = FN + '.c'
+    elif ext.lower() != '.c':
+      raise Exception, "Unrecognised extension: " + FN
+    
+    try:
+      open(FN, 'r')
+    except:
+      raise Exception, "Could not open: " + FN
+
+    if not object_files: root1 = root  #Remember first filename        
+    object_files += root + '.o '  
+  
+  
+    # Compile
+    #
+    if utilities_include_dir is None:    
+      s = '%s -c %s -I"%s" -o "%s.o" -Wall -O3'\
+          %(compiler, FN, python_include, root)
+    else:
+      if FN == "triangle.c" or FN == "mesh_engine_c_layer.c":
+        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -O3 -DTRILIBRARY=1 -DNO_TIMER=1'\
+            % (compiler, FN, I_dirs, python_include, utilities_include_dir, root)
+      elif FN == "polygon_ext.c":
+        # gcc 4.3.x is problematic with this file if '-O3' is used
+        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -Wall'\
+            %(compiler, FN, I_dirs, python_include, utilities_include_dir, root)
+      else:
+        s = '%s -c %s %s -I"%s" -I"%s" -o "%s.o" -Wall -O3'\
+            %(compiler, FN, I_dirs, python_include, utilities_include_dir, root)
+
+    if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
+      # Extra flags for 64 bit architectures.
+      #   AMD Opteron will give x86_64
+      #   Itanium will give ia64
+      #
+      # Second clause will always fail on Win32 because uname is UNIX specific
+      # but won't get past first clause
+
+      s += ' -fPIC' # Use position independent code for 64 bit archs
+      #s += ' -m64' # Used to be necessary for AMD Opteron
+      
+      
+    if verbose:
+      print s
+
+    # Doesn't work on Windows anyway  
+    #else:
+    #  s = s + ' 2> /dev/null' #Suppress errors
+    
+    try:
+      err = os.system(s)
+      if err != 0:
+          raise 'Attempting to compile %s failed - please try manually' %FN 
+    except:
+      raise 'Could not compile %s - please try manually' %FN  
+
+  
+  # Make shared library (*.so or *.dll)
+  if libs is "":
+    s = '%s -%s %s -o %s.%s -lm' %(loader, sharedflag, object_files, root1, libext)
+  else:
+    s = '%s -%s %s -o %s.%s "%s" -lm' %(loader, sharedflag, object_files, root1, libext, libs)
+    
+  if verbose:
+    print s
+
+  # Doesn't work on Windows anyway      
+  #else:
+  #  s = s + ' 2> /dev/null' #Suppress warnings
+  
+  try:  
+    err=os.system(s)
+    if err != 0:        
+        raise 'Attempting to link %s failed - please try manually' %root1     
+  except:
+    raise 'Could not link %s - please try manually' %root1
+    
+
+def can_use_C_extension(filename):
+    """Determine whether specified C-extension
+    can and should be used.
+    """
+
+    from anuga.config import use_extensions
+
+    from os.path import splitext
+
+    root, ext = splitext(filename)
+    
+    C=False
+    if use_extensions:
+        try:
+            s = 'import %s' %root
+            #print s
+            exec(s)
+        except:
+            try:
+                open(filename)
+            except:
+                msg = 'C extension %s cannot be opened' %filename
+                print msg                
+            else:    
+                print '------- Trying to compile c-extension %s' %filename
+            
+                try:
+                    compile(filename)
+                except:
+                    print 'WARNING: Could not compile C-extension %s'\
+                          %filename
+                else:
+                    try:
+                        exec('import %s' %root)
+                    except:
+                        msg = 'C extension %s seems to compile OK, '
+                        msg += 'but it can still not be imported.'
+                        raise msg
+                    else:
+                        C=True
+        else:
+            C=True
+            
+    if not C:
+        pass
+        print 'NOTICE: C-extension %s not used' %filename
+
+    return C
+
+
+def set_python_dll_path():
+  """ On windows, find which of the two usual hiding places the python
+  dll is located.
+
+  If the file can't be found, return None.
+  """
+  import sys
+  from os import access, F_OK
+  
+  version = sys.version[:3]
+  v = version.replace('.','')
+  dllfilename = 'python%s.dll' %(v)
+  libs = os.path.join(sys.exec_prefix,dllfilename)
+
+  is_found = True   
+  if access(libs,F_OK) == 0 :
+    # Sometimes python dll is hidden in %WINDIR%/system32
+    libs = os.path.join(os.environ.get('WINDIR', 'C:\WINNT'), 'system32',
+                        dllfilename)
+    if access(libs,F_OK) == 0 :
+      # could not find the dll
+      libs = os.path.join(sys.exec_prefix,dllfilename)
+      is_found = False
+  return libs, is_found
+
+def check_python_dll():
+  libs, is_found = set_python_dll_path()
+  if not is_found:
+    print "%s not found.\nPlease install.\nIt is available on the web." \
+              %(libs)
+    import sys; sys.exit()
+    
+      
+if __name__ == '__main__':
+  from os.path import splitext
+  
+  if sys.platform == 'win32':
+    check_python_dll()
+  if len(sys.argv) > 1:
+      files = sys.argv[1:]
+      for filename in files:
+          root, ext = splitext(filename)
+
+          if ext <> '.c':
+              print 'WARNING (compile.py): Skipping %s. I only compile C-files.' %filename
+      
+  else:  
+      #path = os.path.split(sys.argv[0])[0] or os.getcwd()
+      path = '.'
+      files = os.listdir(path)
+      
+      
+
+  for filename in files:
+      root, ext = splitext(filename)
+
+      if ext == '.c':
+          for x in ['.dll', '.so']:
+              try:
+                  os.remove(root + x)
+              except:
+                  pass
+
+          print separation_line
+          print 'Trying to compile c-extension %s in %s'\
+                %(filename, os.getcwd())
+          try:
+            if filename == 'triang.c': 
+              compile(['triang.c','triangle.c'])
+            elif  filename == 'mesh_engine_c_layer.c': 
+              compile(['mesh_engine_c_layer.c','triangle.c'])
+            else:
+              compile(filename)
+          except Exception, e:
+              msg = 'Could not compile C extension %s\n' %filename
+              msg += str(e)
+              raise Exception, msg
+          else:
+              print 'C extension %s OK' %filename
+          print    
+        
+
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 8051)
@@ -0,0 +1,230 @@
+from anuga.abstract_2d_finite_volumes.domain import Domain
+#from anuga.utilities.sparse import Sparse, Sparse_CSR
+from sparse import Sparse, Sparse_CSR #the new import
+from anuga.utilities.cg_solve import conjugate_gradient
+import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+import numpy as num
+import kinematic_viscosity_ext
+import anuga.utilities.log as log
+
+class Kinematic_Viscosity_Operator:
+
+	def __init__(self, domain, triangle_areas=True, verbose=False):
+		if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
+		#Expose the domain attributes
+		self.domain = domain
+		self.mesh = domain.mesh
+		self.boundary = domain.boundary
+		self.n = len(self.domain)
+		self.dt = 1e-6 #Need to set to domain.timestep
+		self.boundary_len = len(domain.boundary)
+		self.tot_len = self.n+self.boundary_len
+		self.boundary_enum = self.enumerate_boundary()
+		self.verbose = verbose
+		
+		#Geometric Information
+		if verbose: log.critical('Kinematic Viscosity: Building geometric structure')
+		self.geo_structure_indices = num.zeros((self.n,3),num.int)
+		self.geo_structure_values = num.zeros((self.n,3),num.float)
+		kinematic_viscosity_ext.build_geo_structure(self,self.n,self.tot_len)
+		self.apply_triangle_areas = triangle_areas
+		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]
+		self.triangle_areas = Sparse_CSR(self.triangle_areas)
+		
+		#Application Information
+		self.qty_considered = 1 #1 or 2 (uh or vh respectively)
+		self.operator_matrix = Sparse(self.n,self.tot_len)
+		self.operator_data = num.zeros((4*self.n,),num.float) #Sparse_CSR.data
+		self.operator_colind = num.zeros((4*self.n,),num.int) #Sparse_CSR.colind
+		self.operator_rowptr = 4*num.indices((self.n+1,))[0,:] #Sparse_CSR.rowptr (4 entries in every row, we know this already) = [0,4,8,...,4*n]
+		self.stage_heights = num.zeros((self.n,1),num.float)
+		self.stage_heights_scaling = Sparse(self.n,self.n)
+		self.boundary_vector = num.zeros((self.n,),num.float)
+		
+		self.parabolic_solve = False #Are we doing a parabolic solve at the moment?
+		
+		if verbose: log.critical('Kinematic Viscosity: Initialisation Done')
+    
+	# Enumerate the boundary conditions in some way
+	def enumerate_boundary(self):
+		#Just enumerate by the triangle number, then edge number
+		enumeration = {}
+		enum = 0
+		for i in range(self.n):
+			for edge in range(3):
+				if self.mesh.neighbours[i,edge] == -1:
+					enumeration[(i,edge)] = enum
+					enum += 1
+		return enumeration
+	
+	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
+		
+		self.operator_matrix = Sparse(self.n,self.tot_len)
+		
+		#Evaluate the boundary stage heights - use generic_domain.update_boundary? (check enumeration)
+		boundary_heights = num.zeros((self.boundary_len,1),num.float)
+		for key in self.boundary_enum.keys():
+			boundary_heights[self.boundary_enum[key]] = self.boundary[key].evaluate()[0]
+		
+		kinematic_viscosity_ext.build_operator_matrix(self,self.n,self.tot_len,h,boundary_heights)
+		self.operator_matrix = Sparse_CSR(None,self.operator_data,self.operator_colind,self.operator_rowptr,self.n,self.tot_len)
+		
+		self.stage_heights = h
+		#Set up the scaling matrix
+		data = h
+		num.putmask(data, data!=0, 1/data) #take the reciprocal of each entry unless it is zero
+		self.stage_heights_scaling = Sparse_CSR(None,data,num.arange(self.n),num.append(num.arange(self.n),self.n),self.n,self.n)
+		
+		self.build_boundary_vector()
+		
+		return self.operator_matrix
+	
+	def build_boundary_vector(self):
+        #Require stage heights applied
+		X = num.zeros((self.tot_len,2),num.float)
+		for key in self.boundary_enum.keys():
+			quantities = self.boundary[key].evaluate()
+			h_i = quantities[0]
+			if h_i == 0.0:
+				X[self.n + self.boundary_enum[key],:] = num.array([0.0, 0.0])
+			else:
+				X[self.n + self.boundary_enum[key],:] = quantities[1:] / 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 elliptic_multiply(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 qty_considered!=None:
+			print "Quantity Considered changed!"
+			self.set_qty_considered(qty_considered)
+		
+		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
+		
+		X = num.append(X,num.zeros((self.boundary_len,)),axis=0)
+		D = self.operator_matrix * X
+		
+		if include_boundary:
+			D += self.boundary_vector[:,self.qty_considered-1]
+		
+		D = self.stage_heights * D #make sure we return (uh) not (u), for the solver's benefit
+		
+		return num.array(D).reshape(self.n,)
+	
+	#parabolic_multiply(V) = identity - dt*elliptic_multiply(V)
+	#We do not need include boundary, as we will only use this method for solving (include_boundary = False)
+	#Here V is already either (u) or (v), not (uh) or (vh)
+	def parabolic_multiply(self,V):
+		msg = "(KV_Operator.parabolic_multiply) V vector has incorrect dimensions"
+		assert V.shape == (self.n,) or V.shape == (self.n,1), msg
+		
+		D = V #The return value
+		X = V #a temporary array
+		
+		#Apply triangle areas
+		if self.apply_triangle_areas:
+			X = self.triangle_areas * X
+		
+		#Multiply out
+		X = num.append(X,num.zeros((self.boundary_len,)),axis=0)
+		D = D - self.dt * (self.operator_matrix * X)
+		
+		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
+			if self.parabolic_solve:
+				R = self.parabolic_multiply(other)
+			else:
+				#include_boundary=False is this is *only* used for cg_solve()
+				R = self.elliptic_multiply(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)
+		x = conjugate_gradient(self,b - self.boundary_vector[:,self.qty_considered-1],iprint=1) #Call the ANUGA conjugate gradient utility
+		return x
+	
+	#Solve the parabolic equation to find u^(n+1), where u = u^n
+	#Here B = [uh vh] where uh,vh = n*1 column vectors
+	def parabolic_solver(self,B):
+		assert B.shape == (self.n,2), "Input matrix has incorrect dimensions"
+		
+		next_B = num.zeros((self.n,2),num.float)
+		#The equation is self.parabolic_multiply*u^(n+1) = u^n + (dt)*self.boundary_vector
+		self.parabolic_solve = True
+		#Change (uh) and (vh) to (u) and (v)
+		b = self.stage_heights_scaling * B
+
+		next_B = conjugate_gradient(self, b + (self.dt * self.boundary_vector), iprint=1)
+		
+		#Return (u) and (v) back to (uh) and (vh)
+		next_B = self.stage_heights_scaling * next_B
+		
+		self.parabolic_solve = False
+		
+		return next_B
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 8051)
@@ -0,0 +1,214 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include <math.h>
+#include <stdio.h>
+#include <unistd.h>
+#include "util_ext.h"
+
+//Rough quicksort implementation (for build_operator_matrix)
+int *quicksort(int *array, int n) {
+	int *less, *more, pivot, i, num_less, num_more, *sorted;
+	less = malloc(sizeof(int)*n);
+	more = malloc(sizeof(int)*n);
+	sorted = malloc(sizeof(int)*n);
+	if (n<2) return array;
+	pivot = array[n-1]; //choose the last element as a pivot
+	num_less = 0; num_more = 0;
+	for (i=0; i<n-1; i++) {
+		if (array[i] <= pivot) {
+			less[num_less] = array[i];
+			num_less++;
+		} else {
+			more[num_more] = array[i];
+			num_more++;
+		}
+	}
+	
+	less = quicksort(less,num_less);
+	more = quicksort(more,num_more);
+	
+	for (i=0; i<n; i++) {
+		if (i<num_less) {
+			sorted[i] = less[i];
+		} else if (i==num_less) {
+			sorted[i] = pivot;
+		} else {
+			sorted[i] = more[i-num_less-1];
+		}
+	}
+	return sorted;
+}
+
+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 i, edge, edge_counted, j;
+	double dist, this_x, this_y, other_x, other_y, edge_length;
+	edge_counted = 0;
+	for (i=0; i<n; i++) {
+		//The centroid coordinates of triangle i
+		this_x = centroids[2*i];
+		this_y = centroids[2*i+1];
+		for (edge=0; edge<3; edge++) {
+		
+			j = neighbours[3*i+edge];
+			
+			//Get the index and the coordinates of the interacting point
+			if (j == -1) {
+				geo_indices[3*i+edge] = n + edge_counted;
+				edge_counted++;
+				
+				other_x = edge_midpoints[2*(3*i+edge)];
+				other_y = edge_midpoints[2*(3*i+edge)+1];
+			} else {
+				geo_indices[3*i+edge] = j;
+				
+				other_x = centroids[2*j];
+				other_y = centroids[2*j+1];
+			}
+			
+			//Compute the interaction
+			edge_length = edgelengths[3*i+edge];
+			dist = sqrt((this_x-other_x)*(this_x-other_x) + (this_y-other_y)*(this_y-other_y));
+			geo_values[3*i+edge] = - edge_length / dist;
+		}
+	}
+	return 0;
+}
+
+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 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
+	for (i=0; i<n; i++) {
+		v_i = 0.0;
+		j[3] = i;
+		//Get the values of each interaction, and the column index at which they occur
+		for (edge=0; edge<3; edge++) {
+			j[edge] = geo_indices[3*i+edge];
+			if (j[edge]<n) { //interior
+				h_j = h[j[edge]];
+			} else { //boundary
+				h_j = boundary_heights[j[edge]-n];
+			}
+			v[edge] = -0.5*(h[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
+			v_i += 0.5*(h[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
+		}
+		//Organise the set of 4 values/indices into the data and colind arrays
+		sorted_j = quicksort((int *)j,4);
+		for (k=0; k<4; k++) { //loop through the nonzero indices
+			this_index = sorted_j[k];
+			if (this_index == i) {
+				data[4*i+k] = v_i;
+				colind[4*i+k] = i;
+			} else if (this_index == j[0]) {
+				data[4*i+k] = v[0];
+				colind[4*i+k] = j[0];
+			} else if (this_index == j[1]) {
+				data[4*i+k] = v[1];
+				colind[4*i+k] = j[1];
+			} else { //this_index == j[2]
+				data[4*i+k] = v[2];
+				colind[4*i+k] = j[2];
+			}
+		}
+	}
+	return 0;
+}
+
+/**
+* Python wrapper methods
+*/
+
+static PyObject *py_build_geo_structure(PyObject *self, PyObject *args) {
+    PyObject *kv_operator, *mesh;
+    int n, tot_len, err;
+    PyArrayObject
+            *centroid_coordinates,
+            *neighbours,
+            *edgelengths,
+            *edge_midpoint_coordinates,
+			*geo_indices,
+			*geo_values;
+    
+    //Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "Oii", &kv_operator, &n, &tot_len)) {
+        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input");
+        return NULL;
+    }
+    mesh = PyObject_GetAttrString(kv_operator, "mesh"); //kv_operator.mesh
+    if (!mesh) {
+        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not obtain mesh object from kv_operator");
+        return NULL;
+    }
+
+    //Extract parameters
+    centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
+    neighbours = get_consecutive_array(mesh,"neighbours");
+    edgelengths = get_consecutive_array(mesh,"edgelengths");
+    edge_midpoint_coordinates = get_consecutive_array(mesh,"edge_midpoint_coordinates");
+	geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
+	geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
+    
+    //Release
+    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);
+    if (err != 0) {
+        PyErr_SetString(PyExc_RuntimeError, "Could not build geo structure");
+        return NULL;
+    }
+    
+    //Release the arrays
+    Py_DECREF(centroid_coordinates);
+    Py_DECREF(neighbours);
+    Py_DECREF(edgelengths);
+    Py_DECREF(edge_midpoint_coordinates);
+	Py_DECREF(geo_indices);
+	Py_DECREF(geo_values);
+    
+    return Py_BuildValue("");
+}
+
+static PyObject *py_build_operator_matrix(PyObject *self, PyObject *args) {
+	PyObject *kv_operator;
+	int n, tot_len, err;
+	PyArrayObject
+		*h,
+		*boundary_heights,
+		*geo_indices,
+		*geo_values,
+		*_data,
+		*colind;
+	
+	//Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OiiOO", &kv_operator, &n, &tot_len, &h, &boundary_heights)) {
+        PyErr_SetString(PyExc_RuntimeError, "get_stage_height_interactions could not parse input");
+        return NULL;
+    }
+	
+	geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
+	geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
+	_data = get_consecutive_array(kv_operator,"operator_data");
+	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);
+    if (err != 0) {
+        PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
+        return NULL;
+    }
+	
+	Py_DECREF(geo_indices);
+	Py_DECREF(geo_values);
+	Py_DECREF(_data);
+	Py_DECREF(colind);
+    
+    return Py_BuildValue("");
+}
+
+static struct PyMethodDef MethodTable[] = {
+        {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
+		{"build_operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},
+        {NULL,NULL,0,NULL} // sentinel
+};
+void initkinematic_viscosity_ext(){
+    (void) Py_InitModule("kinematic_viscosity_ext", MethodTable);
+    import_array();
+}
Index: trunk/anuga_work/development/2010-projects/kv/make.bat
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/make.bat	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/make.bat	(revision 8051)
@@ -0,0 +1,2 @@
+gcc -I"C:\Python26\lib\site-packages\numpy\core\include" -I"C:\Python26\include" -Wall -c kinematic_viscosity_ext.c
+gcc -shared kinematic_viscosity_ext.o -o kinematic_viscosity_ext.pyd "C:\Windows\system32\python26.dll" -lm
Index: trunk/anuga_work/development/2010-projects/kv/profiling.py
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/profiling.py	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/profiling.py	(revision 8051)
@@ -0,0 +1,43 @@
+from anuga.abstract_2d_finite_volumes.domain import Domain
+import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory
+import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+from kinematic_viscosity import Kinematic_Viscosity_Operator
+import numpy as num
+import anuga.utilities.log as log
+
+def evolve():
+    print "Setting up"
+    grid_rows = 40
+    grid_cols = 40
+    #Set up a rectangular grid
+    points,elements,boundary = mesh_factory.rectangular_cross(grid_rows,grid_cols)
+    mesh = neighbour_mesh.Mesh(points,elements)
+    n = len(mesh)
+    print "n =",n
+    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 = Domain(source=points,triangles=elements,boundary=boundary_map)
+    KV_Operator = Kinematic_Viscosity_Operator(domain)
+    print "Applying stage heights"
+    #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
+    print "Solving"
+    #Solve
+    velocities = KV_Operator.cg_solve(rhs)
+    new_velocities = KV_Operator.parabolic_solver(velocities)
+
+#Do the profiling
+import profile, pstats
+FN = 'kv_profile.dat'
+print "Starting..."
+profile.run('evolve()', FN)
+S = pstats.Stats(FN)
+s = S.sort_stats('cumulative').print_stats(30)
Index: trunk/anuga_work/development/2010-projects/kv/readme.txt
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/readme.txt	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/readme.txt	(revision 8051)
@@ -0,0 +1,66 @@
+ANUGA Kinematic Viscosity Module
+Lindon Roberts, 2010
+
+Summary
+============================================================================
+This module applies kinematic viscosity smoothing to the velocities in a 
+domain, given by the parabolic equations
+  du/dt = (h*u_x)_x + (h*u_y)_y
+  dv/dt = (h*v_x)_x + (h*v_y)_y
+for u(x,y) and v(x,y) being water velocities, and h(x,y) is the stage height.
+For more details on the implementation, see the report. In particular, the
+code implements a pure elliptic solver. Both solvers are implicit.
+
+Files Included
+============================================================================
+- kinematic_viscosity.py: The main python module.
+- kinematic_viscosity_ext.c: The C extension to the python module.
+- test_kinematic_viscosity.py: A unit test for the python module.
+- util_ext.*: Copies of the files in anuga.utilities (for kinematic_viscosity_ext.c).
+- sparse.py: A modified version of anuga.utilities.sparse (which will supersede it).
+- test_sparse.py: A unit test for the sparse module.
+- make.bat: A makefile script for MinGW on Windows (this may have to be modified).
+
+Usage
+============================================================================
+Some sample code can be found in the unit tests. In particular, one 
+requires an anuga.abstract_2d_finite_volumes.domain domain, a numpy 
+vector h with h.shape == (n,) of stage heights and momentum numpy matrix 
+p with p == [(hu) | (hv)] where u and v are x- and y-velocity and
+p.shape == (n,2).
+
+The stage height vector is the stage heights for the current time step.
+
+To initialise the kinematic viscosity operator, one uses
+kv_operator = kinematic_viscosity.Kinematic_Viscosity_Operator(domain)
+with optional parameters
+
+1. apply_triangle_areas (default True)
+Whether or not to include the factor of 1/|Triangle area| in the numerical 
+approximation (see the report for more details).
+
+2. verbose (default False)
+Whether to include detailed logs for anuga.log.
+
+Parabolic Solver:
+We use the commands
+> kv_operator.apply_stage_heights(h)
+> p_new = kv_operator.parabolic_solver(p)
+
+Elliptic Solver:
+This solves
+  (h*u_x)_x + (h*u_y)_y = f_1
+  (h*v_y)_y + (h*v_y)_y = f_2
+We also need a right-hand side numpy matrix f with f.shape == (n,2)
+and f == [f_1 | f_2]. Then
+> kv_operator.apply_stage_heights(h)
+> p = kv_operator.cg_solve(f)
+You can also solve for in only one dimension with
+> kv_operator.apply_stage_heights(h)
+> kv_operator.set_qty_considered('u')
+> p_x = kv_operator.cg_solve(f1)
+(and similarly for v)
+
+============================================================================
+
+Last Modified: 25 October 2010
Index: trunk/anuga_work/development/2010-projects/kv/sparse.py
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/sparse.py	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/sparse.py	(revision 8051)
@@ -0,0 +1,401 @@
+"""Proof of concept sparse matrix code
+"""
+#Exactly as in the ANUGA code *except* for the initialiser for Sparse_CSR
+#Note: Sparse_CSR works *only* if there is a nonzero entry in every row!
+
+import numpy as num
+
+
+class Sparse:
+
+    def __init__(self, *args):
+        """Create sparse matrix.
+        There are two construction forms
+        Usage:
+
+        Sparse(A)     #Creates sparse matrix from dense matrix A 
+        Sparse(M, N)  #Creates empty MxN sparse matrix
+        """
+
+        self.Data = {}
+            
+        if len(args) == 1:
+            try:
+                A = num.array(args[0])
+            except:
+                raise 'Input must be convertable to a numeric array'
+
+            assert len(A.shape) == 2, 'Input must be a 2d matrix'
+            
+            self.M, self.N = A.shape
+            for i in range(self.M):
+                for j in range(self.N):
+                    if A[i, j] != 0.0:
+                        self.Data[i, j] = A[i, j]
+                
+            
+        elif len(args) == 2:
+            self.M = args[0]
+            self.N = args[1]
+        else:
+            raise 'Invalid construction'
+            
+        self.shape = (self.M, self.N) 
+
+
+    def __repr__(self):
+        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.Data`
+
+    def __len__(self):
+        """Return number of nonzeros of A
+        """
+        return len(self.Data)
+
+    def nonzeros(self):
+        """Return number of nonzeros of A
+        """        
+        return len(self)
+    
+    def __setitem__(self, key, x):
+
+        i,j = key
+        # removing these asserts will not speed things up
+        assert 0 <= i < self.M
+        assert 0 <= j < self.N        
+
+        if x != 0:
+            self.Data[key] = float(x)
+        else:
+            if self.Data.has_key( key ):            
+                del self.Data[key]
+
+    def __getitem__(self, key):
+        
+        i,j = key
+        # removing these asserts will not speed things up
+        assert 0 <= i < self.M
+        assert 0 <= j < self.N                
+
+        if self.Data.has_key( key ):
+            return self.Data[ key ]
+        else:
+            return 0.0
+
+    def copy(self):
+        #FIXME: Use the copy module instead
+        new = Sparse(self.M,self.N)
+
+        for key in self.Data.keys():
+            i, j = key
+
+            new[i,j] = self.Data[i,j]
+
+        return new
+
+
+    def todense(self):
+        D = num.zeros( (self.M, self.N), num.float)
+        
+        for i in range(self.M):
+            for j in range(self.N):
+                if self.Data.has_key( (i,j) ):                
+                    D[i, j] = self.Data[ (i,j) ]
+        return D
+
+
+    
+    def __mul__(self, other):
+        """Multiply this matrix onto 'other' which can either be
+        a numeric vector, a numeric matrix or another sparse matrix.
+        """
+
+        try:
+            B = num.array(other)
+        except:
+            msg = 'FIXME: Only numeric types implemented so far'
+            raise msg
+            
+
+        # Assume numeric types from now on
+	
+        if len(B.shape) == 0:
+            # Scalar - use __rmul__ method
+            R = B*self
+	    
+        elif len(B.shape) == 1:
+            # Vector
+            msg = 'Mismatching dimensions: You cannot multiply (%d x %d) matrix onto %d-vector'\
+                  %(self.M, self.N, B.shape[0])
+            assert B.shape[0] == self.N, msg
+
+            R = num.zeros(self.M, num.float) #Result
+	    
+            # Multiply nonzero elements
+            for key in self.Data.keys():
+                i, j = key
+
+                R[i] += self.Data[key]*B[j]
+        elif len(B.shape) == 2:
+	
+            
+            R = num.zeros((self.M, B.shape[1]), num.float) #Result matrix
+
+            # Multiply nonzero elements
+	    for col in range(R.shape[1]):
+	        # For each column
+		
+                for key in self.Data.keys():
+                    i, j = key
+
+                    R[i, col] += self.Data[key]*B[j, col]
+	    
+	    
+        else:
+            raise ValueError, 'Dimension too high: d=%d' %len(B.shape)
+
+        return R
+    
+
+    def __add__(self, other):
+        """Add this matrix onto 'other' 
+        """
+
+        new = other.copy()
+        for key in self.Data.keys():
+            i, j = key
+
+            new[i,j] += self.Data[key]
+
+        return new
+
+
+    def __rmul__(self, other):
+        """Right multiply this matrix with scalar
+        """
+
+        try:
+            other = float(other)
+        except:
+            msg = 'Sparse matrix can only "right-multiply" onto a scalar'
+            raise TypeError, msg
+        else:
+            new = self.copy()
+            #Multiply nonzero elements
+            for key in new.Data.keys():
+                i, j = key
+
+                new.Data[key] = other*new.Data[key]
+
+        return new
+
+
+    def trans_mult(self, other):
+        """Multiply the transpose of matrix with 'other' which can be
+        a numeric vector.
+        """
+
+        try:
+            B = num.array(other)
+        except:
+            print 'FIXME: Only numeric types implemented so far'
+
+
+        #Assume numeric types from now on
+        if len(B.shape) == 1:
+            #Vector
+
+            assert B.shape[0] == self.M, 'Mismatching dimensions'
+
+            R = num.zeros((self.N,), num.float) #Result
+
+            #Multiply nonzero elements
+            for key in self.Data.keys():
+                i, j = key
+
+                R[j] += self.Data[key]*B[i]
+
+        else:
+            raise 'Can only multiply with 1d array'
+
+        return R
+
+class Sparse_CSR:
+
+    def __init__(self, A=None, data=None, Colind=None, rowptr=None, m=None, n=None):
+        """Create sparse matrix in csr format.
+
+        Sparse_CSR(A) #creates csr sparse matrix from sparse matrix
+        Matrices are not built using this format, since it's painful to
+        add values to an existing sparse_CSR instance (hence there are no
+        objects to do this.)
+
+        Rather, build a matrix, and convert it to this format for a speed
+        increase.
+
+        data - a 1D array of the data
+        Colind - The ith item in this 1D array is the column index of the
+                 ith data in the data array
+        rowptr - 1D array, with the index representing the row of the matrix.
+                 The item in the row represents the index into colind of the
+                 first data value of this row.
+                 Regard it as a pointer into the colind array, for the ith row.
+
+                 
+        """
+
+        if isinstance(A,Sparse):
+
+            keys = A.Data.keys()
+            keys.sort()
+            nnz = len(keys)
+            data    = num.zeros ( (nnz,), num.float)
+            colind  = num.zeros ( (nnz,), num.int)
+            row_ptr = num.zeros ( (A.M+1,), num.int)
+            current_row = -1
+            k = 0
+            for key in keys:
+                ikey0 = int(key[0])
+                ikey1 = int(key[1])
+                if ikey0 != current_row:
+                    current_row = ikey0
+                    row_ptr[ikey0] = k
+                data[k] = A.Data[key]
+                colind[k] = ikey1
+                k += 1
+            for row in range(current_row+1, A.M+1):
+                row_ptr[row] = nnz
+            #row_ptr[-1] = nnz
+        
+            self.data    = data
+            self.colind  = colind
+            self.row_ptr = row_ptr
+            self.M       = A.M
+            self.N       = A.N
+        elif isinstance(data,num.ndarray) and isinstance(Colind,num.ndarray) and isinstance(rowptr,num.ndarray) and isinstance(m,int) and isinstance(n,int):
+            msg = "Sparse_CSR: data is array of wrong dimensions"
+            #assert len(data.shape) == 1, msg
+            nnz = data.size
+
+            msg = "Sparse_CSR: Colind is array of wrong dimensions"
+            assert Colind.shape == (nnz,), msg
+
+            msg = "Sparse_CSR: rowptr is array of wrong dimensions"
+            assert rowptr.shape == (m+1,), msg
+
+            self.data = data
+            self.colind = Colind
+            self.row_ptr = rowptr
+            self.M = m
+            self.N = n
+        else:
+            raise ValueError, "Sparse_CSR(A) expects A == Sparse Matrix *or* data==array,colind==array,rowptr==array,m==int,n==int"
+            
+    def __repr__(self):
+        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.data`
+
+    def __len__(self):
+        """Return number of nonzeros of A
+        """
+        return self.row_ptr[-1]
+
+    def nonzeros(self):
+        """Return number of nonzeros of A
+        """        
+        return len(self)
+
+    def todense(self):
+        D = num.zeros( (self.M, self.N), num.float)
+        
+        for i in range(self.M):
+            for ckey in range(self.row_ptr[i],self.row_ptr[i+1]):
+                j = self.colind[ckey]
+                D[i, j] = self.data[ckey]
+        return D
+
+    def __mul__(self, other):
+        """Multiply this matrix onto 'other' which can either be
+        a numeric vector, a numeric matrix or another sparse matrix.
+        """
+
+        try:
+            B = num.array(other)
+        except:
+            print 'FIXME: Only numeric types implemented so far'
+
+        return csr_mv(self,B) 
+
+
+# Setup for C extensions
+from anuga.utilities import compile
+if compile.can_use_C_extension('sparse_ext.c'):
+    # Access underlying c implementations
+    from sparse_ext import csr_mv
+
+
+if __name__ == '__main__':
+    # A little selftest
+    
+    A = Sparse(3,3)
+
+    A[1,1] = 4
+
+
+    print A
+    print A.todense()
+
+    A[1,1] = 0
+
+    print A
+    print A.todense()    
+
+    A[1,2] = 0
+
+
+    A[0,0] = 3
+    A[1,1] = 2
+    A[1,2] = 2
+    A[2,2] = 1
+
+    print A
+    print A.todense()
+
+
+    #Right hand side vector
+    v = [2,3,4]
+
+    u = A*v
+    print u
+    assert num.allclose(u, [6,14,4])
+
+    u = A.trans_mult(v)
+    print u
+    assert num.allclose(u, [6,6,10])
+
+    #Right hand side column
+    v = num.array([[2,4],[3,4],[4,4]])
+
+    u = A*v[:,0]
+    assert num.allclose(u, [6,14,4])
+
+    #u = A*v[:,1]
+    #print u
+    print A.shape
+
+    B = 3*A
+    print B.todense()
+
+    B[1,0] = 2
+
+    C = A+B
+
+    print C.todense()
+
+    C = Sparse_CSR(C)
+
+    y = C*[6,14,4]
+
+    print y
+
+    y2 = C*[[6,4],[4,28],[4,8]]
+
+    print y2
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 8051)
@@ -0,0 +1,153 @@
+from anuga.abstract_2d_finite_volumes.domain import Domain
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+from kinematic_viscosity import Kinematic_Viscosity_Operator
+import numpy as num
+from math import sqrt
+import unittest
+
+class Test_Kinematic_Viscosity(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])
+		domain = 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])
+		domain = Domain(source=points,triangles=elements,boundary=boundary_map)
+		return Kinematic_Viscosity_Operator(domain)
+
+	def test_enumerate_boundary(self):
+		operator1 = self.operator1()
+		boundary_enum = operator1.boundary_enum
+		assert boundary_enum[(0,0)] == 0
+		assert boundary_enum[(0,1)] == 1
+		assert boundary_enum[(0,2)] == 2
+		
+		operator2 = self.operator2()
+		boundary_enum = operator2.boundary_enum
+		assert boundary_enum[(0,1)] == 0
+		assert boundary_enum[(0,2)] == 1
+		assert boundary_enum[(1,0)] == 2
+		assert boundary_enum[(1,2)] == 3
+
+	def test_geo_structure(self):
+		operator1 = self.operator1()
+		indices = operator1.geo_structure_indices
+		values = operator1.geo_structure_values
+		assert num.allclose(indices, num.array([[1, 2, 3]]))
+		assert num.allclose(values, num.array([[-6.0, -6.0/sqrt(5), -6.0/sqrt(5)]]))
+		
+		operator2 = self.operator2()
+		indices = operator2.geo_structure_indices
+		values = operator2.geo_structure_values
+		assert num.allclose(indices, num.array([[1,2,3],[4,0,5]]))
+		assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]]))
+		
+	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,  6.0/sqrt(5), 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, 6.0/sqrt(5), 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,6.0/sqrt(5.0),0.0,0.0],\
+						[0.0, -6.0/sqrt(5.0), 0.0, 0.0, 0.0, 6.0/sqrt(5.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_elliptic_multiply(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 <- Centriod value
+		V2 = num.array([2.0]) #(vh)=2 <- Centroid value
+		A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+		U1 = num.array([[2.0],[2.0],[1.0],[1.0]]) #(uh) for centroid, 3 edges
+		U2 = num.array([[2.0],[1.0],[2.0],[0.0]]) #(vh) for centroid, 3 edge
+		
+		operator1.set_qty_considered('u')
+		D1 = operator1.elliptic_multiply(V1)
+		assert num.allclose(D1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) #2* for triangle_areas
+		
+		operator1.set_qty_considered(2)
+		D2 = operator1.elliptic_multiply(V2)
+		assert num.allclose(D2, 2*num.array(num.mat(A)*num.mat(U2)).reshape(1,)) #2* for triangle_areas
+	
+	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,  6.0/sqrt(5), 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,  6.0/sqrt(5), 6.0/sqrt(5)])
+		U = num.array([[2.0,2.0],[2.0,1.0],[1.0,2.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)
+	
+	def test_parabolic_solve(self):
+		operator1 = self.operator1()
+		operator1.apply_stage_heights(num.array([[1.0]])) #h=1
+		A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+		U = num.array([[2.0,1.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
+		u = num.array([[2.0,1.0]])
+		U_new = operator1.parabolic_solver(u)
+		U_mod = num.array([[0.0,0.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
+		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')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: trunk/anuga_work/development/2010-projects/kv/test_sparse.py
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/test_sparse.py	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/test_sparse.py	(revision 8051)
@@ -0,0 +1,214 @@
+#!/usr/bin/env python
+#Same as ANUGA code, except testing the new Sparse_CSR initialisation routine
+import unittest
+from math import sqrt
+
+from sparse import *
+import numpy as num
+
+
+class Test_Sparse(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+    def test_init1(Self):
+        """Test initialisation from dimensions
+        """
+        A = Sparse(3,3)
+        A[1,1] = 4
+
+        for i in range(3):
+            for j in range(3):
+                if i==1 and j==1:
+                    assert A[i,j] == 4.0
+                else:
+                    assert A[i,j] == 0.0
+
+    def test_init2(Self):
+        """Test initialisation from dense matrix
+        """
+
+        A = Sparse(4,3)
+        A[1,1] = 4
+        A[2,0] = 6
+        A[0,1] = 13
+        A[0,2] = -6
+        A[2,2] = 1
+
+        B = A.todense()
+
+        C = Sparse(B)
+
+        assert num.allclose(C.todense(), B)
+
+
+    def test_dense(self):
+        A = Sparse(4,3)
+        A[1,1] = 4
+
+        assert num.allclose(A.todense(), [[0,0,0], [0,4,0], [0,0,0], [0,0,0]])
+
+
+    def test_reset_to_zero_possible(self):
+        """Test that nonzero element can be reset to zero
+        (This was not possible when we tried sparse from scipy).
+        """
+
+        A = Sparse(3,3)
+        A[1,1] = 4
+        A[1,1] = 0
+
+        assert len(A) == 0
+        assert num.allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
+
+        #Set an existing zero element to zero
+        A[1,2] = 0
+        assert len(A) == 0
+        assert num.allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]])
+
+    def test_sparse_multiplication_vector(self):
+        A = Sparse(3,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+        #Right hand side vector
+        v = [2,3,4]
+
+        u = A*v
+        assert num.allclose(u, [6,14,4])
+
+        #Right hand side column
+        v = num.array([[2,4],[3,4],[4,4]])
+
+        u = A*v[:,0]
+        assert num.allclose(u, [6,14,4])
+
+        u = A*v[:,1]
+        assert num.allclose(u, [12,16,4])
+
+
+    def test_sparse_multiplication_matrix(self):
+        A = Sparse(3,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+        #Right hand side matrix
+        v = num.array([[2,4],[3,4],[4,4]])
+
+        u = A*v
+        assert num.allclose(u, [[6,12], [14,16], [4,4]])
+
+
+
+    def test_sparse_transpose_multiplication(self):
+        A = Sparse(3,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+        #Right hand side vector
+        v = [2,3,4]
+
+        u = A.trans_mult(v)
+        assert num.allclose(u, [6,6,10])
+
+
+    def test_scalar_multiplication(self):
+        """Test method __rmul__
+        """
+
+        A = Sparse(3,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+        B = 3*A
+        assert num.allclose(B.todense(), 3*A.todense())
+
+        B = A*3
+        assert num.allclose(B.todense(), 3*A.todense())
+
+        try:
+            B = 'a'*A
+        except TypeError:
+            pass
+        else:
+            raise 'Should have failed'
+
+    def test_sparse_addition(self):
+        """ Test sparse addition with dok format
+        """
+
+        A = Sparse(3,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+
+
+        B = 3*A
+        B[1,0] = 2
+
+        C = A+B
+
+        assert num.allclose(C.todense(), [[12,0,0], [2,8,8], [0,0,4]])
+
+    def test_sparse_tocsr(self):
+        """ Test conversion to csr format
+        """
+
+        A = Sparse(4,3)
+
+        A[0,0] = 3
+        A[1,1] = 2
+        A[1,2] = 2
+        A[2,2] = 1
+        A[0,2] = 4
+        A[2,0] = 5
+
+        #print ' '
+        #print A.todense()
+
+        B = Sparse_CSR(A)
+
+        #print B.todense()
+
+        C = [1, 2, 3]
+
+        assert num.allclose(B*C, [15.0, 10.0 ,8.0, 0.0])
+
+        C2 = [[1,2],[2,4],[3,6]]
+
+        #print B*C2
+
+        assert num.allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]])
+
+    def test_sparse_csr_init(self):
+        A = num.array([[1.0,0.0,-1.0,0.0],[0.0,2.0,0.0,0.0],[0.0,0.0,0.0,-3.0],[0.0,0.0,4.0,0.0]])
+        data = num.array([1.0,-1.0,2.0,-3.0,4.0])
+        Colind = num.array([0,2,1,3,2])
+        rowptr = num.array([0,2,3,4,5]) #the 5 does not correspond to any row, it is just there so we know how many nonzero entries there are!
+        A_CSR = Sparse_CSR(None,data,Colind,rowptr,4,4)
+        A_dense = A_CSR.todense()
+        assert num.allclose(A, A_dense)
+################################################################################
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Sparse, 'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: trunk/anuga_work/development/2010-projects/kv/util_ext.c
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/util_ext.c	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/util_ext.c	(revision 8051)
@@ -0,0 +1,108 @@
+// Python - C extension for finite_volumes util module.
+//
+// To compile (Python2.3):
+//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
+//  gcc -shared util_ext.o  -o util_ext.so
+//
+// See the module util.py
+//
+//
+// Ole Nielsen, GA 2004
+//
+//NOTE: On 64 bit systems use long* instead of int* for numeric arrays
+//this will also work on 32 bit systems
+
+#include <float.h>
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+
+//Shared snippets
+#include "util_ext.h"
+
+
+
+
+PyObject *gradient(PyObject *self, PyObject *args) {
+  //
+  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+  //
+
+  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
+  PyObject *result;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
+			&q0, &q1, &q2)) {
+    
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "gradient could not parse input");    
+    return NULL;
+  }
+
+
+  // Call underlying routine
+  _gradient(x0, y0, x1, y1, x2, y2,
+	    q0, q1, q2, &a, &b);
+
+  // Return values a and b
+  result = Py_BuildValue("dd", a, b);
+  return result;
+}
+
+PyObject *gradient2(PyObject *self, PyObject *args) {
+  //
+  // a,b = gradient2(x0, y0, x1, y1, q0, q1)
+  //
+
+  double x0, y0, x1, y1, q0, q1, a, b;
+  PyObject *result;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddddd", &x0, &y0, &x1, &y1, &q0, &q1)) {
+    PyErr_SetString(PyExc_RuntimeError, 
+		    "gradient2 could not parse input");      
+    return NULL;
+  }
+
+
+  // Call underlying routine
+  _gradient2(x0, y0, x1, y1, q0, q1, &a, &b);
+
+  // Return values a and b
+  result = Py_BuildValue("dd", a, b);
+  return result;
+}
+
+PyObject * double_precision(PyObject * self, PyObject * args){
+  // Get the precision of the double datatype on this system.
+  return Py_BuildValue("i", DBL_DIG);
+}
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+
+  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"gradient", gradient, METH_VARARGS, "Print out"},
+  {"gradient2", gradient2, METH_VARARGS, "Print out"},  
+  {"double_precision", double_precision, METH_VARARGS, "Precision of this machine\'s \'double\' type"},
+  {NULL, NULL, 0, NULL}   /* sentinel */
+};
+
+
+
+// Module initialisation
+void initutil_ext(void){
+  Py_InitModule("util_ext", MethodTable);
+
+  import_array();     //Necessary for handling of NumPY structures
+}
+
+
+
+
Index: trunk/anuga_work/development/2010-projects/kv/util_ext.h
===================================================================
--- trunk/anuga_work/development/2010-projects/kv/util_ext.h	(revision 8051)
+++ trunk/anuga_work/development/2010-projects/kv/util_ext.h	(revision 8051)
@@ -0,0 +1,367 @@
+// Python - C extension for finite_volumes util module.
+//
+// To compile (Python2.3):
+//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
+//  gcc -shared util_ext.o  -o util_ext.so	
+//
+// See the module util.py
+//
+//
+// Ole Nielsen, GA 2004
+	
+#include "Python.h"	
+#include "numpy/arrayobject.h"
+#include "math.h"
+
+
+double max(double x, double y) {  
+  //Return maximum of two doubles
+  
+  if (x > y) return x;
+  else return y;
+}
+
+
+double min(double x, double y) {  
+  //Return minimum of two doubles
+  
+  if (x < y) return x;
+  else return y;
+}
+
+
+int _gradient(double x0, double y0, 
+	      double x1, double y1, 
+	      double x2, double y2, 
+	      double q0, double q1, double q2, 
+	      double *a, double *b) {
+	      
+  /*Compute gradient (a,b) based on three points (x0,y0), (x1,y1) and (x2,y2) 
+  with values q0, q1 and q2.
+  
+  Extrapolation formula (q0 is selected as an arbitrary origin)
+    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
+  
+  Substituting the known values for q1 and q2 into (1) yield the 
+  equations for a and b 
+  
+      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
+      q2-q0 = a*(x2-x0) + b*(y2-y0)                      (3)      
+      
+  or in matrix form
+  
+  /               \  /   \   /       \  
+  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
+  |               |  |   | = |       | 
+  |  x2-x0  y2-y0 |  | b |   | q2-q0 |
+  \               /  \   /   \       /
+   
+  which is solved using the standard determinant technique    
+      
+  */
+	      
+
+  double det;
+  
+  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
+
+  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
+  *a /= det;
+
+  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
+  *b /= det;
+
+  return 0;
+}
+
+
+int _gradient2(double x0, double y0, 
+	       double x1, double y1, 
+	       double q0, double q1, 
+	       double *a, double *b) {
+  /*Compute gradient (a,b) between two points (x0,y0) and (x1,y1) 
+  with values q0 and q1 such that the plane is constant in the direction 
+  orthogonal to (x1-x0, y1-y0).
+  
+  Extrapolation formula
+    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
+  
+  Substituting the known values for q1 into (1) yields an 
+  under determined  equation for a and b 
+      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
+      
+      
+  Now add the additional requirement that the gradient in the direction 
+  orthogonal to (x1-x0, y1-y0) should be zero. The orthogonal direction 
+  is given by the vector (y0-y1, x1-x0).
+  
+  Define the point (x2, y2) = (x0 + y0-y1, y0 + x1-x0) on the orthognal line. 
+  Then we know that the corresponding value q2 should be equal to q0 in order 
+  to obtain the zero gradient, hence applying (1) again    
+      q0 = q2 = q(x2, y2) = q0 + a*(x2-x0) + b*(y2-y0)
+                          = q0 + a*(x0 + y0-y1-x0) + b*(y0 + x1-x0 - y0)
+			  = q0 + a*(y0-y1) + b*(x1-x0)
+			  
+  leads to the orthogonality constraint
+     a*(y0-y1) + b*(x1-x0) = 0                           (3) 
+     
+  which closes the system and yields
+  
+  /               \  /   \   /       \  
+  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
+  |               |  |   | = |       | 
+  |  y0-y1  x1-x0 |  | b |   |   0   |
+  \               /  \   /   \       /
+   
+  which is solved using the standard determinant technique    
+      
+  */
+
+  double det, xx, yy, qq;
+  
+  xx = x1-x0;
+  yy = y1-y0;
+  qq = q1-q0;
+    
+  det = xx*xx + yy*yy;  //FIXME  catch det == 0
+  *a = xx*qq/det;
+  *b = yy*qq/det;
+        
+  return 0;
+}
+
+
+void _limit_old(int N, double beta, double* qc, double* qv, 
+	    double* qmin, double* qmax) { 
+
+  //N are the number of elements
+  int k, i, k3;
+  double dq, dqa[3], phi, r;
+  
+  //printf("INSIDE\n");
+  for (k=0; k<N; k++) {
+    k3 = k*3;
+    
+    //Find the gradient limiter (phi) across vertices  
+    phi = 1.0;
+    for (i=0; i<3; i++) {    
+      r = 1.0;
+      
+      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
+      dqa[i] = dq;              //Save dq for use in the next loop
+      
+      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
+      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
+  
+  
+      phi = min( min(r*beta, 1.0), phi);    
+    }
+    
+    //Then update using phi limiter
+    for (i=0; i<3; i++) {    
+      qv[k3+i] = qc[k] + phi*dqa[i];
+    }
+  }
+}
+
+
+void  print_double_array(char* name, double* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	    printf("%g ",array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+void  print_int_array(char* name, int* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	    printf("%i ",array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+
+void  print_long_array(char* name, long* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+	km = m*k;
+	printf("[");
+	for (i=0; i<m ; i++){
+	  printf("%i ",(int) array[km+i]);
+	}
+	if (k==(n-1))
+	    printf("]");
+	else
+	    printf("]\n");
+    }
+    printf("]\n");
+}
+
+void print_numeric_array(PyArrayObject *x) {  
+  int i, j;
+  for (i=0; i<x->dimensions[0]; i++) {  
+    for (j=0; j<x->dimensions[1]; j++) {
+      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
+    }
+    printf("\n");  
+  }
+  printf("\n");    
+}
+
+void print_numeric_vector(PyArrayObject *x) {  
+  int i;
+  for (i=0; i<x->dimensions[0]; i++) {
+    printf("%f ", *(double*) (x->data + i*x->strides[0]));  
+  }
+  printf("\n");  
+}
+
+PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
+  PyArrayObject *A, *B;
+  
+
+  //Get array object from attribute
+  
+  /*
+  //FIXME: THE TEST DOESN't WORK
+  printf("Err = %d\n", PyObject_HasAttrString(O, name));
+  if (PyObject_HasAttrString(O, name) == 1) {
+    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
+    if (!B) return NULL;
+  } else {
+    return NULL;
+    } 
+  */
+    
+  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
+
+  //printf("B = %p\n",(void*)B);
+  if (!B) {
+    printf("util_ext.h: get_consecutive_array could not obtain python object");
+    printf(" %s\n",name);
+    fflush(stdout);
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain python object");
+    return NULL;
+  }     
+  
+  //Convert to consecutive array
+  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
+						    B -> descr -> type, 0, 0);
+  
+  Py_DECREF(B); //FIXME: Is this really needed??
+  
+  if (!A) {
+    printf("util_ext.h: get_consecutive_array could not obtain array object");
+    printf(" %s \n",name);
+    fflush(stdout);
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }
+
+
+  return A;
+}
+
+double get_python_double(PyObject *O, char *name) {
+  PyObject *TObject;
+  #define BUFFER_SIZE 80
+  char buf[BUFFER_SIZE];
+  double tmp;
+  int n;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+	n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
+	//printf("name = %s",name);
+    PyErr_SetString(PyExc_RuntimeError, buf);
+
+    return 0.0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
+
+
+
+int get_python_integer(PyObject *O, char *name) {
+  PyObject *TObject;
+  #define BUFFER_SIZE 80
+  char buf[BUFFER_SIZE];
+  long tmp;
+  int n;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+  	n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
+	//printf("name = %s",name);
+    PyErr_SetString(PyExc_RuntimeError, buf);
+    return 0;
+  }  
+  
+  tmp = PyInt_AsLong(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
+
+PyObject *get_python_object(PyObject *O, char *name) {
+  PyObject *Oout;
+
+  Oout = PyObject_GetAttrString(O, name);
+  if (!Oout) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
+    return NULL;
+  }
+
+  return Oout;
+
+}
+
+
+// check that numpy array objects are C contiguous memory
+#define CHECK_C_CONTIG(varname)	if (!PyArray_ISCONTIGUOUS(varname)) { \
+				    char msg[1024]; \
+				    sprintf(msg, \
+					    "%s(): file %s, line %d: " \
+				            "'%s' object is not C contiguous memory", \
+				             __func__, __FILE__, __LINE__, #varname); \
+				    PyErr_SetString(PyExc_RuntimeError, msg); \
+				    return NULL; \
+				}
