# source:trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py@8200

Last change on this file since 8200 was 8200, checked in by steve, 11 years ago

Pulled the number_of_full_nodes and triangle out of generic_mesh up into generic_domain. These "numbers" are now calculated from ghost_recv and full_send data

File size: 22.6 KB
Line
1import copy
2import numpy as num
3
4from anuga.coordinate_transforms.geo_reference import Geo_reference
5import anuga.utilities.log as log
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,
63                 nodes,
64                 triangles,
65                 geo_reference=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        """
83
84        if verbose: log.critical('General_mesh: Building basic mesh structure '
85                                 'in ANUGA domain')
86
87        self.triangles = num.array(triangles, num.int)
88        self.nodes = num.array(nodes, num.float)
89
90        # Register number of elements and nodes
91        self.number_of_triangles = N = self.triangles.shape[0]
92        self.number_of_nodes = self.nodes.shape[0]
93
94
95        # FIXME: this stores a geo_reference, but when coords are returned
96        # This geo_ref is not taken into account!
97        if geo_reference is None:
98            self.geo_reference = Geo_reference()    # Use defaults
99        else:
100            self.geo_reference = geo_reference
101
102        # Input checks
103        msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
104               'The supplied array has the shape: %s'
105               % str(self.triangles.shape))
106        assert len(self.triangles.shape) == 2, msg
107
108        msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples'
109               'The supplied array has the shape: %s' % str(self.nodes.shape))
110        assert len(self.nodes.shape) == 2, msg
111
112        msg = 'Vertex indices reference non-existing coordinate sets'
113        assert num.max(self.triangles) < self.nodes.shape[0], msg
114
115        # FIXME: Maybe move to statistics?
116        # Or use with get_extent
117        xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]),
118                     max(self.nodes[:,0]), max(self.nodes[:,1])]
119
120        self.xy_extent = num.array(xy_extent, num.float)
121
122        # Allocate space for geometric quantities
123        self.normals = num.zeros((N, 6), num.float)
124        self.areas = num.zeros(N, num.float)
125        self.edgelengths = num.zeros((N, 3), num.float)
126
127        # Get x,y coordinates for all triangle vertices and store
128        self.vertex_coordinates = V = self.compute_vertex_coordinates()
129
130        # Get x,y coordinates for all triangle edge midpoints and store
131        self.edge_midpoint_coordinates  = self.compute_edge_midpoint_coordinates()
132
133        # Initialise each triangle
134        if verbose:
135            log.critical('General_mesh: Computing areas, normals '
136                         'and edgelengths')
137
138        for i in range(N):
139            if verbose and i % ((N+10)/10) == 0: log.critical('(%d/%d)' % (i, N))
140
141            x0, y0 = V[3*i, :]
142            x1, y1 = V[3*i+1, :]
143            x2, y2 = V[3*i+2, :]
144
145
146            i0 = self.triangles[i][0]
147            i1 = self.triangles[i][1]
148            i2 = self.triangles[i][2]
149
150            assert x0 == self.nodes[i0][0]
151            assert y0 == self.nodes[i0][1]
152
153            assert x1 == self.nodes[i1][0]
154            assert y1 == self.nodes[i1][1]
155
156            assert x2 == self.nodes[i2][0]
157            assert y2 == self.nodes[i2][1]
158
159            # Area
160            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
161
162            msg = 'Triangle %g (%f,%f), (%f,%f), (%f, %f)' % (i,x0,y0,x1,y1,x2,y2)
163            msg += ' is degenerate:  area == %f' % self.areas[i]
164            assert self.areas[i] > 0.0, msg
165
166            # Normals
167            # The normal vectors
168            #   - point outward from each edge
169            #   - are orthogonal to the edge
170            #   - have unit length
171            #   - Are enumerated according to the opposite corner:
172            #     (First normal is associated with the edge opposite
173            #     the first vertex, etc)
174            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
175            n0 = num.array([x2-x1, y2-y1], num.float)
176            l0 = num.sqrt(num.sum(n0**2))
177
178            n1 = num.array([x0-x2, y0-y2], num.float)
179            l1 = num.sqrt(num.sum(n1**2))
180
181            n2 = num.array([x1-x0, y1-y0], num.float)
182            l2 = num.sqrt(num.sum(n2**2))
183
184            # Normalise
185            n0 /= l0
186            n1 /= l1
187            n2 /= l2
188
189            # Compute and store
190            self.normals[i, :] = [n0[1], -n0[0],
191                                  n1[1], -n1[0],
192                                  n2[1], -n2[0]]
193
194            # Edgelengths
195            self.edgelengths[i, :] = [l0, l1, l2]
196
197        # Build structure listing which triangles belong to which node.
198        if verbose: log.critical('Building inverted triangle structure')
199        self.build_inverted_triangle_structure()
200
201    def __len__(self):
202        return self.number_of_triangles
203
204    def __repr__(self):
205        return ('Mesh: %d vertices, %d triangles'
206                % (self.nodes.shape[0], len(self)))
207
208    def get_normals(self):
209        """Return all normal vectors.
210
211        Return normal vectors for all triangles as an Nx6 array
212        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
213        """
214
215        return self.normals
216
217    def get_normal(self, i, j):
218        """Return normal vector j of the i'th triangle.
219
220        Return value is the numeric array slice [x, y]
221        """
222
223        return self.normals[i, 2*j:2*j+2]
224
225    def get_edgelength(self, i, j):
226        """Return length of j'th edge of the i'th triangle.
227        Return value is the numeric array slice [x, y]
228        """
229        return self.edgelengths[i, j]
230
231
232    def get_number_of_triangles(self):
233        return self.number_of_triangles
234
235
236    def get_number_of_nodes(self):
237        return self.number_of_nodes
238
239
240    def get_nodes(self, absolute=False):
241        """Return all nodes in mesh.
242
243        The nodes are ordered in an Nx2 array where N is the number of nodes.
244        This is the same format they were provided in the constructor
245        i.e. without any duplication.
246
247        Boolean keyword argument absolute determines whether coordinates
248        are to be made absolute by taking georeference into account
249        Default is False as many parts of ANUGA expects relative coordinates.
250        (To see which, switch to default absolute=True and run tests).
251        """
252
253        N = self.number_of_nodes
254        V = self.nodes[:N,:]
255        if absolute is True:
256            if not self.geo_reference.is_absolute():
257                V = self.geo_reference.get_absolute(V)
258
259        return V
260
261    def get_node(self, i, absolute=False):
262        """Return node coordinates for triangle i.
263
264        Boolean keyword argument absolute determines whether coordinates
265        are to be made absolute by taking georeference into account
266        Default is False as many parts of ANUGA expects relative coordinates.
267        (To see which, switch to default absolute=True and run tests).
268
269        Note: This method returns a modified _copy_ of the nodes slice if
270              absolute is True.  If absolute is False, just return the slice.
271              This is related to the ensure_numeric() returning a copy problem.
272        """
273
274        V = self.nodes[i,:]
275        if absolute is True:
276            if not self.geo_reference.is_absolute():
277                # get a copy so as not to modify the internal self.nodes array
278                V = copy.copy(V)
279                V += num.array([self.geo_reference.get_xllcorner(),
280                                self.geo_reference.get_yllcorner()], num.float)
281        return V
282
283    def get_vertex_coordinates(self, triangle_id=None, absolute=False):
284        """Return vertex coordinates for all triangles.
285
286        Return all vertex coordinates for all triangles as a 3*M x 2 array
287        where the jth vertex of the ith triangle is located in row 3*i+j and
288        M the number of triangles in the mesh.
289
290        if triangle_id is specified (an integer) the 3 vertex coordinates
291        for triangle_id are returned.
292
293        Boolean keyword argument absolute determines whether coordinates
294        are to be made absolute by taking georeference into account
295        Default is False as many parts of ANUGA expects relative coordinates.
296        """
297
298        V = self.vertex_coordinates
299
300        if triangle_id is None:
301            if absolute is True:
302                if not self.geo_reference.is_absolute():
303                    V = self.geo_reference.get_absolute(V)
304            return V
305        else:
306            i = triangle_id
307            msg = 'triangle_id must be an integer'
308            assert int(i) == i, msg
309            assert 0 <= i < self.number_of_triangles
310
311            i3 = 3*i
312            if absolute is True and not self.geo_reference.is_absolute():
313                offset=num.array([self.geo_reference.get_xllcorner(),
314                                  self.geo_reference.get_yllcorner()], num.float)
315
316                return V[i3:i3+3,:] + offset
317            else:
318                return V[i3:i3+3,:]
319
320    def get_vertex_coordinate(self, i, j, absolute=False):
321        """Return coordinates for vertex j of the i'th triangle.
322        Return value is the numeric array slice [x, y]
323        """
324
325        msg = 'vertex id j must be an integer in [0,1,2]'
326        assert j in [0,1,2], msg
327
328        V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute)
329        return V[j,:]
330
331
332
333    def compute_vertex_coordinates(self):
334        """Return all vertex coordinates for all triangles as a 3*M x 2 array
335        where the jth vertex of the ith triangle is located in row 3*i+j.
336
337        This function is used to precompute this important structure. Use
338        get_vertex coordinates to retrieve the points.
339        """
340
341        M = self.number_of_triangles
342        vertex_coordinates = num.zeros((3*M, 2), num.float)
343
344        for i in range(M):
345            for j in range(3):
346                k = self.triangles[i,j] # Index of vertex j in triangle i
347                vertex_coordinates[3*i+j,:] = self.nodes[k]
348
349        return vertex_coordinates
350
351
352    def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
353        """Return edge midpoint coordinates for all triangles or from particular triangle.
354
355        Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
356        where the jth midpoint of the ith triangle is located in row 3*i+j and
357        M the number of triangles in the mesh.
358
359        if triangle_id is specified (an integer) the 3 midpoint coordinates
360        for triangle_id are returned.
361
362        Boolean keyword argument absolute determines whether coordinates
363        are to be made absolute by taking georeference into account
364        Default is False as many parts of ANUGA expects relative coordinates.
365        """
366
367        E = self.edge_midpoint_coordinates
368
369        if triangle_id is None:
370            if absolute is True:
371                if not self.geo_reference.is_absolute():
372                    E = self.geo_reference.get_absolute(E)
373            return E
374        else:
375            i = triangle_id
376            msg = 'triangle_id must be an integer'
377            assert int(i) == i, msg
378            assert 0 <= i < self.number_of_triangles
379
380            i3 = 3*i
381            if absolute is True and not self.geo_reference.is_absolute():
382                offset=num.array([self.geo_reference.get_xllcorner(),
383                                  self.geo_reference.get_yllcorner()], num.float)
384
385                return E[i3:i3+3,:] + offset
386            else:
387                return E[i3:i3+3,:]
388
389
390    def get_edge_midpoint_coordinate(self, i, j, absolute=False):
391        """Return coordinates for edge midpoint j of the i'th triangle.
392        Return value is the numeric array slice [x, y]
393        """
394
395        msg = 'edge midpoint id j must be an integer in [0,1,2]'
396        assert j in [0,1,2], msg
397
398        E = self.get_edge_midpoint_coordinates(triangle_id=i, absolute=absolute)
399        return E[j,:] # Return (x, y) for edge mid point
400
401
402    def compute_edge_midpoint_coordinates(self):
403        """Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
404        where the jth edge midpoint of the ith triangle is located in row 3*i+j.
405
406        This function is used to precompute this important structure. Use
407        get_edge_midpoint_coordinates to retrieve the points.
408
409        Assumes that vertex_coordinates have been computed
410        """
411
412        M = self.number_of_triangles
413        E = num.zeros((3*M, 2), num.float)
414
415        V = self.vertex_coordinates
416
417        V0 = V[0:3*M:3, :]
418        V1 = V[1:3*M:3, :]
419        V2 = V[2:3*M:3, :]
420
421
422        #print V.shape, V0.shape, V1.shape, V2.shape
423
424        #print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
425        E[0:3*M:3, :] = 0.5*(V1+V2)
426        E[1:3*M:3, :] = 0.5*(V2+V0)
427        E[2:3*M:3, :] = 0.5*(V0+V1)
428
429        return E
430
431
432
433    def get_triangles(self, indices=None):
434        """Get mesh triangles.
435
436        Return Mx3 integer array where M is the number of triangles.
437        Each row corresponds to one triangle and the three entries are
438        indices into the mesh nodes which can be obtained using the method
439        get_nodes()
440
441        Optional argument, indices is the set of triangle ids of interest.
442        """
443
444
445        if indices is None:
446            return self.triangles
447
448        return num.take(self.triangles, indices, axis=0)
449
450    def get_disconnected_triangles(self):
451        """Get mesh based on nodes obtained from get_vertex_coordinates.
452
453        Return array Mx3 array of integers where each row corresponds to
454        a triangle. A triangle is a triplet of indices into
455        point coordinates obtained from get_vertex_coordinates and each
456        index appears only once
457
458        This provides a mesh where no triangles share nodes
459        (hence the name disconnected triangles) and different
460        nodes may have the same coordinates.
461
462        This version of the mesh is useful for storing meshes with
463        discontinuities at each node and is e.g. used for storing
464        data in sww files.
465
466        The triangles created will have the format
467        [[0,1,2],
468         [3,4,5],
469         [6,7,8],
470         ...
471         [3*M-3 3*M-2 3*M-1]]
472        """
473
474        M = len(self) # Number of triangles
475        K = 3*M       # Total number of unique vertices
476        return num.reshape(num.arange(K, dtype=num.int), (M,3))
477
478    def get_unique_vertices(self, indices=None):
479        """FIXME(Ole): This function needs a docstring"""
480
481        triangles = self.get_triangles(indices=indices)
482        unique_verts = {}
483        for triangle in triangles:
484            unique_verts[triangle[0]] = 0
485            unique_verts[triangle[1]] = 0
486            unique_verts[triangle[2]] = 0
487        return unique_verts.keys()
488
489    def get_triangles_and_vertices_per_node(self, node=None):
490        """Get triangles associated with given node.
491
492        Return list of triangle_ids, vertex_ids for specified node.
493        If node in None or absent, this information will be returned
494        for all nodes in a list L where L[v] is the triangle
495        list for node v.
496        """
497
498        triangle_list = []
499        if node is not None:
500            # Get index for this node
501            first = num.sum(self.number_of_triangles_per_node[:node])
502
503            # Get number of triangles for this node
504            count = self.number_of_triangles_per_node[node]
505
506            for i in range(count):
507                index = self.vertex_value_indices[first+i]
508
509                volume_id = index / 3
510                vertex_id = index % 3
511
512                triangle_list.append( (volume_id, vertex_id) )
513
514            triangle_list = num.array(triangle_list, num.int)    #array default#
515        else:
516            # Get info for all nodes recursively.
517            # If need be, we can speed this up by
518            # working directly with the inverted triangle structure
519            for i in range(self.number_of_nodes):
520                L = self.get_triangles_and_vertices_per_node(node=i)
521                triangle_list.append(L)
522
523        return triangle_list
524
525    def build_inverted_triangle_structure(self):
526        """Build structure listing triangles belonging to each node
527
528        Two arrays are created and store as mesh attributes
529
530        number_of_triangles_per_node: An integer array of length N
531        listing for each node how many triangles use it. N is the number of
532        nodes in mesh.
533
534        vertex_value_indices: An array of length M listing indices into
535        triangles ordered by node number. The (triangle_id, vertex_id)
536        pairs are obtained from each index as (index/3, index%3) or each
537        index can be used directly into a flattened triangles array. This
538        is for example the case in the quantity.c where this structure is
539        used to average vertex values efficiently.
540
541        Example:
542        a = [0.0, 0.0] # node 0
543        b = [0.0, 2.0] # node 1
544        c = [2.0, 0.0] # node 2
545        d = [0.0, 4.0] # node 3
546        e = [2.0, 2.0] # node 4
547        f = [4.0, 0.0] # node 5
548        nodes = array([a, b, c, d, e, f])
549
550        #                    bac,     bce,     ecf,     dbe
551        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
552
553        For this structure:
554        number_of_triangles_per_node = [1 3 3 1 3 1]
555        which means that node a has 1 triangle associated with it, node b
556        has 3, node has 3 and so on.
557
558        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
559        which reflects the fact that
560        node 0 is used by triangle 0, vertex 1 (index = 1)
561        node 1 is used by triangle 0, vertex 0 (index = 0)
562                   and by triangle 1, vertex 0 (index = 3)
563                   and by triangle 3, vertex 1 (index = 10)
564        node 2 is used by triangle 0, vertex 2 (index = 2)
565                   and by triangle 1, vertex 1 (index = 4)
566                   and by triangle 2, vertex 1 (index = 7)
567        node 3 is used by triangle 3, vertex 0 (index = 9)
568        node 4 is used by triangle 1, vertex 2 (index = 5)
569                   and by triangle 2, vertex 0 (index = 6)
570                   and by triangle 3, vertex 2 (index = 11)
571        node 5 is used by triangle 2, vertex 2 (index = 8)
572
573        Preconditions:
574          self.nodes and self.triangles are defined
575
576        Postcondition:
577          self.number_of_triangles_per_node is built
578          self.vertex_value_indices is built
579        """
580
581        # Count number of triangles per node
582        number_of_triangles_per_node = num.zeros(self.number_of_nodes,
583                                                 num.int)       #array default#
584        for volume_id, triangle in enumerate(self.get_triangles()):
585            for vertex_id in triangle:
586                number_of_triangles_per_node[vertex_id] += 1
587
588        # Allocate space for inverted structure
589        number_of_entries = num.sum(number_of_triangles_per_node)
590        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
591
592        # Register (triangle, vertex) indices for each node
593        vertexlist = [None] * self.number_of_nodes
594        for volume_id in range(self.number_of_triangles):
595            a = self.triangles[volume_id, 0]
596            b = self.triangles[volume_id, 1]
597            c = self.triangles[volume_id, 2]
598
599            for vertex_id, node_id in enumerate([a, b, c]):
600                if vertexlist[node_id] is None:
601                    vertexlist[node_id] = []
602                vertexlist[node_id].append((volume_id, vertex_id))
603
604        # Build inverted triangle index array
605        k = 0
606        for vertices in vertexlist:
607            if vertices is not None:
608                for volume_id, vertex_id in vertices:
609                    vertex_value_indices[k] = 3*volume_id + vertex_id
610                    k += 1
611
612        # Save structure
613        self.number_of_triangles_per_node = number_of_triangles_per_node
614        self.vertex_value_indices = vertex_value_indices
615
616    def get_extent(self, absolute=False):
617        """Return min and max of all x and y coordinates
618
619        Boolean keyword argument absolute determines whether coordinates
620        are to be made absolute by taking georeference into account
621        """
622
623        C = self.get_vertex_coordinates(absolute=absolute)
624        X = C[:,0:6:2].copy()
625        Y = C[:,1:6:2].copy()
626
627        xmin = num.min(X)
628        xmax = num.max(X)
629        ymin = num.min(Y)
630        ymax = num.max(Y)
631
632        return xmin, xmax, ymin, ymax
633
634    def get_areas(self):
635        """Get areas of all individual triangles."""
636
637        return self.areas
638
639    def get_area(self):