source: inundation/pyvolution/general_mesh.py @ 2532

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

Added mesh statistics

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