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

Last change on this file since 6410 was 6410, checked in by rwilson, 14 years ago

numpy changes.

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