Changeset 322

Timestamp:

Sep 20, 2004, 10:18:06 AM (20 years ago)

Author:

duncan

Message:

mesh inherits from general_mesh

Location:

inundation/ga/storm_surge/pyvolution

Files:

• general_mesh.py (modified) (3 diffs)
• mesh.py (modified) (7 diffs)
• test_mesh.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -69 +69 @@
         
         self.normals = zeros((N, 6), Float)
+        self.areas = zeros(N, Float)
+        self.edgelengths = zeros((N, 3), Float)
 
         #Get x,y coordinates for all triangles and store
@@ -82 +84 @@
 
             #Area
-            area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
 
             msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
-            msg += ' is degenerate:  area == %f' %area
-            assert area > 0.0, msg
+            msg += ' is degenerate:  area == %f' %self.areas[i]
+            assert self.areas[i] > 0.0, msg
         
 
@@ -117 +119 @@
                                   n1[1], -n1[0],
                                   n2[1], -n2[0]]
+            
+            #Edgelengths
+            self.edgelengths[i, :] = [l0, l1, l2]
             
     def __len__(self):

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -5 +5 @@
 """
 
-
-class Mesh:
+from general_mesh import General_Mesh
+
+class Mesh(General_Mesh):
     """Collection of triangular elements (purely geometric)
 
@@ -65 +66 @@
         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
 
-        self.triangles = array(triangles).astype(Int)
-        self.coordinates = array(coordinates)
-
-        #Input checks
-        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
-        assert len(self.triangles.shape) == 2, msg
-
-        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
-        assert len(self.coordinates.shape) == 2, msg        
-
-        msg = 'Vertex indices reference non-existing coordinate sets'
-        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
-
-
-        #Register number of elements (N)
-        self.number_of_elements = N = self.triangles.shape[0]
-
+        General_Mesh.__init__(self,coordinates, triangles)
+
+        N = self.number_of_elements
+        
         #Allocate space for geometric quantities
         self.centroid_coordinates = zeros((N, 2), Float)
 
-        
-        self.areas = zeros(N, Float)
         self.radii = zeros(N, Float)
-        self.edgelengths = zeros((N, 3), Float)
 
         self.neighbours = zeros((N, 3), Int)
@@ -95 +80 @@
         self.surrogate_neighbours = zeros((N, 3), Int)        
         
-        self.normals = zeros((N, 6), Float)
-
         #Get x,y coordinates for all triangles and store
-        self.vertex_coordinates = V = self.compute_vertex_coordinates()
+        V = self.vertex_coordinates
         
         #Initialise each triangle
@@ -112 +95 @@
             self.centroid_coordinates[i] = centroid
 
-            #Area
-            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-
-            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
-            msg += ' is degenerate:  area == %f' %self.areas[i]
-            assert self.areas[i] > 0.0, msg
-        
             #Midpoints        
             m0 = array([(x1 + x2)/2, (y1 + y2)/2])
@@ -131 +107 @@
             d2 = sqrt(sum( (centroid-m2)**2 ))
         
-            self.radii[i] = min(d0, d1, d2)
-
-            #Normals
-            #The normal vectors
-            # - point outward from each edge
-            # - are orthogonal to the edge
-            # - have unit length
-            # - Are enumerated according to the opposite corner:
-            #   (First normal is associated with the edge opposite
-            #    the first vertex, etc)
-            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
-
-            n0 = array([x2 - x1, y2 - y1])
-            l0 = sqrt(sum(n0**2))
-
-            n1 = array([x0 - x2, y0 - y2])
-            l1 = sqrt(sum(n1**2))
-
-            n2 = array([x1 - x0, y1 - y0])
-            l2 = sqrt(sum(n2**2))
-
-            #Normalise
-            n0 /= l0
-            n1 /= l1
-            n2 /= l2
-
-            #Compute and store
-            self.normals[i, :] = [n0[1], -n0[0],
-                                  n1[1], -n1[0],
-                                  n2[1], -n2[0]]             
-
-            #Edgelengths
-            self.edgelengths[i, :] = [l0, l1, l2]        
+            self.radii[i] = min(d0, d1, d2)               
 
             #Initialise Neighbours (-1 means that it is a boundary neighbour)
@@ -189 +133 @@
 
         
-    def __len__(self):
-        return self.number_of_elements
-
     def __repr__(self):
         return 'Mesh: %d triangles, %d elements, %d boundary segments'\
@@ -498 +439 @@
         #for id, edge in self.boundary:
         #    assert self.neighbours[id,edge] < 0
-        
-
-    
-
-    def get_normals(self):
-        """Return all normal vectors.
-        Return normal vectors for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
-        """
-        return self.normals
-
-
-    def get_normal(self, i, j):
-        """Return normal vector j of the i'th triangle.
-        Return value is the numeric array slice [x, y]
-        """
-        return self.normals[i, 2*j:2*j+2]     
-    
-
-            
-    def get_vertex_coordinates(self):
-        """Return all vertex coordinates.
-        Return all vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
-        """
-        return self.vertex_coordinates
-
-
-    def get_vertex_coordinate(self, i, j):
-        """Return coordinates for vertex j of the i'th triangle.
-        Return value is the numeric array slice [x, y]
-        """
-        return self.vertex_coordinates[i, 2*j:2*j+2]     
-
-
-    def compute_vertex_coordinates(self):
-        """Return vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-        """
-
-        #FIXME: Perhaps they should be ordered as in obj files??
-        #See quantity.get_vertex_values
-        
-        from Numeric import zeros, Float
-
-        N = self.number_of_elements
-        vertex_coordinates = zeros((N, 6), Float) 
-
-        for i in range(N):
-            for j in range(3):
-                k = self.triangles[i,j]  #Index of vertex 0
-                v_k = self.coordinates[k]
-                vertex_coordinates[i, 2*j+0] = v_k[0]
-                vertex_coordinates[i, 2*j+1] = v_k[1]
-
-        return vertex_coordinates
-  
-    def get_triangles(self, unique=False):
-        """Get connectivity
-        If unique is True give them only once as stored internally.
-        If unique is False return
-        """
-
-        if unique is True:
-            return self.triangles
-        else:
-            from Numeric import reshape, array, Int
-                
-            m = len(self)  #Number of volumes
-            M = 3*m        #Total number of unique vertices
-            return reshape(array(range(M)).astype(Int), (m,3))
-
-
 
 #FIXME: May get rid of

inundation/ga/storm_surge/pyvolution/test_mesh.py

@@ -1 +1 @@
 #!/usr/bin/env python
+
+ #FIXME: Seperate the tests for mesh and general_mesh
 
 import unittest
@@ -48 +50 @@
         #Area
         assert mesh.areas[0] == 6.0,\
-               'Area was %f, should have been 6.0' %triangle.get_area()        
+               'Area was %f, should have been 6.0' %mesh.areas[0]
 
         #Normals