source: inundation/pyvolution/general_mesh.py @ 1916

Last change on this file since 1916 was 1894, checked in by ole, 18 years ago

Comments and allow for duplicate timesteps

File size: 10.1 KB
Line 
1
2from coordinate_transforms.geo_reference import Geo_reference
3
4class General_mesh:
5    """Collection of triangular elements (purely geometric)
6
7    A triangular element is defined in terms of three vertex ids,
8    ordered counter clock-wise,
9    each corresponding to a given coordinate set.
10    Vertices from different elements can point to the same
11    coordinate set.
12
13    Coordinate sets are implemented as an N x 2 Numeric array containing
14    x and y coordinates.
15
16
17    To instantiate:
18       Mesh(coordinates, triangles)
19
20    where
21
22      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
23      floats representing all x, y coordinates in the mesh.
24
25      triangles is either a list of 3-tuples or an Nx3 Numeric array of
26      integers representing indices of all vertices in the mesh.
27      Each vertex is identified by its index i in [0, M-1].
28
29
30    Example:
31        a = [0.0, 0.0]
32        b = [0.0, 2.0]
33        c = [2.0,0.0]
34        e = [2.0, 2.0]
35
36        points = [a, b, c, e]
37        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
38        mesh = Mesh(points, triangles)
39
40        #creates two triangles: bac and bce
41
42    Other:
43
44      In addition mesh computes an Nx6 array called vertex_coordinates.
45      This structure is derived from coordinates and contains for each
46      triangle the three x,y coordinates at the vertices.
47
48
49        This is a cut down version of mesh from pyvolution mesh.py
50    """
51
52    def __init__(self, coordinates, triangles,
53                 geo_reference=None):
54        """
55        Build triangles from x,y coordinates (sequence of 2-tuples or
56        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
57        or Nx3 Numeric array of non-negative integers).
58
59        origin is a 3-tuple consisting of UTM zone, easting and northing.
60        If specified coordinates are assumed to be relative to this origin.
61        """
62
63        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
64
65        self.triangles = array(triangles).astype(Int)
66        self.coordinates = array(coordinates)
67
68        if geo_reference is None:
69            self.geo_reference = Geo_reference() #Use defaults
70        else:
71            self.geo_reference = geo_reference
72
73        #Input checks
74        msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples'
75        assert len(self.triangles.shape) == 2, msg
76
77        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
78        assert len(self.coordinates.shape) == 2, msg
79
80        msg = 'Vertex indices reference non-existing coordinate sets'
81        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
82
83
84        #Register number of elements (N)
85        self.number_of_elements = N = self.triangles.shape[0]
86
87        xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
88                      max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
89
90        self.xy_extent = array(xy_extent, Float)
91
92
93        #Allocate space for geometric quantities
94        self.normals = zeros((N, 6), Float)
95        self.areas = zeros(N, Float)
96        self.edgelengths = zeros((N, 3), Float)
97
98        #Get x,y coordinates for all triangles and store
99        self.vertex_coordinates = V = self.compute_vertex_coordinates()
100
101        #Initialise each triangle
102        for i in range(N):
103            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
104
105            x0 = V[i, 0]; y0 = V[i, 1]
106            x1 = V[i, 2]; y1 = V[i, 3]
107            x2 = V[i, 4]; y2 = V[i, 5]
108
109            #Area
110            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
111
112            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
113            msg += ' is degenerate:  area == %f' %self.areas[i]
114            assert self.areas[i] > 0.0, msg
115
116
117            #Normals
118            #The normal vectors
119            # - point outward from each edge
120            # - are orthogonal to the edge
121            # - have unit length
122            # - Are enumerated according to the opposite corner:
123            #   (First normal is associated with the edge opposite
124            #    the first vertex, etc)
125            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
126
127            n0 = array([x2 - x1, y2 - y1])
128            l0 = sqrt(sum(n0**2))
129
130            n1 = array([x0 - x2, y0 - y2])
131            l1 = sqrt(sum(n1**2))
132
133            n2 = array([x1 - x0, y1 - y0])
134            l2 = sqrt(sum(n2**2))
135
136            #Normalise
137            n0 /= l0
138            n1 /= l1
139            n2 /= l2
140
141            #Compute and store
142            self.normals[i, :] = [n0[1], -n0[0],
143                                  n1[1], -n1[0],
144                                  n2[1], -n2[0]]
145
146            #Edgelengths
147            self.edgelengths[i, :] = [l0, l1, l2]
148
149        #Build vertex list
150        self.build_vertexlist()
151
152    def __len__(self):
153        return self.number_of_elements
154
155    def __repr__(self):
156        return 'Mesh: %d triangles, %d elements'\
157               %(self.coordinates.shape[0], len(self))
158
159    def get_normals(self):
160        """Return all normal vectors.
161        Return normal vectors for all triangles as an Nx6 array
162        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
163        """
164        return self.normals
165
166
167    def get_normal(self, i, j):
168        """Return normal vector j of the i'th triangle.
169        Return value is the numeric array slice [x, y]
170        """
171        return self.normals[i, 2*j:2*j+2]
172
173
174
175    def get_vertex_coordinates(self, obj = False):
176        """Return all vertex coordinates.
177        Return all vertex coordinates for all triangles as an Nx6 array
178        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
179
180        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
181        FIXME, we might make that the default.
182        FIXME Maybe use keyword: continuous = False for this condition?
183
184       
185        """
186
187        if obj is True:
188            from Numeric import concatenate, reshape
189            V = self.vertex_coordinates
190            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
191
192            N = V.shape[0]
193            return reshape(V, (3*N, 2))
194        else:   
195            return self.vertex_coordinates
196
197
198    def get_vertex_coordinate(self, i, j):
199        """Return coordinates for vertex j of the i'th triangle.
200        Return value is the numeric array slice [x, y]
201        """
202        return self.vertex_coordinates[i, 2*j:2*j+2]
203
204
205    def compute_vertex_coordinates(self):
206        """Return vertex coordinates for all triangles as an Nx6 array
207        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
208        """
209
210        #FIXME (Ole): Perhaps they should be ordered as in obj files??
211        #See quantity.get_vertex_values
212        #FIXME (Ole) - oh yes they should
213
214        from Numeric import zeros, Float
215
216        N = self.number_of_elements
217        vertex_coordinates = zeros((N, 6), Float)
218
219        for i in range(N):
220            for j in range(3):
221                k = self.triangles[i,j]  #Index of vertex 0
222                v_k = self.coordinates[k]
223                vertex_coordinates[i, 2*j+0] = v_k[0]
224                vertex_coordinates[i, 2*j+1] = v_k[1]
225
226        return vertex_coordinates
227
228    def get_vertices(self, indices=None):
229        """Get connectivity
230        indices is the set of element ids of interest
231        """
232
233        from Numeric import take
234
235        if (indices ==  None):
236            indices = range(len(self))  #len(self)=number of elements
237
238        return  take(self.triangles, indices)
239
240    #FIXME - merge these two (get_vertices and get_triangles)
241    def get_triangles(self, obj = False):
242        """Get connetivity
243        Return triangles (triplets of indices into point coordinates)
244       
245        If obj is True return structure commensurate with replicated
246        points, allowing for discontinuities
247        (FIXME: Need good name for this concept)
248        """
249
250        if obj is True:
251            from Numeric import array, reshape, Int       
252            m = len(self)  #Number of triangles
253            M = 3*m        #Total number of unique vertices
254            T = reshape(array(range(M)).astype(Int), (m,3))
255        else:
256            T = self.triangles
257
258        return T     
259
260   
261
262    def get_unique_vertices(self,  indices=None):
263        triangles = self.get_vertices(indices=indices)
264        unique_verts = {}
265        for triangle in triangles:
266           unique_verts[triangle[0]] = 0
267           unique_verts[triangle[1]] = 0
268           unique_verts[triangle[2]] = 0
269        return unique_verts.keys()
270
271    def build_vertexlist(self):
272        """Build vertexlist index by vertex ids and for each entry (point id)
273        build a list of (triangles, vertex_id) pairs that use the point
274        as vertex.
275
276        Preconditions:
277          self.coordinates and self.triangles are defined
278
279        Postcondition:
280          self.vertexlist is built
281        """
282
283        vertexlist = [None]*len(self.coordinates)
284        for i in range(self.number_of_elements):
285
286            a = self.triangles[i, 0]
287            b = self.triangles[i, 1]
288            c = self.triangles[i, 2]
289
290            #Register the vertices v as lists of
291            #(triangle_id, vertex_id) tuples associated with them
292            #This is used for smoothing
293            for vertex_id, v in enumerate([a,b,c]):
294                if vertexlist[v] is None:
295                    vertexlist[v] = []
296
297                vertexlist[v].append( (i, vertex_id) )
298
299        self.vertexlist = vertexlist
300
301
302    def get_extent(self):
303        """Return min and max of all x and y coordinates
304        """
305
306
307        C = self.get_vertex_coordinates()
308        X = C[:,0:6:2].copy()
309        Y = C[:,1:6:2].copy()
310
311        xmin = min(X.flat)
312        xmax = max(X.flat)
313        ymin = min(Y.flat)
314        ymax = max(Y.flat)
315
316        return xmin, xmax, ymin, ymax
317
318
319    def get_area(self):
320        """Return total area of mesh
321        """
322
323        return sum(self.areas)
Note: See TracBrowser for help on using the repository browser.