source: inundation/pyvolution/general_mesh.py @ 3072

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

Made boundary polygon return absolute UTM coordinates as per ticket:81

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