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

Last change on this file since 8279 was 8279, checked in by neweyv, 11 years ago

Additional information written to the log file. This information will assist in creating an algorithm to predict the runtime of ANUGA given various parameters.

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