# source:anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py@3956

Last change on this file since 3956 was 3956, checked in by ole, 17 years ago

More cleanup and refactoring. Also adhered to style guide for comparisons to
singletons such as None.

File size: 13.8 KB
Line
1
2from Numeric import concatenate, reshape, take, allclose
3from Numeric import array, zeros, Int, Float, sqrt, sum
4
5from anuga.coordinate_transforms.geo_reference import Geo_reference
6
7class General_mesh:
8    """Collection of 2D triangular elements
9
10    A triangular element is defined in terms of three vertex ids,
11    ordered counter clock-wise, each corresponding to a given node
12    which is represented as a coordinate set (x,y).
13    Vertices from different triangles can point to the same node.
14    The nodes are implemented as an Nx2 Numeric array containing the
15    x and y coordinates.
16
17
18    To instantiate:
19       Mesh(nodes, triangles)
20
21    where
22
23      nodes is either a list of 2-tuples or an Nx2 Numeric array of
24      floats representing all x, y coordinates in the mesh.
25
26      triangles is either a list of 3-tuples or an Mx3 Numeric array of
27      integers representing indices of all vertices in the mesh.
28      Each vertex is identified by its index i in [0, N-1].
29
30
31    Example:
32
33        a = [0.0, 0.0]
34        b = [0.0, 2.0]
35        c = [2.0,0.0]
36        e = [2.0, 2.0]
37
38        nodes = [a, b, c, e]
39        triangles = [ [1,0,2], [1,2,3] ]   # bac, bce
40
41        # Create mesh with two triangles: bac and bce
42        mesh = Mesh(nodes, triangles)
43
44
45
46    Other:
47
48      In addition mesh computes an Mx6 array called vertex_coordinates.
49      This structure is derived from coordinates and contains for each
50      triangle the three x,y coordinates at the vertices.
51
52      See neighbourmesh.py for a specialisation of the general mesh class
53      which includes information about neighbours and the mesh boundary.
54
55      The mesh object is purely geometrical and contains no information
56      about quantities defined on the mesh.
57
58    """
59
60    #FIXME: It would be a good idea to use geospatial data as an alternative
61    #input
62    def __init__(self, nodes, triangles,
63                 geo_reference=None,
64                 number_of_full_nodes=None,
65                 number_of_full_triangles=None,
66                 verbose=False):
67        """Build triangular 2d mesh from nodes and triangle information
68
69        Input:
70
71          nodes: x,y coordinates represented as a sequence of 2-tuples or
72                 a Nx2 Numeric array of floats.
73
74          triangles: sequence of 3-tuples or Mx3 Numeric array of
75                     non-negative integers representing indices into
76                     the nodes array.
77
78          georeference (optional): If specified coordinates are
79          assumed to be relative to this origin.
80
81
82        number_of_full_nodes and number_of_full_triangles relate to
83        parallelism when each mesh has an extra layer of ghost points and
84        ghost triangles attached to the end of the two arrays.
85        In this case it is usefull to specify the number of real (called full)
86        nodes and triangles. If omitted they will default to all.
87
88        """
89
90        if verbose: print 'General_mesh: Building basic mesh structure'
91
92        self.triangles = array(triangles, Int)
93        self.nodes = array(nodes, Float)
94
95
96        # Register number of elements and nodes
97        self.number_of_triangles = N = self.triangles.shape[0]
98        self.number_of_nodes = self.nodes.shape[0]
99
100
101
102        if number_of_full_nodes is None:
103            self.number_of_full_nodes = self.number_of_nodes
104        else:
105            assert int(number_of_full_nodes)
106            self.number_of_full_nodes = number_of_full_nodes
107
108
109        if number_of_full_triangles is None:
110            self.number_of_full_triangles = self.number_of_triangles
111        else:
112            assert int(number_of_full_triangles)
113            self.number_of_full_triangles = number_of_full_triangles
114
115
116        #print self.number_of_full_nodes, self.number_of_nodes
117        #print self.number_of_full_triangles, self.number_of_triangles
118
119
120
121        # FIXME: this stores a geo_reference, but when coords are returned
122        # This geo_ref is not taken into account!
123        if geo_reference is None:
124            self.geo_reference = Geo_reference() #Use defaults
125        else:
126            self.geo_reference = geo_reference
127
128        # Input checks
129        msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
130        msg += 'The supplied array has the shape: %s'\
131               %str(self.triangles.shape)
132        assert len(self.triangles.shape) == 2, msg
133
134        msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
135        msg += 'The supplied array has the shape: %s'\
136               %str(self.nodes.shape)
137        assert len(self.nodes.shape) == 2, msg
138
139        msg = 'Vertex indices reference non-existing coordinate sets'
140        assert max(max(self.triangles)) <= self.nodes.shape[0], msg
141
142
143        # FIXME: Maybe move to statistics?
144        # Or use with get_extent
145        xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) ,
146                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
147
148        self.xy_extent = array(xy_extent, Float)
149
150
151        # Allocate space for geometric quantities
152        self.normals = zeros((N, 6), Float)
153        self.areas = zeros(N, Float)
154        self.edgelengths = zeros((N, 3), Float)
155
156        # Get x,y coordinates for all triangles and store
157        self.vertex_coordinates = V = self.compute_vertex_coordinates()
158
159
160        # Initialise each triangle
161        if verbose:
162            print 'General_mesh: Computing areas, normals and edgelenghts'
163
164        for i in range(N):
165            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
166
167            x0, y0 = V[3*i, :]
168            x1, y1 = V[3*i+1, :]
169            x2, y2 = V[3*i+2, :]
170
171            # Area
172            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
173
174            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
175            msg += ' is degenerate:  area == %f' %self.areas[i]
176            assert self.areas[i] > 0.0, msg
177
178
179            # Normals
180            # The normal vectors
181            #   - point outward from each edge
182            #   - are orthogonal to the edge
183            #   - have unit length
184            #   - Are enumerated according to the opposite corner:
185            #     (First normal is associated with the edge opposite
186            #     the first vertex, etc)
187            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
188
189            n0 = array([x2 - x1, y2 - y1])
190            l0 = sqrt(sum(n0**2))
191
192            n1 = array([x0 - x2, y0 - y2])
193            l1 = sqrt(sum(n1**2))
194
195            n2 = array([x1 - x0, y1 - y0])
196            l2 = sqrt(sum(n2**2))
197
198            # Normalise
199            n0 /= l0
200            n1 /= l1
201            n2 /= l2
202
203            # Compute and store
204            self.normals[i, :] = [n0[1], -n0[0],
205                                  n1[1], -n1[0],
206                                  n2[1], -n2[0]]
207
208            # Edgelengths
209            self.edgelengths[i, :] = [l0, l1, l2]
210
211
212        # Build vertex list
213        if verbose: print 'Building vertex list'
214        self.build_vertexlist()
215
216
217
218    def __len__(self):
219        return self.number_of_triangles
220
221
222    def __repr__(self):
223        return 'Mesh: %d vertices, %d triangles'\
224               %(self.nodes.shape[0], len(self))
225
226    def get_normals(self):
227        """Return all normal vectors.
228        Return normal vectors for all triangles as an Nx6 array
229        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
230        """
231        return self.normals
232
233
234    def get_normal(self, i, j):
235        """Return normal vector j of the i'th triangle.
236        Return value is the numeric array slice [x, y]
237        """
238        return self.normals[i, 2*j:2*j+2]
239
240
241    def get_nodes(self, absolute=False):
242        """Return all nodes in mesh.
243
244        The nodes are ordered in an Nx2 array where N is the number of nodes.
245        This is the same format they were provided in the constructor
246        i.e. without any duplication.
247
248        Boolean keyword argument absolute determines whether coordinates
249        are to be made absolute by taking georeference into account
250        Default is False as many parts of ANUGA expects relative coordinates.
251        (To see which, switch to default absolute=True and run tests).
252        """
253
254        N = self.number_of_full_nodes
255        V = self.nodes[:N,:]
256        if absolute is True:
257            if not self.geo_reference.is_absolute():
258                V = self.geo_reference.get_absolute(V)
259
260        return V
261
262
263
264    def get_vertex_coordinates(self, absolute=False):
265        """Return vertex coordinates for all triangles.
266
267        Return all vertex coordinates for all triangles as a 3*M x 2 array
268        where the jth vertex of the ith triangle is located in row 3*i+j and
269        M the number of triangles in the mesh.
270
271        Boolean keyword argument absolute determines whether coordinates
272        are to be made absolute by taking georeference into account
273        Default is False as many parts of ANUGA expects relative coordinates.
274        """
275
276        V = self.vertex_coordinates
277        if absolute is True:
278            if not self.geo_reference.is_absolute():
279                V = self.geo_reference.get_absolute(V)
280
281        return V
282
283
284
285    def get_vertex_coordinate(self, i, j, absolute=False):
286        """Return coordinates for vertex j of the i'th triangle.
287        Return value is the numeric array slice [x, y]
288        """
289
290        V = self.get_vertex_coordinates(absolute=absolute)
291        return V[3*i+j, :]
292
293
294    def compute_vertex_coordinates(self):
295        """Return all vertex coordinates for all triangles as a 3*M x 2 array
296        where the jth vertex of the ith triangle is located in row 3*i+j.
297
298        This function is used to precompute this important structure. Use
299        get_vertex coordinates to retrieve the points.
300        """
301
302        M = self.number_of_triangles
303        vertex_coordinates = zeros((3*M, 2), Float)
304
305        for i in range(M):
306            for j in range(3):
307                k = self.triangles[i,j] # Index of vertex j in triangle i
308                vertex_coordinates[3*i+j,:] = self.nodes[k]
309
310        return vertex_coordinates
311
312
313
314    def get_triangles(self, indices=None):
315        """Get mesh triangles.
316
317        Return Mx3 integer array where M is the number of triangles.
318        Each row corresponds to one triangle and the three entries are
319        indices into the mesh nodes which can be obtained using the method
320        get_nodes()
321
322        Optional argument, indices is the set of triangle ids of interest.
323        """
324
325        M = self.number_of_full_triangles
326
327        if indices is None:
328            indices = range(M)
329
330        return take(self.triangles, indices)
331
332
333
334    def get_disconnected_triangles(self):
335        """Get mesh based on nodes obtained from get_vertex_coordinates.
336
337        Return array Mx3 array of integers where each row corresponds to
338        a triangle. A triangle is a triplet of indices into
339        point coordinates obtained from get_vertex_coordinates and each
340        index appears only once
341
342        This provides a mesh where no triangles share nodes
343        (hence the name disconnected triangles) and different
344        nodes may have the same coordinates.
345
346        This version of the mesh is useful for storing meshes with
347        discontinuities at each node and is e.g. used for storing
348        data in sww files.
349
350        The triangles created will have the format
351
352        [[0,1,2],
353         [3,4,5],
354         [6,7,8],
355         ...
356         [3*M-3 3*M-2 3*M-1]]
357        """
358
359        M = len(self) # Number of triangles
360        K = 3*M       # Total number of unique vertices
361        T = reshape(array(range(K)).astype(Int), (M,3))
362
363        return T
364
365
366
367    def get_unique_vertices(self,  indices=None):
368        """FIXME(Ole): This function needs a docstring
369        """
370        triangles = self.get_triangles(indices=indices)
371        unique_verts = {}
372        for triangle in triangles:
373           unique_verts[triangle[0]] = 0
374           unique_verts[triangle[1]] = 0
375           unique_verts[triangle[2]] = 0
376        return unique_verts.keys()
377
378
379    def build_vertexlist(self):
380        """Build vertexlist indexed by vertex ids and for each entry (point id)
381        build a list of (triangles, vertex_id) pairs that use the point
382        as vertex.
383
384        The vertex list will have length N, where N is the number of nodes
385        in the mesh.
386
387        Preconditions:
388          self.nodes and self.triangles are defined
389
390        Postcondition:
391          self.vertexlist is built
392        """
393
394        vertexlist = [None]*self.number_of_nodes
395        for i in range(self.number_of_triangles):
396
397            a = self.triangles[i, 0]
398            b = self.triangles[i, 1]
399            c = self.triangles[i, 2]
400
401            #Register the vertices v as lists of
402            #(triangle_id, vertex_id) tuples associated with them
403            #This is used for averaging multiple vertex values.
404            for vertex_id, v in enumerate([a,b,c]):
405                if vertexlist[v] is None:
406                    vertexlist[v] = []
407
408                vertexlist[v].append( (i, vertex_id) )
409
410        self.vertexlist = vertexlist
411
412
413    def get_extent(self, absolute=False):
414        """Return min and max of all x and y coordinates
415
416        Boolean keyword argument absolute determines whether coordinates
417        are to be made absolute by taking georeference into account
418        """
419
420
421
422        C = self.get_vertex_coordinates(absolute=absolute)
423        X = C[:,0:6:2].copy()
424        Y = C[:,1:6:2].copy()
425
426        xmin = min(X.flat)
427        xmax = max(X.flat)
428        ymin = min(Y.flat)
429        ymax = max(Y.flat)
430
431        return xmin, xmax, ymin, ymax
432
433
434    def get_area(self):