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

Last change on this file since 7776 was 7519, checked in by steve, 14 years ago

Have been playing with using a slope limited velocity to calculate
fluxes (hence the addition of evolved_quantities as well as conserved
quantities.

But the commit is to fix a problem Rudy found with sww2dem. Seems
numpy.array2string is a little too clever, in that it summarizes
output if there is a long sequence of zeros to
[0.0, 0.0, 0.0, ... 0.0, 0.0 ,0.0] To get around this I have added
a call to numpy.set_options(threshold=sys.max_int) to turn this
behaviour off!

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