Changeset 1742 for inundation/validation/LWRU2/lwru2.py

Timestamp:

Aug 23, 2005, 2:56:21 PM (20 years ago)

Author:

ole

Message:

Recovered validation stuff from 18 August + wrapping up

File:

• inundation/validation/LWRU2/lwru2.py (modified) (2 diffs)

inundation/validation/LWRU2/lwru2.py

@@ -77 +77 @@
 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 +162 @@
 
 #Set boundary conditions
-
 if read_mesh is True:
     domain.set_boundary({'wave': Bts, 'wall': Br})