Changeset 1909

Timestamp:

Oct 13, 2005, 6:10:41 PM (19 years ago)

Author:

steve

Message:

 

Location:

inundation

Files:

• analytical solutions/Analytical solution_oblique_shock_01.py (added)
• pyvolution/mesh_factory.py (modified) (1 diff)

inundation/pyvolution/mesh_factory.py

@@ -676 +676 @@
 
     return points, elements, boundary
+
+
+
+
+def oblique_cross(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] = len(points)
+            points.append([x + origin[0], y + 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
+            v5 = len(points)
+            points.append([x, y])
+
+            #Update boundary dictionary and create elements
+                        #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