# source:anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/general_mesh.py@5899

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

Initial NumPy? changes (again!).

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