Changeset 540

Timestamp:

Nov 15, 2004, 1:31:53 PM (19 years ago)

Author:

ole

Message:

Got Smapson to work

Location:

inundation/ga/storm_surge

Files:

• analytical solutions/Analytical_solution_Sampson_import_mesh.py (modified) (4 diffs)
• analytical solutions/Analytical_solution_Yoon_import_mesh.py (modified) (2 diffs)
• pyvolution/mesh_factory.py (modified) (2 diffs)

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

@@ -23 +23 @@
 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
-
-
-
-######################
-# Boundary conditions
-# Not required for this problem
-
-from domain import Strang_domain
+from mesh_factory import strang_mesh
 
 
@@ -43 +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
 domain.smooth = True
 
+# Provide file name for storing output
+domain.store = True
+domain.format = 'sww'
+domain.filename = 'sampson_strang_second_order'
 
-# Provide file name for storing output
-domain.filename="yoon_strang_second_order"
+print 'Number of triangles = ', len(domain)
 
-print "Number of triangles = ", len(Volume.instances)
-
-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
+from util import mean
 domain.reduction = mean 
 #domain.reduction = min  #Looks better near steep slopes
@@ -91 +80 @@
     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]
@@ -107 +95 @@
     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/analytical solutions/Analytical_solution_Yoon_import_mesh.py

@@ -20 +20 @@
 from shallow_water import Domain, Dirichlet_boundary
 from math import sqrt, cos, sin, pi
-
-
-
-######################
-# Boundary conditions
-# Not required for this problem
-
 from mesh_factory import strang_mesh
 
@@ -45 +38 @@
 domain.store = True
 domain.format = 'sww'
-domain.filename='yoon_strang_second_order'
+domain.filename = 'yoon_strang_second_order'
 
 print 'Number of triangles = ', len(domain)

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -62 +62 @@
     return points, elements, boundary       
 
-
-
-def from_polyfile(name):
-    """Read mesh from .poly file, an obj like file format
-    listing first vertex coordinates and then connectivity
-    """
-
-    from util import anglediff
-    from math import pi
-    import os.path
-    root, ext = os.path.splitext(name)
-
-    if ext == 'poly':
-        filename = name
-    else:
-        filename = name + '.poly'
-
-    
-    fid = open(filename)
-
-    points = []    #x, y
-    values = []    #z
-    ##vertex_values = []    #Repeated z    
-    triangles = [] #v0, v1, v2
-    
-    lines = fid.readlines()
-
-    keyword = lines[0].strip()
-    msg = 'First line in .poly file must contain the keyword: POINTS'
-    assert keyword == 'POINTS', msg
-
-    offending = 0
-    i = 1
-    while keyword == 'POINTS':
-        line = lines[i].strip()
-        i += 1
-
-        if line == 'POLYS':
-            keyword = line
-            break
-
-        fields = line.split(':')
-        assert int(fields[0]) == i-1, 'Point indices not consecutive'
-
-        #Split the three floats
-        xyz = fields[1].split()
-
-        x = float(xyz[0])
-        y = float(xyz[1])
-        z = float(xyz[2])        
-
-        points.append([x, y])
-        values.append(z)
-        
-
-    k = i    
-    while keyword == 'POLYS':
-        line = lines[i].strip()
-        i += 1
-
-        if line == 'END':
-            keyword = line
-            break
-        
-
-        fields = line.split(':')
-        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
-
-        #Split the three indices
-        vvv = fields[1].split()
-
-        i0 = int(vvv[0])-1
-        i1 = int(vvv[1])-1
-        i2 = int(vvv[2])-1
-
-        #Check for and exclude degenerate areas
-        x0 = points[i0][0]
-        y0 = points[i0][1]
-        x1 = points[i1][0]
-        y1 = points[i1][1]
-        x2 = points[i2][0]
-        y2 = points[i2][1]                
-
-        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-        if area > 0:
-
-            #Ensure that points are arranged in counter clock-wise order
-            v0 = [x1-x0, y1-y0]
-            v1 = [x2-x1, y2-y1]
-            v2 = [x0-x2, y0-y2]
-
-            a0 = anglediff(v1, v0)
-            a1 = anglediff(v2, v1)
-            a2 = anglediff(v0, v2)        
-
-            
-            if a0 < pi and a1 < pi and a2 < pi:
-                #all is well
-                j0 = i0
-                j1 = i1
-                j2 = i2
-            else:
-                #Swap two vertices
-                j0 = i1
-                j1 = i0
-                j2 = i2
-                
-            triangles.append([j0, j1, j2])
-            ##vertex_values.append([values[j0], values[j1], values[j2]])
-        else:
-            offending +=1 
-
-    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),
@@ -388 +161 @@
     
 #     return M
+
+def from_polyfile(name):
+    """Read mesh from .poly file, an obj like file format
+    listing first vertex coordinates and then connectivity
+    """
+
+    from util import anglediff
+    from math import pi
+    import os.path
+    root, ext = os.path.splitext(name)
+
+    if ext == 'poly':
+        filename = name
+    else:
+        filename = name + '.poly'
+
+    
+    fid = open(filename)
+
+    points = []    #x, y
+    values = []    #z
+    ##vertex_values = []    #Repeated z    
+    triangles = [] #v0, v1, v2
+    
+    lines = fid.readlines()
+
+    keyword = lines[0].strip()
+    msg = 'First line in .poly file must contain the keyword: POINTS'
+    assert keyword == 'POINTS', msg
+
+    offending = 0
+    i = 1
+    while keyword == 'POINTS':
+        line = lines[i].strip()
+        i += 1
+
+        if line == 'POLYS':
+            keyword = line
+            break
+
+        fields = line.split(':')
+        assert int(fields[0]) == i-1, 'Point indices not consecutive'
+
+        #Split the three floats
+        xyz = fields[1].split()
+
+        x = float(xyz[0])
+        y = float(xyz[1])
+        z = float(xyz[2])        
+
+        points.append([x, y])
+        values.append(z)
+        
+
+    k = i    
+    while keyword == 'POLYS':
+        line = lines[i].strip()
+        i += 1
+
+        if line == 'END':
+            keyword = line
+            break
+        
+
+        fields = line.split(':')
+        assert int(fields[0]) == i-k, 'Poly indices not consecutive'
+
+        #Split the three indices
+        vvv = fields[1].split()
+
+        i0 = int(vvv[0])-1
+        i1 = int(vvv[1])-1
+        i2 = int(vvv[2])-1
+
+        #Check for and exclude degenerate areas
+        x0 = points[i0][0]
+        y0 = points[i0][1]
+        x1 = points[i1][0]
+        y1 = points[i1][1]
+        x2 = points[i2][0]
+        y2 = points[i2][1]                
+
+        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+        if area > 0:
+
+            #Ensure that points are arranged in counter clock-wise order
+            v0 = [x1-x0, y1-y0]
+            v1 = [x2-x1, y2-y1]
+            v2 = [x0-x2, y0-y2]
+
+            a0 = anglediff(v1, v0)
+            a1 = anglediff(v2, v1)
+            a2 = anglediff(v0, v2)        
+
+            
+            if a0 < pi and a1 < pi and a2 < pi:
+                #all is well
+                j0 = i0
+                j1 = i1
+                j2 = i2
+            else:
+                #Swap two vertices
+                j0 = i1
+                j1 = i0
+                j2 = i2
+                
+            triangles.append([j0, j1, j2])
+            ##vertex_values.append([values[j0], values[j1], values[j2]])
+        else:
+            offending +=1 
+
+    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 contracting_channel_mesh(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0,