source: storm_surge/pyvolution/general_mesh.py @ 1

Last change on this file since 1 was 1, checked in by duncan, 20 years ago

initial import

File size: 6.5 KB
Line 
1   
2class General_Mesh:
3    """Collection of triangular elements (purely geometric)
4
5    A triangular element is defined in terms of three vertex ids,
6    ordered counter clock-wise,
7    each corresponding to a given coordinate set.
8    Vertices from different elements can point to the same
9    coordinate set.
10
11    Coordinate sets are implemented as an N x 2 Numeric array containing
12    x and y coordinates.
13   
14
15    To instantiate:
16       Mesh(coordinates, triangles)
17
18    where
19
20      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
21      floats representing all x, y coordinates in the mesh.
22
23      triangles is either a list of 3-tuples or an Nx3 Numeric array of
24      integers representing indices of all vertices in the mesh.
25      Each vertex is identified by its index i in [0, M-1].
26
27
28    Example:
29        a = [0.0, 0.0]
30        b = [0.0, 2.0]
31        c = [2.0,0.0]
32        e = [2.0, 2.0]
33       
34        points = [a, b, c, e]
35        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
36        mesh = Mesh(points, triangles)
37   
38        #creates two triangles: bac and bce
39
40        This is a cut down version of mesh from pyvolution\mesh.py
41    """
42
43    def __init__(self, coordinates, triangles):
44        """
45        Build triangles from x,y coordinates (sequence of 2-tuples or
46        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
47        or Nx3 Numeric array of non-negative integers).
48        """
49
50        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
51
52        self.triangles = array(triangles).astype(Int)
53        self.coordinates = array(coordinates)
54
55        #Input checks
56        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
57        assert len(self.triangles.shape) == 2, msg
58
59        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
60        assert len(self.coordinates.shape) == 2, msg       
61
62        msg = 'Vertex indices reference non-existing coordinate sets'
63        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
64
65
66        #Register number of elements (N)
67        self.number_of_elements = N = self.triangles.shape[0]
68   
69        #Allocate space for geometric quantities       
70        self.normals = zeros((N, 6), Float)
71        self.areas = zeros(N, Float)
72        self.edgelengths = zeros((N, 3), Float)
73
74        #Get x,y coordinates for all triangles and store
75        self.vertex_coordinates = V = self.compute_vertex_coordinates()
76       
77        #Initialise each triangle
78        for i in range(N):
79            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
80           
81            x0 = V[i, 0]; y0 = V[i, 1]
82            x1 = V[i, 2]; y1 = V[i, 3]
83            x2 = V[i, 4]; y2 = V[i, 5]           
84
85            #Area
86            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
87
88            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
89            msg += ' is degenerate:  area == %f' %self.areas[i]
90            assert self.areas[i] > 0.0, msg
91       
92
93            #Normals
94            #The normal vectors
95            # - point outward from each edge
96            # - are orthogonal to the edge
97            # - have unit length
98            # - Are enumerated according to the opposite corner:
99            #   (First normal is associated with the edge opposite
100            #    the first vertex, etc)
101            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
102
103            n0 = array([x2 - x1, y2 - y1])
104            l0 = sqrt(sum(n0**2))
105
106            n1 = array([x0 - x2, y0 - y2])
107            l1 = sqrt(sum(n1**2))
108
109            n2 = array([x1 - x0, y1 - y0])
110            l2 = sqrt(sum(n2**2))
111
112            #Normalise
113            n0 /= l0
114            n1 /= l1
115            n2 /= l2
116
117            #Compute and store
118            self.normals[i, :] = [n0[1], -n0[0],
119                                  n1[1], -n1[0],
120                                  n2[1], -n2[0]]
121           
122            #Edgelengths
123            self.edgelengths[i, :] = [l0, l1, l2]
124           
125    def __len__(self):
126        return self.number_of_elements
127
128    def __repr__(self):
129        return 'Mesh: %d triangles, %d elements'\
130               %(self.coordinates.shape[0], len(self))
131
132    def get_normals(self):
133        """Return all normal vectors.
134        Return normal vectors for all triangles as an Nx6 array
135        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
136        """
137        return self.normals
138
139
140    def get_normal(self, i, j):
141        """Return normal vector j of the i'th triangle.
142        Return value is the numeric array slice [x, y]
143        """
144        return self.normals[i, 2*j:2*j+2]     
145   
146
147           
148    def get_vertex_coordinates(self):
149        """Return all vertex coordinates.
150        Return all vertex coordinates for all triangles as an Nx6 array
151        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
152        """
153        return self.vertex_coordinates
154
155
156    def get_vertex_coordinate(self, i, j):
157        """Return coordinates for vertex j of the i'th triangle.
158        Return value is the numeric array slice [x, y]
159        """
160        return self.vertex_coordinates[i, 2*j:2*j+2]     
161
162
163    def compute_vertex_coordinates(self):
164        """Return vertex coordinates for all triangles as an Nx6 array
165        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
166        """
167
168        #FIXME: Perhaps they should be ordered as in obj files??
169        #See quantity.get_vertex_values
170       
171        from Numeric import zeros, Float
172
173        N = self.number_of_elements
174        vertex_coordinates = zeros((N, 6), Float) 
175
176        for i in range(N):
177            for j in range(3):
178                k = self.triangles[i,j]  #Index of vertex 0
179                v_k = self.coordinates[k]
180                vertex_coordinates[i, 2*j+0] = v_k[0]
181                vertex_coordinates[i, 2*j+1] = v_k[1]
182
183        return vertex_coordinates
184 
185    def get_triangles(self, unique=False):
186        """Get connectivity
187        If unique is True give them only once as stored internally.
188        If unique is False return
189        """
190
191        if unique is True:
192            return self.triangles
193        else:
194            from Numeric import reshape, array, Int
195               
196            m = len(self)  #Number of volumes
197            M = 3*m        #Total number of unique vertices
198            return reshape(array(range(M)).astype(Int), (m,3))
Note: See TracBrowser for help on using the repository browser.