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

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

Work on ticket:175 with special consideration to mesh intersections
Also reinstated disabled Okada test which currently fails but is faster.

File size: 20.0 KB
Line
1
2from Numeric import concatenate, reshape, take, allclose
3from Numeric import array, zeros, Int, Float, sqrt, sum, arange
4
5from anuga.coordinate_transforms.geo_reference import Geo_reference
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, nodes, triangles,
63                 geo_reference=None,
64                 number_of_full_nodes=None,
65                 number_of_full_triangles=None,
66                 verbose=False):
67        """Build triangular 2d mesh from nodes and triangle information
68
69        Input:
70
71          nodes: x,y coordinates represented as a sequence of 2-tuples or
72                 a Nx2 Numeric array of floats.
73
74          triangles: sequence of 3-tuples or Mx3 Numeric array of
75                     non-negative integers representing indices into
76                     the nodes array.
77
78          georeference (optional): If specified coordinates are
79          assumed to be relative to this origin.
80
81
82        number_of_full_nodes and number_of_full_triangles relate to
83        parallelism when each mesh has an extra layer of ghost points and
84        ghost triangles attached to the end of the two arrays.
85        In this case it is usefull to specify the number of real (called full)
86        nodes and triangles. If omitted they will default to all.
87
88        """
89
90        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
91
92        self.triangles = array(triangles, Int)
93        self.nodes = array(nodes, Float)
94
95
96        # Register number of elements and nodes
97        self.number_of_triangles = N = self.triangles.shape
98        self.number_of_nodes = self.nodes.shape
99
100
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
109        if number_of_full_triangles is None:
110            self.number_of_full_triangles = self.number_of_triangles
111        else:
112            assert int(number_of_full_triangles)
113            self.number_of_full_triangles = number_of_full_triangles
114
115
116        #print self.number_of_full_nodes, self.number_of_nodes
117        #print self.number_of_full_triangles, self.number_of_triangles
118
119
120
121        # FIXME: this stores a geo_reference, but when coords are returned
122        # This geo_ref is not taken into account!
123        if geo_reference is None:
124            self.geo_reference = Geo_reference() #Use defaults
125        else:
126            self.geo_reference = geo_reference
127
128        # Input checks
129        msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
130        msg += 'The supplied array has the shape: %s'\
131               %str(self.triangles.shape)
132        assert len(self.triangles.shape) == 2, msg
133
134        msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
135        msg += 'The supplied array has the shape: %s'\
136               %str(self.nodes.shape)
137        assert len(self.nodes.shape) == 2, msg
138
139        msg = 'Vertex indices reference non-existing coordinate sets'
140        assert max(max(self.triangles)) <= self.nodes.shape, msg
141
142
143        # FIXME: Maybe move to statistics?
144        # Or use with get_extent
145        xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) ,
146                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
147
148        self.xy_extent = array(xy_extent, Float)
149
150
151        # Allocate space for geometric quantities
152        self.normals = zeros((N, 6), Float)
153        self.areas = zeros(N, Float)
154        self.edgelengths = zeros((N, 3), Float)
155
156        # Get x,y coordinates for all triangles and store
157        self.vertex_coordinates = V = self.compute_vertex_coordinates()
158
159
160        # Initialise each triangle
161        if verbose:
162            print 'General_mesh: Computing areas, normals and edgelenghts'
163
164        for i in range(N):
165            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
166
167            x0, y0 = V[3*i, :]
168            x1, y1 = V[3*i+1, :]
169            x2, y2 = V[3*i+2, :]
170
171            # Area
172            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
173
174            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
175            msg += ' is degenerate:  area == %f' %self.areas[i]
176            assert self.areas[i] > 0.0, msg
177
178
179            # Normals
180            # The normal vectors
181            #   - point outward from each edge
182            #   - are orthogonal to the edge
183            #   - have unit length
184            #   - Are enumerated according to the opposite corner:
185            #     (First normal is associated with the edge opposite
186            #     the first vertex, etc)
187            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
188
189            n0 = array([x2 - x1, y2 - y1])
190            l0 = sqrt(sum(n0**2))
191
192            n1 = array([x0 - x2, y0 - y2])
193            l1 = sqrt(sum(n1**2))
194
195            n2 = array([x1 - x0, y1 - y0])
196            l2 = sqrt(sum(n2**2))
197
198            # Normalise
199            n0 /= l0
200            n1 /= l1
201            n2 /= l2
202
203            # Compute and store
204            self.normals[i, :] = [n0, -n0,
205                                  n1, -n1,
206                                  n2, -n2]
207
208            # Edgelengths
209            self.edgelengths[i, :] = [l0, l1, l2]
210
211
212        # Build structure listing which trianglse belong to which node.
213        if verbose: print 'Building inverted triangle structure'
214        self.build_inverted_triangle_structure()
215
216
217
218    def __len__(self):
219        return self.number_of_triangles
220
221
222    def __repr__(self):
223        return 'Mesh: %d vertices, %d triangles'\
224               %(self.nodes.shape, len(self))
225
226    def get_normals(self):
227        """Return all normal vectors.
228        Return normal vectors for all triangles as an Nx6 array
229        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
230        """
231        return self.normals
232
233
234    def get_normal(self, i, j):
235        """Return normal vector j of the i'th triangle.
236        Return value is the numeric array slice [x, y]
237        """
238        return self.normals[i, 2*j:2*j+2]
239
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 + array([self.geo_reference.get_xllcorner(),
277                                  self.geo_reference.get_yllcorner()])
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=array([self.geo_reference.get_xllcorner(),
319                                  self.geo_reference.get_yllcorner()])
320                return array([V[i3,:]+offset,
321                              V[i3+1,:]+offset,
322                              V[i3+2,:]+offset])
323            else:
324                return array([V[i3,:], V[i3+1,:], V[i3+2,:]])
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 = zeros((3*M, 2), 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            indices = range(M)
377
378        return take(self.triangles, indices)
379
380
381
382    def get_disconnected_triangles(self):
383        """Get mesh based on nodes obtained from get_vertex_coordinates.
384
385        Return array Mx3 array of integers where each row corresponds to
386        a triangle. A triangle is a triplet of indices into
387        point coordinates obtained from get_vertex_coordinates and each
388        index appears only once
389
390        This provides a mesh where no triangles share nodes
391        (hence the name disconnected triangles) and different
392        nodes may have the same coordinates.
393
394        This version of the mesh is useful for storing meshes with
395        discontinuities at each node and is e.g. used for storing
396        data in sww files.
397
398        The triangles created will have the format
399
400        [[0,1,2],
401         [3,4,5],
402         [6,7,8],
403         ...
404         [3*M-3 3*M-2 3*M-1]]
405        """
406
407        M = len(self) # Number of triangles
408        K = 3*M       # Total number of unique vertices
409
410        #T = reshape(array(range(K)).astype(Int), (M,3))
411        T = reshape(arange(K).astype(Int), (M,3))  # Faster
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 = 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 = array(triangle_list)
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 = 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 = sum(number_of_triangles_per_node)
538        vertex_value_indices = 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):