# source:anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py@6191

Last change on this file since 6191 was 6191, checked in by ole, 15 years ago

Refactored Mesh-Domain inheritance pattern into a composition pattern, thereby
paving the way for reuse of Mesh instance in fitting as per ticket:242.
Had to disable test for isinstance(domain, Domain) in quantity.py as it
failed for unknown reasons. All other tests and validation suite passes.

File size: 20.4 KB
Line
1import Numeric 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)
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, typecode=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           unique_verts[triangle] = 0
424           unique_verts[triangle] = 0
425           unique_verts[triangle] = 0
426        return unique_verts.keys()
427
428
429    def get_triangles_and_vertices_per_node(self, node=None):
430        """Get triangles associated with given node.
431
432        Return list of triangle_ids, vertex_ids for specified node.
433        If node in None or absent, this information will be returned
434        for all (full) nodes in a list L where L[v] is the triangle
435        list for node v.
436        """
437
438        triangle_list = []
439        if node is not None:
440            # Get index for this node
441            first = num.sum(self.number_of_triangles_per_node[:node])
442
443            # Get number of triangles for this node
444            count = self.number_of_triangles_per_node[node]
445
446            for i in range(count):
447                index = self.vertex_value_indices[first+i]
448
449                volume_id = index / 3
450                vertex_id = index % 3
451
452                triangle_list.append( (volume_id, vertex_id) )
453
454            triangle_list = num.array(triangle_list, num.Int)    #array default#
455        else:
456            # Get info for all nodes recursively.
457            # If need be, we can speed this up by
458            # working directly with the inverted triangle structure
459            for i in range(self.number_of_full_nodes):
460
461                L = self.get_triangles_and_vertices_per_node(node=i)
462
463                triangle_list.append(L)
464
465        return triangle_list
466
467
468    def build_inverted_triangle_structure(self):
469        """Build structure listing triangles belonging to each node
470
471        Two arrays are created and store as mesh attributes
472
473        number_of_triangles_per_node: An integer array of length N
474        listing for each node how many triangles use it. N is the number of
475        nodes in mesh.
476
477        vertex_value_indices: An array of length M listing indices into
478        triangles ordered by node number. The (triangle_id, vertex_id)
479        pairs are obtained from each index as (index/3, index%3) or each
480        index can be used directly into a flattened triangles array. This
481        is for example the case in the quantity.c where this structure is
482        used to average vertex values efficiently.
483
484
485        Example:
486
487        a = [0.0, 0.0] # node 0
488        b = [0.0, 2.0] # node 1
489        c = [2.0, 0.0] # node 2
490        d = [0.0, 4.0] # node 3
491        e = [2.0, 2.0] # node 4
492        f = [4.0, 0.0] # node 5
493
494        nodes = array([a, b, c, d, e, f])
495
496        #bac, bce, ecf, dbe, daf, dae
497        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
498
499
500        For this structure
501
502        number_of_triangles_per_node = [1 3 3 1 3 1]
503        which means that node a has 1 triangle associated with it, node b
504        has 3, node has 3 and so on.
505
506        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
507        which reflects the fact that
508        node 0 is used by triangle 0, vertex 1 (index = 1)
509        node 1 is used by triangle 0, vertex 0 (index = 0)
510                   and by triangle 1, vertex 0 (index = 3)
511                   and by triangle 3, vertex 1 (index = 10)
512        node 2 is used by triangle 0, vertex 2 (index = 2)
513                   and by triangle 1, vertex 1 (index = 4)
514                   and by triangle 2, vertex 1 (index = 7)
515        node 3 is used by triangle 3, vertex 0 (index = 9)
516        node 4 is used by triangle 1, vertex 2 (index = 5)
517                   and by triangle 2, vertex 0 (index = 6)
518                   and by triangle 3, vertex 2 (index = 11)
519        node 5 is used by triangle 2, vertex 2 (index = 8)
520
521
522        Preconditions:
523          self.nodes and self.triangles are defined
524
525        Postcondition:
526          self.number_of_triangles_per_node is built
527          self.vertex_value_indices is built
528        """
529
530        # Count number of triangles per node
531        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes)
532        for volume_id, triangle in enumerate(self.get_triangles()):
533            for vertex_id in triangle:
534                number_of_triangles_per_node[vertex_id] += 1
535
536        # Allocate space for inverted structure
537        number_of_entries = num.sum(number_of_triangles_per_node)
538        vertex_value_indices = num.zeros(number_of_entries)
539
540        # Register (triangle, vertex) indices for each node
541        vertexlist = [None]*self.number_of_full_nodes
542        for volume_id in range(self.number_of_full_triangles):
543
544            a = self.triangles[volume_id, 0]
545            b = self.triangles[volume_id, 1]
546            c = self.triangles[volume_id, 2]
547
548            for vertex_id, node_id in enumerate([a,b,c]):
549                if vertexlist[node_id] is None:
550                    vertexlist[node_id] = []
551
552                vertexlist[node_id].append( (volume_id, vertex_id) )
553
554
555        # Build inverted triangle index array
556        k = 0
557        for vertices in vertexlist:
558            if vertices is not None:
559                for volume_id, vertex_id in vertices:
560                    vertex_value_indices[k] = 3*volume_id + vertex_id
561
562                    k += 1
563
564        # Save structure
565        self.number_of_triangles_per_node = number_of_triangles_per_node
566        self.vertex_value_indices = vertex_value_indices
567
568
569
570
571    def get_extent(self, absolute=False):
572        """Return min and max of all x and y coordinates
573
574        Boolean keyword argument absolute determines whether coordinates
575        are to be made absolute by taking georeference into account
576        """
577
578
579
580        C = self.get_vertex_coordinates(absolute=absolute)
581        X = C[:,0:6:2].copy()
582        Y = C[:,1:6:2].copy()
583
584        xmin = min(X.flat)
585        xmax = max(X.flat)
586        ymin = min(Y.flat)
587        ymax = max(Y.flat)
588        #print "C",C
589        return xmin, xmax, ymin, ymax
590
591    def get_areas(self):
592        """Get areas of all individual triangles.
593        """
594        return self.areas
595
596    def get_area(self):