Changeset 534

Timestamp:

Nov 15, 2004, 11:08:25 AM (20 years ago)

Author:

ole

Message:

Migrayed analytical solution to current version of Pyvolution

Location:

inundation/ga/storm_surge

Files:

• analytical solutions/Analytical_solution_Yoon_import_mesh.py (modified) (7 diffs)
• pyvolution/mesh_factory.py (modified) (2 diffs)
• pyvolution/shallow_water.py (modified) (2 diffs)

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

@@ -6 +6 @@
    Geoscience Australia
    
-Specific methods pertaining to the 2D shallow water equation
-are imported from shallow_water
-for use with the generic finite volume framework
-
-Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
-numerical vector named conserved_quantities.
 """
 
@@ -22 +16 @@
 sys.path.append('..'+sep+'pyvolution')
 
-from shallow_water import Transmissive_boundary, Reflective_boundary, \
-     Dirichlet_boundary
-from domain import Domain
-from Numeric import array
+from shallow_water import Domain, Dirichlet_boundary
 from math import sqrt, cos, sin, pi
 
@@ -34 +25 @@
 # Not required for this problem
 
-from domain import Strang_domain
+from mesh_factory import strang_mesh
 
 
@@ -42 +33 @@
 # two or three entries. Two entries are for points and three for triangles.
 
-domain = Strang_domain('yoon_circle.pt')
+points, elements = strang_mesh('yoon_circle.pt')
+domain = Domain(points, elements) 
 
 domain.default_order = 2
@@ -49 +41 @@
 
 # Provide file name for storing output
-domain.filename="yoon_strang_second_order"
+domain.store = True
+domain.format = 'sww'
+domain.filename='yoon_strang_second_order'
 
-print "Number of triangles = ", len(Volume.instances)
+print 'Number of triangles = ', len(domain)
 
-domain.visualise = False
-domain.checkpoint = False #
-domain.store = True    #Store for visualisation purposes
-domain.format = 'sww'   #Native netcdf visualisation format
-
-#Reduction operation for get_vertex_values              
-from pytools.stats import mean
+#Reduction operation for get_vertex_values
+from util import mean
 domain.reduction = mean 
 #domain.reduction = min  #Looks better near steep slopes
@@ -90 +79 @@
     return z
 
-domain.set_field_values(x_slope, None)
+domain.set_quantity('elevation', x_slope)
 
-
-#Set the water depth
-def height(x,y):
+#Set the water level (w = z+h)
+def level(x,y):
     z = x_slope(x,y)
     n = x.shape[0]
@@ -106 +94 @@
     return h   
 
-domain.set_conserved_quantities(height, None, None)
+domain.set_quantity('level', level)
 
+
+############
+#Boundary
+domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
     
-######################
-#Visualisation
-if domain.visualise:
-    visualise.create_surface(domain)
-  
 
 ######################
 #Evolution
-
 import time
 t0 = time.time()
-for t in domain.evolve(max_timestep = 1.0, yieldstep = 10.0, finaltime = 300):
+for t in domain.evolve(yieldstep = 10.0, finaltime = 300):
     domain.write_time()
 
-# Remove the following if Visual Python not required       
-    if domain.visualise:  
-        visualise.update(domain, index=0, smooth=True)
-    
 print 'That took %.2f seconds' %(time.time()-t0)
 

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -62 +62 @@
     return points, elements, boundary       
 
-
-#FIXME: Get oblique and contracting and circular meshes from Chris
 
 
@@ -178 +176 @@
     print 'Removed %d offending triangles out of %d' %(offending, len(lines))
     return points, triangles, values
+
+
+
+def strang_mesh(filename):
+    """Read Strang generated mesh.
+    """
+
+    from math import pi
+    from util import anglediff
+    
+    
+    fid = open(filename)
+    points = []    # List of x, y coordinates
+    triangles = [] # List of vertex ids as listed in the file
+    
+    for line in fid.readlines():
+        fields = line.split()
+        if len(fields) == 2:
+            # we are reading vertex coordinates
+            points.append([float(fields[0]), float(fields[1])])
+        elif len(fields) == 3:
+            # we are reading triangle point id's (format ae+b)
+            triangles.append([int(float(fields[0]))-1,
+                              int(float(fields[1]))-1,
+                              int(float(fields[2]))-1])
+        else:
+            raise 'wrong format in ' + filename
+
+    elements = [] #Final list of elements
+
+    for t in triangles:
+        #Get vertex coordinates
+        v0 = t[0]
+        v1 = t[1]
+        v2 = t[2]
+        
+        x0 = points[v0][0]
+        y0 = points[v0][1]
+        x1 = points[v1][0]
+        y1 = points[v1][1]
+        x2 = points[v2][0]
+        y2 = points[v2][1]                
+
+        #Check that points are arranged in counter clock-wise order
+        vec0 = [x1-x0, y1-y0]
+        vec1 = [x2-x1, y2-y1]
+        vec2 = [x0-x2, y0-y2]
+
+        a0 = anglediff(vec1, vec0)
+        a1 = anglediff(vec2, vec1)
+        a2 = anglediff(vec0, vec2)        
+
+        if a0 < pi and a1 < pi and a2 < pi:
+            elements.append([v0, v1, v2])
+        else:
+            elements.append([v0, v2, v1])       
+
+    return points, elements                             
+
+
+
+#FIXME: These haven't been done yet
+# def circular_mesh(m, n, radius=1.0, center = (0.0, 0.0),  Triangle=Triangle,
+#                   Mesh=Mesh, Point=Point):
+#     """Setup a circular grid of triangles
+#     with m+1 radial segments by n+1 points around the circuference
+#     and radius. If radius is are omitted the mesh defaults to the unit circle
+#     radius.
+ 
+#     radius: radius of circle
+
+#     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
+#     rectangular_mesh(...,Triangle=Volume, Mesh=Domain)"""
+    
+
+#     from Numeric import array
+#     from visual import rate
+#     import math
+
+#     delta = radius/float(m)
+
+#     #Dictionary of vertex objects
+#     vertices = {}
+#     for i in range(n+1):
+#         theta = float(i)*(2*math.pi)/float(n)  
+#         for j in range(m+1):
+#             delta = float(j)*radius/float(m)
+#             vertices[i,j] = Point(delta*math.cos(theta),delta*math.sin(theta))
+
+#     #Construct 2 triangles per element        
+#     elements = []
+#     for i in range(n):
+#         for j in range(1,m):
+#             v1 = vertices[i,j+1]
+#             v2 = vertices[i,j]            
+#             v3 = vertices[i+1,j+1]            
+#             v4 = vertices[i+1,j]
+
+#             elements.append(Triangle(v4,v2,v3)) #Lower                
+#             elements.append(Triangle(v1,v3,v2)) #Upper    
+
+#     #Do the center
+#     v1 = vertices[0,0]
+#     for i in range(n):
+#         v2 = vertices[i,1]
+#         v3 = vertices[i+1,1]
+
+#         T = Triangle(v1,v2,v3) #center
+#         elements.append(T)
+
+#     return Mesh(elements)
+
+# 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
+                
+#             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
+
+#     for id, face in M.boundary:
+
+#         e = element_class.instances[id]
+#         x0, y0, x1, y1, x2, y2 = e.get_instance_vertex_coordinates()
+
+#         if face==2:
+#             #Left or right#
+#             if abs(x0-origin[0]) < epsilon:
+#                 M.boundary[(id,face)] = 'left'            
+#             elif abs(origin[0]+lenx-x0) < epsilon:
+#                 M.boundary[(id,face)] = 'right'
+#             else:
+#                 print face, id, id%m, m, n
+#                 raise 'Left or Right Unknown boundary'                            
+#         elif face==1:
+#             #Top or bottom
+#             if x0 > 0.25*lenx and abs(y0-a1-a2*x0-origin[1]) < epsilon or\
+#                x0 <= 0.25*lenx and abs(y0-origin[1]) < epsilon:            
+#                    M.boundary[(id,face)] = 'bottom'
+#             elif abs(origin[1]+leny-y0) < epsilon:
+#                 M.boundary[(id,face)] = 'top'
+#             else:
+#                 print face, id, id%m, m, n
+#                 raise 'Top or Bottom Unknown boundary'  
+#         else:
+#             print face, id, id%m, m, n
+#             raise 'Unknown boundary'          
+    
+#     return M
+
+# def contracting_channel_mesh(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0,
+#                             y1 =0.0, y4 = -1./4., y8 = 1.0,
+#                              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 = (x4 - x1)/float(m)
+#     deltay = (y8 - y1)/float(n)
+#     a = y4 - y1
+    
+#     if y8 - a <= y1 + a:
+#         print a,y1,y4
+#         raise 'Constriction is too large reduce a'
+            
+#     y2 = y1
+#     y3 = y4
+#     x5 = x4
+#     y5 = y8 - a
+#     x6 = x3
+#     y6 = y5
+#     x7 = x2
+#     y7 = y8
+#     x8 = x1
+    
+#     a2 = a/(x3 - x2)
+#     a1 = a - a*x3/(x3 - x2)
+#     a4 = (-a + a2*(x7 - x6))/(x6 - x7)/y7
+#     a3 = (y7 - a1 - x7*a2 - a4*x7*y7)/y7
+    
+#     # 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 > x2:
+#                 if x < x3:
+#                     y = a1 + a2*x + a3*y + a4*x*y
+#                 else:
+#                     y = a + y*(y5 - y4)/(y8 - y1)   
+                
+#             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
+
+#     for id, face in M.boundary:
+
+#         e = element_class.instances[id]
+#         x_0, y_0, x_1, y_1, x_2, y_2 = e.get_instance_vertex_coordinates()
+#         lenx = x4 - x1
+#         if face==2:
+#             #Left or right#
+#             if abs(x_0-origin[0]) < epsilon:
+#                 M.boundary[(id,face)] = 'left'            
+#             elif abs(origin[0]+lenx-x_0) < epsilon:
+#                 M.boundary[(id,face)] = 'right'
+#             else:
+#                 print face, id, id%m, m, n
+#                 raise 'Left or Right Unknown boundary'                            
+#         elif face==1:
+#             #Top or bottom
+#             if x_0 <= x2 and abs(y_0-y1-origin[1]) < epsilon or\
+#                x_0 > x3 and abs(y_0-y3-origin[1]) < epsilon or\
+#                x_0 >  x2 and x_0 <= x3 and abs(y_0-(y2+(y3-y2)*(x_0-x2)/(x3-x2)+origin[1])) < epsilon:            
+#                 M.boundary[(id,face)] = 'bottom'
+#             elif x_0 <= x7 and abs(y_0-y8-origin[1]) < epsilon or\
+#                x_0 > x6 and abs(y_0-y6-origin[1]) < epsilon or\
+#                x_0 > x7 and x_0 <= x6 and abs(y_0-(y7+(y6-y7)*(x_0-x7)/(x6-x7)+origin[1])) < epsilon:            
+#                 M.boundary[(id,face)] = 'top'
+#             else:
+#                 print face, id, id%m, m, n 
+#                 raise 'Top or Bottom Unknown boundary'  
+#         else:
+#             print face, id, id%m, m, n
+#             raise 'Unknown boundary'          
+
+            
+#       #  print id, face, M.boundary[(id,face)],x_0,y_0,x_1,y_1,x_2,y_2
+    
+#     return M
+
+
+
+# #Map from edge number to indices of associated vertices
+# edge_map = ((1,2), (0,2), (0,1))
+
+
+                                         
+    

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -48 +48 @@
 
         #Reduction operation for get_vertex_values              
-        #from pytools.stats import mean
+        #from util import mean
         #self.reduction = mean 
         self.reduction = min  #Looks better near steep slopes
@@ -155 +155 @@
         #Store model data, e.g. for visualisation    
         if self.store is True and self.time == 0.0:
-            self.initialise_storage()   
+            self.initialise_storage()
+            print 'Storing results in ' + self.writer.filename
+        else:
+            print 'Results will not be stored.'
+            print 'To store results set domain.store = True'