Changeset 541

Timestamp:

Nov 15, 2004, 1:51:59 PM (20 years ago)

Author:

ole

Message:

Got oblique analytical solution to work

Location:

inundation/ga/storm_surge

Files:

• analytical solutions/oblique_test_chris.py (modified) (4 diffs)
• analytical solutions/visualise2_chris.py (deleted)
• pyvolution/mesh_factory.py (modified) (2 diffs)

inundation/ga/storm_surge/analytical solutions/oblique_test_chris.py

@@ -12 +12 @@
 # Module imports 
 #
+
+#Were these used?
+#import visualise2_chris as visualise
+#import Image, ImageGrab
+
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+
+from shallow_water import Domain, Constant_height
 from shallow_water import Transmissive_boundary, Reflective_boundary,\
      Dirichlet_boundary
-from domain import oblique_domain, rectangular_domain
 
-from shallow_water import Friction, Constant_height
-import visualise2_chris as visualise
-import Image, ImageGrab
+from math import sqrt, cos, sin, pi
+from mesh_factory import oblique
+
 
 ######################
@@ -28 +38 @@
 lenx = 40.
 
-domain = oblique_domain(m, n, lenx, leny )
+points, elements, boundary = oblique(m, n, lenx, leny)
+domain = Domain(points, elements, boundary) 
 
 # Order of solver
@@ -51 +62 @@
     return 0*x
 
-domain.set_field_values(x_slope, Friction(0))
+domain.set_quantity('elevation', x_slope)
+domain.set_quantity('friction', 0.0)
 
 ######################
 # Boundary conditions
 #
-R = Reflective_boundary()
-T = Transmissive_boundary()
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
 D = Dirichlet_boundary([1.0, 8.57, 0.0])
 
@@ -65 +77 @@
 #Initial condition
 h = 0.5
-domain.set_conserved_quantities(Constant_height(x_slope, h), None, None)
+domain.set_quantity('level', Constant_height(x_slope, h) )
+
 
 ######################
-# Visualisation using Visual Python, comment out if not required
-#visualise.create_surface(domain)
-  
-######################
 #Evolution
-#raw_input('Press key to start simulation')
 import time
 t0 = time.time()
-for t in domain.evolve(max_timestep = 0.1, yieldstep = 0.5, finaltime = 50):
+for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
     domain.write_time()
-
-# Remove the following if Visual Python not required    
-#    visualise.update(domain, index=0, smooth=True)
     
 print 'That took %.2f seconds' %(time.time()-t0)
 
-print "saving file?"
-im = ImageGrab.grab()
-im.save("ccube.eps")
+#print "saving file?"
+#im = ImageGrab.grab()
+#im.save("ccube.eps")
 
 

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -63 +63 @@
 
 
-# def oblique_mesh(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0),
-#                   point_class=Point, element_class=Triangle, mesh_class=Mesh):
-#     """Setup a oblique grid of triangles
-#     with m segments in the x-direction and n segments in the y-direction
-
-#     Triangle refers to the actual class or subclass to be instantiated:
-#     e.g. if Volume is a subclass of Triangle,
-#     this function can be invoked with the keywords
-#       oblique_mesh(...,Triangle=Volume, Mesh=Domain)
-#     """
-
-#     from Numeric import array
-#     from visual import rate
-#     import math
-
-#     from config import epsilon
-
-#     E = m*n*2        #Number of triangular elements
-#     P = (m+1)*(n+1)  #Number of initial vertices
-
-#     initialise_consecutive_datastructure(P+E, E,
-#                                          point_class,
-#                                          element_class,
-#                                          mesh_class)
-
-#     deltax = lenx/float(m)
-#     deltay = leny/float(n)
-#     a = 0.75*lenx*math.tan(theta/180.*math.pi)
-#     x1 = lenx
-#     y1 = 0
-#     x2 = lenx
-#     y2 = leny
-#     x3 = 0.25*lenx
-#     y3 = leny
-#     x4 = x3
-#     y4 = 0
-#     a2 = a/(x1-x4)
-#     a1 = -a2*x4
-#     a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
-#     a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
-
-#     # Dictionary of vertex objects
-#     vertices = {}
-
-#     for i in range(m+1):
-#         x = deltax*i
-#         for j in range(n+1):     
-#             y = deltay*j
-#             if x > 0.25*lenx:
-#                 y = a1 + a2*x + a3*y + a4*x*y
+def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
+    """Setup a oblique grid of triangles
+    with m segments in the x-direction and n segments in the y-direction
+
+    """
+
+    from Numeric import array
+    import math
+
+    from config import epsilon
+
+
+    deltax = lenx/float(m)
+    deltay = leny/float(n)
+    a = 0.75*lenx*math.tan(theta/180.*math.pi)
+    x1 = lenx
+    y1 = 0
+    x2 = lenx
+    y2 = leny
+    x3 = 0.25*lenx
+    y3 = leny
+    x4 = x3
+    y4 = 0
+    a2 = a/(x1-x4)
+    a1 = -a2*x4
+    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
+    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
+
+    # Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        x = deltax*i
+        for j in range(n+1):     
+            y = deltay*j
+            if x > 0.25*lenx:
+                y = a1 + a2*x + a3*y + a4*x*y
                 
-#             vertices[i,j] = Point(x + origin[0],y + origin[1])
-
-#     # Construct 2 elements per element        
-#     elements = []
-#     for i in range(m):
-#         for j in range(n):
-#             v1 = vertices[i,j+1]
-#             v2 = vertices[i,j]            
-#             v3 = vertices[i+1,j+1]            
-#             v4 = vertices[i+1,j]
- 
-#             elements.append(element_class(v4,v3,v2)) #Lower                
-#             elements.append(element_class(v1,v2,v3)) #Upper    
-
-#     M = mesh_class(elements)
-
-#     #Set a default tagging
-
+            vertices[i,j] = len(points)
+            points.append([x + origin[0], y + origin[1]])
+
+    # Construct 2 triangles per element        
+    elements = []
+    boundary = {}
+    for i in range(m):
+        for j in range(n):
+            v1 = vertices[i,j+1]
+            v2 = vertices[i,j]            
+            v3 = vertices[i+1,j+1]            
+            v4 = vertices[i+1,j]
+            
+            #Update boundary dictionary and create elements
+            if i == m-1:
+                boundary[(len(elements), 2)] = 'right'
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'                
+            elements.append([v4,v3,v2]) #Lower
+
+            if i == 0:
+                boundary[(len(elements), 2)] = 'left'
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top' 
+             
+            elements.append([v1,v2,v3]) #Upper
+
+    return points, elements, boundary       
+
+#OLD FORMULA FOR SETTING BOUNDARY TAGS - OBSOLETE NOW             
 #     for id, face in M.boundary:
 
@@ -159 +159 @@
 #             print face, id, id%m, m, n
 #             raise 'Unknown boundary'          
-    
+#   
 #     return M
+
+
 
 def from_polyfile(name):