Index: /inundation/README.txt
===================================================================
--- /inundation/README.txt	(revision 1741)
+++ /inundation/README.txt	(revision 1742)
@@ -40,4 +40,9 @@
   
 
+PATH
+
+It is assumed that the root directory for all anuga/inundation modules is on the pythonpath  
+  
+  
 LINKS  
     Python: http://www.python.org
Index: /inundation/validation/LWRU2/README.txt
===================================================================
--- /inundation/validation/LWRU2/README.txt	(revision 1742)
+++ /inundation/validation/LWRU2/README.txt	(revision 1742)
@@ -0,0 +1,15 @@
+VALIDATION CODE
+
+This directory contains code conduction a validation of pyvolution aganst a
+wave tank simulation of the 1993 Okushiri island tsunami as described in the file problem02.doc and URL 
+
+http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
+
+Run the code lwru2.py with the variable read_mesh = False (at the top
+of the code) and then extract_timeseries.py to see the results.
+
+If you wish to look at the a more detailed animation, create a mesh
+using the script create_mesh.py and run lwru2.py again with read_mesh
+= True. The resulting sww file can be viewed with e.g. Swollen.
+
+
Index: /inundation/validation/LWRU2/create_mesh.py
===================================================================
--- /inundation/validation/LWRU2/create_mesh.py	(revision 1741)
+++ /inundation/validation/LWRU2/create_mesh.py	(revision 1742)
@@ -14,5 +14,5 @@
 
     #Basic geometry
-
+    
     xleft   = 0
     xright  = 5.448
@@ -20,27 +20,23 @@
     ytop    = 3.402
 
-    #Outline
+    #Outline 
     point_sw = [xleft, ybottom]
     point_se = [xright, ybottom]
-    point_nw = [xleft, ytop]
+    point_nw = [xleft, ytop]    
     point_ne = [xright, ytop]
 
     #Midway points (left)
-    point_ml1 = [xleft + (xright-xleft)/3+1, ytop]
-    point_ml2 = [xleft + (xright-xleft)/3+1, ybottom]
-
-    #Midway points (right)
-    point_mr1 = [xleft + 2*(xright-xleft)/3+0.5, ytop]
-    point_mr2 = [xleft + 2*(xright-xleft)/3+0.5, ybottom]
+    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
+    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
 
 
     geo = Geo_reference(xllcorner = xleft,
                         yllcorner = ybottom)
-
-
+    
+                        
     print "***********************"
     print "geo ref", geo
     print "***********************"
-
+    
     m = Mesh(geo_reference=geo)
 
@@ -49,18 +45,15 @@
     dict['points'] = [point_se,   #se
                       point_ne,
-                      point_mr1,
-                      point_ml1,
+                      point_mtop,
                       point_nw,
                       point_sw,
-                      point_ml2,
-                      point_mr2]
+                      point_mbottom]
 
-
+    
     dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
-                        [4,5], [5,6], [6,7], [7,0], #The outer border
-                        [2,7], [3,6]]               #Separators
-
+                        [4,5], [5,0],  #The outer border
+                        [2,5]]         #Separator
+    
     dict['segment_tags'] = ['wall',
-                            'wall',
                             'wall',
                             'wall',
@@ -68,10 +61,10 @@
                             'wall',
                             'wall',
-                            'wall',
-                            '',           #Interior
                             '']           #Interior
 
+        
+    m.addVertsSegs(dict)
 
-    m.addVertsSegs(dict)
+
 
 
@@ -87,32 +80,58 @@
     p2 = [xr, yt]
     p3 = [xl, yt]
+    
+    dict['points'] = [p0, p1, p2, p3]
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
+    dict['segment_tags'] = ['', '', '', '']
+    m.addVertsSegs(dict)    
 
-    dict['points'] = [p0, p1, p2, p3]
+    #Island area
+    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
+    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
+    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
+    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
+    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3] 
+    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
+    island_6 = [xl-.01, yb]  #OK
+    island_7 = [xl-.01, yt]  #OK      
+    
+
+    dict['points'] = [island_0, island_1, island_2,
+                      island_3, island_4, island_5,
+                      #p0, p3]                      
+                      island_6, island_7]
+
+                      
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
+                        [4,5], [5,6], [6,7], [7,0]]
+                        
+    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
 
 
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
-    dict['segment_tags'] = ['', '', '', '']
-    m.addVertsSegs(dict)
+    m.addVertsSegs(dict)    
 
+    
+#
+    
+    base_resolution = 1
 
-    base_resolution = 100
+    ocean = m.addRegionEN(xleft+.1, ybottom+.1)
+    ocean.setMaxArea(0.1*base_resolution)
 
-    ocean = m.addRegionEN(xleft+1, ybottom+1)
-    ocean.setMaxArea(0.01*base_resolution)
+    mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
+    mid.setMaxArea(0.0005*base_resolution)
+    
 
-    mid = m.addRegionEN(point_ml2[0]+1, ybottom+1)
-    mid.setMaxArea(0.001*base_resolution)
-
-
-    inner = m.addRegionEN(point_mr2[0]+1, ybottom+1)
-    inner.setMaxArea(0.0001*base_resolution)
+    inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
+    inner.setMaxArea(0.00007*base_resolution)
 
 
     gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
-    gulleys.setMaxArea(0.00001*base_resolution)
-
-
+    gulleys.setMaxArea(0.00002*base_resolution)    
+    
+    
     m.generateMesh('pzq28.0za1000000a')
 
     import filenames
     m.export_mesh_file(filenames.mesh_filename)
+    
Index: /inundation/validation/LWRU2/extract_timeseries.py
===================================================================
--- /inundation/validation/LWRU2/extract_timeseries.py	(revision 1741)
+++ /inundation/validation/LWRU2/extract_timeseries.py	(revision 1742)
@@ -9,8 +9,6 @@
 gauges = [[0.000, 1.696]]  #Boundary
 gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
-#gauges = [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
 
 gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
-#gauge_names = ['ch5', 'ch7', 'ch9']
 
 
Index: /inundation/validation/LWRU2/lwru2.py
===================================================================
--- /inundation/validation/LWRU2/lwru2.py	(revision 1741)
+++ /inundation/validation/LWRU2/lwru2.py	(revision 1742)
@@ -77,76 +77,72 @@
 from Numeric import array, zeros, Float, allclose
 import filenames
+from caching import cache
 
-picklefile = 'domain.pck'
-from cPickle import dump, load
+#Preparing time boundary
+prepare_timeboundary(filenames.boundary_filename)
 
-N = 100
-
-try:
-    #raise Exception
-    fid = open(picklefile)
-    print 'Read pickled domain'
-    domain = load(fid)
-
-except:
-
-    #Preparing time boundary
-    prepare_timeboundary(filenames.boundary_filename)
-
-    #Preparing points
-    from pyvolution.data_manager import xya2pts
-    xya2pts(filenames.bathymetry_filename, verbose = True,
-            z_func = lambda z: -z)
+#Preparing points
+from pyvolution.data_manager import xya2pts
+xya2pts(filenames.bathymetry_filename, verbose = True,
+        z_func = lambda z: -z)
 
 
-    #######################
-    # Domain
+#######################
+# Domain
     
-    if read_mesh is True:
-    	print 'Creating domain from', filenames.mesh_filename
-    	domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)    
-    else:	
-    	print 'Creating regular domain'
-        points, vertices, boundary = rectangular(N, N/5*3,
-    	                                         len1=5.448, len2=3.402)
-    	domain = Domain(points, vertices, boundary)
+if read_mesh is True:
+    print 'Creating domain from', filenames.mesh_filename
+    domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)    
+else:	
+    print 'Creating regular mesh'
+    N = 100
+    points, vertices, boundary = rectangular(N, N/5*3,
+                                             len1=5.448, len2=3.402)
+    print 'Creating domain'
+    
+    #domain = Domain(points, vertices, boundary)
+    domain = cache(Domain, (points, vertices, boundary))
 
 
 
-    print "Number of triangles = ", len(domain)
-    print 'The extent is ', domain.get_extent()
+print "Number of triangles = ", len(domain)
+print 'The extent is ', domain.get_extent()
 
 
-    domain.check_integrity()
-    domain.default_order = 2
+#domain.check_integrity()
+domain.default_order = 2
 
-    print "Number of triangles = ", len(domain)
-    domain.store = True    #Store for visualisation purposes
+print "Number of triangles = ", len(domain)
+domain.store = True    #Store for visualisation purposes
 
-    import sys, os
-    base = os.path.basename(sys.argv[0])
-    domain.filename, _ = os.path.splitext(base)
+import sys, os
+base = os.path.basename(sys.argv[0])
+domain.filename, _ = os.path.splitext(base)
 
-    #LS code to be included in set_quantity
-    from pyvolution import util, least_squares
-    points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
+#LS code to be included in set_quantity
+print 'Reading points'
+from pyvolution import util, least_squares
+points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
+elevation = attributes['elevation']
 
-    #Fit attributes to mesh
-    vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
-                                                  domain.triangles,
-                                                  points,
-                                                  attributes['elevation'],
-                                                  alpha = 0.0001,
-                                                  verbose=True)
+#Fit attributes to mesh
+#vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
+#                                              domain.triangles,
+#                                              points,
+#                                              attributes['elevation'],
+#                                              alpha = 0.0001,
+#                                              verbose=True)
 
-    print 'Initial values'
-    domain.set_quantity('elevation', vertex_attributes)
-    domain.set_quantity('friction', 0.0)
-    domain.set_quantity('stage', 0.0)
+#vertex_attributes = cache(least_squares.fit_to_mesh,
+#                          (domain.coordinates, domain.triangles,
+#                           points, elevation) ,
+#                          {'alpha': 0.0001, 'verbose': True})
 
+print 'Initial values'
+#domain.set_quantity('elevation', vertex_attributes)
+domain.set_quantity('elevation', points=points, attributes=attributes['elevation'], alpha=0.0001)
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('stage', 0.0)
 
-    fid = open('domain.pck', 'w')
-    dump(domain, fid)
-    fid.close()
 
 
@@ -166,5 +162,4 @@
 
 #Set boundary conditions
-
 if read_mesh is True:
     domain.set_boundary({'wave': Bts, 'wall': Br})	
