Changeset 1761

Timestamp:

Aug 25, 2005, 10:06:22 AM (19 years ago)

Author:

steve

Message:

Added a more symmetric rectangle_cross procedure. Can be used in place of rectangle from mesh_factory

Location:

inundation

Files:

• pyvolution/mesh_factory.py (modified) (4 diffs)
• validation/LWRU2/lwru2.py (modified) (5 diffs)

inundation/pyvolution/mesh_factory.py

@@ -20 +20 @@
     from config import epsilon
 
-    #E = m*n*2        #Number of triangular elements
-    #P = (m+1)*(n+1)  #Number of initial vertices
-
-    delta1 = float(len1)/m
-    delta2 = float(len2)/n
+    deltax = float(len1)/m
+    deltay = float(len2)/n
 
     #Dictionary of vertex objects
@@ -34 +31 @@
             vertices[i,j] = len(points)
             points.append([i*delta1 + origin[0], j*delta2 + origin[1]])
-
 
 
@@ -120 +116 @@
             i4 = index(i+1,j)
 
+            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
+            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
+            v5 = len(points)
+            points.append([x, y])
+
+            #Create left triangle
+            if i == 0:
+                boundary[(len(elements), 1)] = 'left'
+            elements.append([v2,v5,v1])
+
+            #Create bottom triangle
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v5,v2])
+
+            #Create right triangle
+            if i == m-1:
+                boundary[(len(elements), 1)] = 'right'
+            elements.append([v3,v5,v4])
+
+            #Create top triangle
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v5,v3])
+
             #Update boundary dictionary and create elements
             if i == m-1:
@@ -133 +154 @@
                 boundary[nt, 1] = 'top'
             elements[nt,:] = [i1,i2,i3] #Upper element
+
+    return points, elements, boundary
+
+
+def rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+    """Setup a rectangular grid of triangles
+    with m+1 by n+1 grid points
+    and side lengths len1, len2. If side lengths are omitted
+    the mesh defaults to the unit square.
+
+    len1: x direction (left to right)
+    len2: y direction (bottom to top)
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements)
+    """
+
+    from config import epsilon
+    from Numeric import zeros, Float, Int
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n
+
+    #Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        for j in range(n+1):
+            vertices[i,j] = len(points)
+            points.append([delta1*i + origin[0], delta2*j  + origin[1]])
+
+    # Construct 4 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]
+            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
+            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
+
+            # Create centre point
+            v5 = len(points)
+            points.append([x, y])
+
+            #Create left triangle
+            if i == 0:
+                boundary[(len(elements), 1)] = 'left'
+            elements.append([v2,v5,v1])
+
+            #Create bottom triangle
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v5,v2])
+
+            #Create right triangle
+            if i == m-1:
+                boundary[(len(elements), 1)] = 'right'
+            elements.append([v3,v5,v4])
+
+            #Create top triangle
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v5,v3])
+
 
     return points, elements, boundary

inundation/validation/LWRU2/lwru2.py

@@ -74 +74 @@
 from pyvolution.shallow_water import Domain, Reflective_boundary,\
             File_boundary, Transmissive_Momentum_Set_Stage_boundary
-from pyvolution.mesh_factory import rectangular
+from pyvolution.mesh_factory import rectangular_cross
 from pyvolution.pmesh2domain import pmesh_to_domain_instance
 from Numeric import array, zeros, Float, allclose
@@ -91 +91 @@
 #######################
 # Domain
-    
+
 if read_mesh is True:
     print 'Creating domain from', filenames.mesh_filename
@@ -101 +101 @@
 
 
-    
-else:   
+
+else:
     print 'Creating regular mesh'
     N = 100
-    points, vertices, boundary = rectangular(N, N/5*3,
+    points, vertices, boundary = rectangular_cross(N, N/5*3,
                                              len1=5.448, len2=3.402)
     print 'Creating domain'
-    
+
     #domain = Domain(points, vertices, boundary)
     domain = cache(Domain, (points, vertices, boundary))
@@ -135 +135 @@
                     verbose = True,
                     use_cache = True)
-                    
+
 domain.set_quantity('friction', 0.0)
 domain.set_quantity('stage', 0.0)
@@ -157 +157 @@
 #Set boundary conditions
 if read_mesh is True:
-    domain.set_boundary({'wave': Bts, 'wall': Br})      
-else:    
+    domain.set_boundary({'wave': Bts, 'wall': Br})
+else:
     domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})