Changeset 317 for inundation/ga/storm_surge/pyvolution/mesh_factory.py

Timestamp:

Sep 16, 2004, 5:48:27 PM (20 years ago)

Author:

ole

Message:

Started playing with the Cornell Room

File:

• inundation/ga/storm_surge/pyvolution/mesh_factory.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -64 +64 @@
 
 
+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, vertex_values