# source:branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py@6428

Last change on this file since 6428 was 6428, checked in by rwilson, 15 years ago

numpy changes

File size: 19.5 KB
Line
1import numpy as num
2
3from anuga.coordinate_transforms.geo_reference import Geo_reference
4
5class General_mesh:
6    """Collection of 2D triangular elements
7
8    A triangular element is defined in terms of three vertex ids,
9    ordered counter clock-wise, each corresponding to a given node
10    which is represented as a coordinate set (x,y).
11    Vertices from different triangles can point to the same node.
12    The nodes are implemented as an Nx2 numeric array containing the
13    x and y coordinates.
14
15
16    To instantiate:
17       Mesh(nodes, triangles)
18
19    where
20
21      nodes is either a list of 2-tuples or an Nx2 numeric array of
22      floats representing all x, y coordinates in the mesh.
23
24      triangles is either a list of 3-tuples or an Mx3 numeric array of
25      integers representing indices of all vertices in the mesh.
26      Each vertex is identified by its index i in [0, N-1].
27
28
29    Example:
30
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        nodes = [a, b, c, e]
37        triangles = [ [1,0,2], [1,2,3] ]   # bac, bce
38
39        # Create mesh with two triangles: bac and bce
40        mesh = Mesh(nodes, triangles)
41
42
43
44    Other:
45
46      In addition mesh computes an Mx6 array called vertex_coordinates.
47      This structure is derived from coordinates and contains for each
48      triangle the three x,y coordinates at the vertices.
49
50      See neighbourmesh.py for a specialisation of the general mesh class
51      which includes information about neighbours and the mesh boundary.
52
53      The mesh object is purely geometrical and contains no information
54      about quantities defined on the mesh.
55
56    """
57
58    # FIXME: It would be a good idea to use geospatial data as an alternative
59    #        input
60    def __init__(self,
61                 nodes,
62                 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 in ANUGA domain'
91
92        self.triangles = num.array(triangles, num.int)
93        self.nodes = num.array(nodes, num.float)
94
95        # Register number of elements and nodes
96        self.number_of_triangles = N = self.triangles.shape
97        self.number_of_nodes = self.nodes.shape
98
99        if number_of_full_nodes is None:
100            self.number_of_full_nodes = self.number_of_nodes
101        else:
102            assert int(number_of_full_nodes)
103            self.number_of_full_nodes = number_of_full_nodes
104
105        if number_of_full_triangles is None:
106            self.number_of_full_triangles = self.number_of_triangles
107        else:
108            assert int(number_of_full_triangles)
109            self.number_of_full_triangles = number_of_full_triangles
110
111        # FIXME: this stores a geo_reference, but when coords are returned
112        # This geo_ref is not taken into account!
113        if geo_reference is None:
114            self.geo_reference = Geo_reference()    # Use defaults
115        else:
116            self.geo_reference = geo_reference
117
118        # Input checks
119        msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
120               'The supplied array has the shape: %s'
121               % str(self.triangles.shape))
122        assert len(self.triangles.shape) == 2, msg
123
124        msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples'
125               'The supplied array has the shape: %s' % str(self.nodes.shape))
126        assert len(self.nodes.shape) == 2, msg
127
128        msg = 'Vertex indices reference non-existing coordinate sets'
129        assert max(self.triangles.flat) < self.nodes.shape, msg
130
131        # FIXME: Maybe move to statistics?
132        # Or use with get_extent
133        xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]),
134                     max(self.nodes[:,0]), max(self.nodes[:,1])]
135
136        self.xy_extent = num.array(xy_extent, num.float)
137
138        # Allocate space for geometric quantities
139        self.normals = num.zeros((N, 6), num.float)
140        self.areas = num.zeros(N, num.float)
141        self.edgelengths = num.zeros((N, 3), num.float)
142
143        # Get x,y coordinates for all triangles and store
144        self.vertex_coordinates = V = self.compute_vertex_coordinates()
145
146        # Initialise each triangle
147        if verbose:
148            print 'General_mesh: Computing areas, normals and edgelengths'
149
150        for i in range(N):
151            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N)
152
153            x0, y0 = V[3*i, :]
154            x1, y1 = V[3*i+1, :]
155            x2, y2 = V[3*i+2, :]
156
157            # Area
158            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
159
160            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2)
161            msg += ' is degenerate:  area == %f' % self.areas[i]
162            assert self.areas[i] > 0.0, msg
163
164            # Normals
165            # The normal vectors
166            #   - point outward from each edge
167            #   - are orthogonal to the edge
168            #   - have unit length
169            #   - Are enumerated according to the opposite corner:
170            #     (First normal is associated with the edge opposite
171            #     the first vertex, etc)
172            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
173            n0 = num.array([x2-x1, y2-y1], num.float)
174            l0 = num.sqrt(num.sum(n0**2))
175
176            n1 = num.array([x0-x2, y0-y2], num.float)
177            l1 = num.sqrt(num.sum(n1**2))
178
179            n2 = num.array([x1-x0, y1-y0], num.float)
180            l2 = num.sqrt(num.sum(n2**2))
181
182            # Normalise
183            n0 /= l0
184            n1 /= l1
185            n2 /= l2
186
187            # Compute and store
188            self.normals[i, :] = [n0, -n0,
189                                  n1, -n1,
190                                  n2, -n2]
191
192            # Edgelengths
193            self.edgelengths[i, :] = [l0, l1, l2]
194
195        # Build structure listing which trianglse belong to which node.
196        if verbose: print 'Building inverted triangle structure'
197        self.build_inverted_triangle_structure()
198
199    def __len__(self):
200        return self.number_of_triangles
201
202    def __repr__(self):
203        return ('Mesh: %d vertices, %d triangles'
204                % (self.nodes.shape, len(self)))
205
206    def get_normals(self):
207        """Return all normal vectors.
208
209        Return normal vectors for all triangles as an Nx6 array
210        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
211        """
212
213        return self.normals
214
215    def get_normal(self, i, j):
216        """Return normal vector j of the i'th triangle.
217
218        Return value is the numeric array slice [x, y]
219        """
220
221        return self.normals[i, 2*j:2*j+2]
222
223    def get_number_of_nodes(self):
224        return self.number_of_nodes
225
226    def get_nodes(self, absolute=False):
227        """Return all nodes in mesh.
228
229        The nodes are ordered in an Nx2 array where N is the number of nodes.
230        This is the same format they were provided in the constructor
231        i.e. without any duplication.
232
233        Boolean keyword argument absolute determines whether coordinates
234        are to be made absolute by taking georeference into account
235        Default is False as many parts of ANUGA expects relative coordinates.
236        (To see which, switch to default absolute=True and run tests).
237        """
238
239        N = self.number_of_full_nodes
240        V = self.nodes[:N,:]
241        if absolute is True:
242            if not self.geo_reference.is_absolute():
243                V = self.geo_reference.get_absolute(V)
244
245        return V
246
247    def get_node(self, i, absolute=False):
248        """Return node coordinates for triangle i.
249
250        Boolean keyword argument absolute determines whether coordinates
251        are to be made absolute by taking georeference into account
252        Default is False as many parts of ANUGA expects relative coordinates.
253        (To see which, switch to default absolute=True and run tests).
254        """
255
256        V = self.nodes[i,:]
257        if absolute is True:
258            if not self.geo_reference.is_absolute():
259                V += num.array([self.geo_reference.get_xllcorner(),
260                                self.geo_reference.get_yllcorner()], num.float)
261        return V
262
263    def get_vertex_coordinates(self, triangle_id=None, absolute=False):
264        """Return vertex coordinates for all triangles.
265
266        Return all vertex coordinates for all triangles as a 3*M x 2 array
267        where the jth vertex of the ith triangle is located in row 3*i+j and
268        M the number of triangles in the mesh.
269
270        if triangle_id is specified (an integer) the 3 vertex coordinates
271        for triangle_id are returned.
272
273        Boolean keyword argument absolute determines whether coordinates
274        are to be made absolute by taking georeference into account
275        Default is False as many parts of ANUGA expects relative coordinates.
276        """
277
278        V = self.vertex_coordinates
279
280        if triangle_id is None:
281            if absolute is True:
282                if not self.geo_reference.is_absolute():
283                    V = self.geo_reference.get_absolute(V)
284            return V
285        else:
286            i = triangle_id
287            msg = 'triangle_id must be an integer'
288            assert int(i) == i, msg
289            assert 0 <= i < self.number_of_triangles
290
291            i3 = 3*i
292            if absolute is True and not self.geo_reference.is_absolute():
293                offset=num.array([self.geo_reference.get_xllcorner(),
294                                  self.geo_reference.get_yllcorner()], num.float)
295                return num.array([V[i3,:]+offset,
296                                  V[i3+1,:]+offset,
297                                  V[i3+2,:]+offset], num.float)
298            else:
299                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
300
301    def get_vertex_coordinate(self, i, j, absolute=False):
302        """Return coordinates for vertex j of the i'th triangle.
303        Return value is the numeric array slice [x, y]
304        """
305
306        msg = 'vertex id j must be an integer in [0,1,2]'
307        assert j in [0,1,2], msg
308
309        V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute)
310        return V[j,:]
311
312    def compute_vertex_coordinates(self):
313        """Return all vertex coordinates for all triangles as a 3*M x 2 array
314        where the jth vertex of the ith triangle is located in row 3*i+j.
315
316        This function is used to precompute this important structure. Use
317        get_vertex coordinates to retrieve the points.
318        """
319
320        M = self.number_of_triangles
321        vertex_coordinates = num.zeros((3*M, 2), num.float)
322
323        for i in range(M):
324            for j in range(3):
325                k = self.triangles[i,j] # Index of vertex j in triangle i
326                vertex_coordinates[3*i+j,:] = self.nodes[k]
327
328        return vertex_coordinates
329
330    def get_triangles(self, indices=None):
331        """Get mesh triangles.
332
333        Return Mx3 integer array where M is the number of triangles.
334        Each row corresponds to one triangle and the three entries are
335        indices into the mesh nodes which can be obtained using the method
336        get_nodes()
337
338        Optional argument, indices is the set of triangle ids of interest.
339        """
340
341        M = self.number_of_full_triangles
342
343        if indices is None:
344            return self.triangles
345
346        return num.take(self.triangles, indices, axis=0)
347
348    def get_disconnected_triangles(self):
349        """Get mesh based on nodes obtained from get_vertex_coordinates.
350
351        Return array Mx3 array of integers where each row corresponds to
352        a triangle. A triangle is a triplet of indices into
353        point coordinates obtained from get_vertex_coordinates and each
354        index appears only once
355
356        This provides a mesh where no triangles share nodes
357        (hence the name disconnected triangles) and different
358        nodes may have the same coordinates.
359
360        This version of the mesh is useful for storing meshes with
361        discontinuities at each node and is e.g. used for storing
362        data in sww files.
363
364        The triangles created will have the format
365        [[0,1,2],
366         [3,4,5],
367         [6,7,8],
368         ...
369         [3*M-3 3*M-2 3*M-1]]
370        """
371
372        M = len(self) # Number of triangles
373        K = 3*M       # Total number of unique vertices
374        return num.reshape(num.arange(K, dtype=num.int), (M,3))
375
376    def get_unique_vertices(self, indices=None):
377        """FIXME(Ole): This function needs a docstring"""
378
379        triangles = self.get_triangles(indices=indices)
380        unique_verts = {}
381        for triangle in triangles:
382            unique_verts[triangle] = 0
383            unique_verts[triangle] = 0
384            unique_verts[triangle] = 0
385        return unique_verts.keys()
386
387    def get_triangles_and_vertices_per_node(self, node=None):
388        """Get triangles associated with given node.
389
390        Return list of triangle_ids, vertex_ids for specified node.
391        If node in None or absent, this information will be returned
392        for all (full) nodes in a list L where L[v] is the triangle
393        list for node v.
394        """
395
396        triangle_list = []
397        if node is not None:
398            # Get index for this node
399            first = num.sum(self.number_of_triangles_per_node[:node])
400
401            # Get number of triangles for this node
402            count = self.number_of_triangles_per_node[node]
403
404            for i in range(count):
405                index = self.vertex_value_indices[first+i]
406
407                volume_id = index / 3
408                vertex_id = index % 3
409
410                triangle_list.append( (volume_id, vertex_id) )
411
412            triangle_list = num.array(triangle_list, num.int)    #array default#
413        else:
414            # Get info for all nodes recursively.
415            # If need be, we can speed this up by
416            # working directly with the inverted triangle structure
417            for i in range(self.number_of_full_nodes):
418                L = self.get_triangles_and_vertices_per_node(node=i)
419                triangle_list.append(L)
420
421        return triangle_list
422
423    def build_inverted_triangle_structure(self):
424        """Build structure listing triangles belonging to each node
425
426        Two arrays are created and store as mesh attributes
427
428        number_of_triangles_per_node: An integer array of length N
429        listing for each node how many triangles use it. N is the number of
430        nodes in mesh.
431
432        vertex_value_indices: An array of length M listing indices into
433        triangles ordered by node number. The (triangle_id, vertex_id)
434        pairs are obtained from each index as (index/3, index%3) or each
435        index can be used directly into a flattened triangles array. This
436        is for example the case in the quantity.c where this structure is
437        used to average vertex values efficiently.
438
439        Example:
440        a = [0.0, 0.0] # node 0
441        b = [0.0, 2.0] # node 1
442        c = [2.0, 0.0] # node 2
443        d = [0.0, 4.0] # node 3
444        e = [2.0, 2.0] # node 4
445        f = [4.0, 0.0] # node 5
446        nodes = array([a, b, c, d, e, f])
447
448        #                    bac,     bce,     ecf,     dbe
449        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
450
451        For this structure:
452        number_of_triangles_per_node = [1 3 3 1 3 1]
453        which means that node a has 1 triangle associated with it, node b
454        has 3, node has 3 and so on.
455
456        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
457        which reflects the fact that
458        node 0 is used by triangle 0, vertex 1 (index = 1)
459        node 1 is used by triangle 0, vertex 0 (index = 0)
460                   and by triangle 1, vertex 0 (index = 3)
461                   and by triangle 3, vertex 1 (index = 10)
462        node 2 is used by triangle 0, vertex 2 (index = 2)
463                   and by triangle 1, vertex 1 (index = 4)
464                   and by triangle 2, vertex 1 (index = 7)
465        node 3 is used by triangle 3, vertex 0 (index = 9)
466        node 4 is used by triangle 1, vertex 2 (index = 5)
467                   and by triangle 2, vertex 0 (index = 6)
468                   and by triangle 3, vertex 2 (index = 11)
469        node 5 is used by triangle 2, vertex 2 (index = 8)
470
471        Preconditions:
472          self.nodes and self.triangles are defined
473
474        Postcondition:
475          self.number_of_triangles_per_node is built
476          self.vertex_value_indices is built
477        """
478
479        # Count number of triangles per node
480        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes,
481                                                 num.int)       #array default#
482        for volume_id, triangle in enumerate(self.get_triangles()):
483            for vertex_id in triangle:
484                number_of_triangles_per_node[vertex_id] += 1
485
486        # Allocate space for inverted structure
487        number_of_entries = num.sum(number_of_triangles_per_node)
488        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
489
490        # Register (triangle, vertex) indices for each node
491        vertexlist = [None] * self.number_of_full_nodes
492        for volume_id in range(self.number_of_full_triangles):
493            a = self.triangles[volume_id, 0]
494            b = self.triangles[volume_id, 1]
495            c = self.triangles[volume_id, 2]
496
497            for vertex_id, node_id in enumerate([a, b, c]):
498                if vertexlist[node_id] is None:
499                    vertexlist[node_id] = []
500                vertexlist[node_id].append((volume_id, vertex_id))
501
502        # Build inverted triangle index array
503        k = 0
504        for vertices in vertexlist:
505            if vertices is not None:
506                for volume_id, vertex_id in vertices:
507                    vertex_value_indices[k] = 3*volume_id + vertex_id
508                    k += 1
509
510        # Save structure
511        self.number_of_triangles_per_node = number_of_triangles_per_node
512        self.vertex_value_indices = vertex_value_indices
513
514    def get_extent(self, absolute=False):
515        """Return min and max of all x and y coordinates
516
517        Boolean keyword argument absolute determines whether coordinates
518        are to be made absolute by taking georeference into account
519        """
520
521        C = self.get_vertex_coordinates(absolute=absolute)
522        X = C[:,0:6:2].copy()
523        Y = C[:,1:6:2].copy()
524
525        xmin = min(X.flat)
526        xmax = max(X.flat)
527        ymin = min(Y.flat)
528        ymax = max(Y.flat)
529
530        return xmin, xmax, ymin, ymax
531
532    def get_areas(self):
533        """Get areas of all individual triangles."""
534
535        return self.areas
536
537    def get_area(self):