source: inundation/pyvolution/general_mesh.py @ 2866

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

Implemented Domain(meshfile) variant widely, updated documentation and added more statistics

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