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

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

and get_edge_midpoint_coordinate into the general mesh. This is so that
boundary classes can get access to the location of the edge midpoint.

I will use this to allow Dirichlet_boundary to take a function of
time and space.

File size: 23.2 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_nodes(self):
238        return self.number_of_nodes
239
240    def get_nodes(self, absolute=False):
241        """Return all nodes in mesh.
242
243        The nodes are ordered in an Nx2 array where N is the number of nodes.
244        This is the same format they were provided in the constructor
245        i.e. without any duplication.
246
247        Boolean keyword argument absolute determines whether coordinates
248        are to be made absolute by taking georeference into account
249        Default is False as many parts of ANUGA expects relative coordinates.
250        (To see which, switch to default absolute=True and run tests).
251        """
252
253        N = self.number_of_full_nodes
254        V = self.nodes[:N,:]
255        if absolute is True:
256            if not self.geo_reference.is_absolute():
257                V = self.geo_reference.get_absolute(V)
258
259        return V
260
261    def get_node(self, i, absolute=False):
262        """Return node coordinates for triangle i.
263
264        Boolean keyword argument absolute determines whether coordinates
265        are to be made absolute by taking georeference into account
266        Default is False as many parts of ANUGA expects relative coordinates.
267        (To see which, switch to default absolute=True and run tests).
268
269        Note: This method returns a modified _copy_ of the nodes slice if
270              absolute is True.  If absolute is False, just return the slice.
271              This is related to the ensure_numeric() returning a copy problem.
272        """
273
274        V = self.nodes[i,:]
275        if absolute is True:
276            if not self.geo_reference.is_absolute():
277                # get a copy so as not to modify the internal self.nodes array
278                V = copy.copy(V)
279                V += num.array([self.geo_reference.get_xllcorner(),
280                                self.geo_reference.get_yllcorner()], num.float)
281        return V
282
283    def get_vertex_coordinates(self, triangle_id=None, absolute=False):
284        """Return vertex coordinates for all triangles.
285
286        Return all vertex coordinates for all triangles as a 3*M x 2 array
287        where the jth vertex of the ith triangle is located in row 3*i+j and
288        M the number of triangles in the mesh.
289
290        if triangle_id is specified (an integer) the 3 vertex coordinates
291        for triangle_id are returned.
292
293        Boolean keyword argument absolute determines whether coordinates
294        are to be made absolute by taking georeference into account
295        Default is False as many parts of ANUGA expects relative coordinates.
296        """
297
298        V = self.vertex_coordinates
299
300        if triangle_id is None:
301            if absolute is True:
302                if not self.geo_reference.is_absolute():
303                    V = self.geo_reference.get_absolute(V)
304            return V
305        else:
306            i = triangle_id
307            msg = 'triangle_id must be an integer'
308            assert int(i) == i, msg
309            assert 0 <= i < self.number_of_triangles
310
311            i3 = 3*i
312            if absolute is True and not self.geo_reference.is_absolute():
313                offset=num.array([self.geo_reference.get_xllcorner(),
314                                  self.geo_reference.get_yllcorner()], num.float)
315                return num.array([V[i3,:]+offset,
316                                  V[i3+1,:]+offset,
317                                  V[i3+2,:]+offset], num.float)
318            else:
319                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
320
321    def get_vertex_coordinate(self, i, j, absolute=False):
322        """Return coordinates for vertex j of the i'th triangle.
323        Return value is the numeric array slice [x, y]
324        """
325
326        msg = 'vertex id j must be an integer in [0,1,2]'
327        assert j in [0,1,2], msg
328
329        V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute)
330        return V[j,:]
331
332
333
334    def compute_vertex_coordinates(self):
335        """Return all vertex coordinates for all triangles as a 3*M x 2 array
336        where the jth vertex of the ith triangle is located in row 3*i+j.
337
338        This function is used to precompute this important structure. Use
339        get_vertex coordinates to retrieve the points.
340        """
341
342        M = self.number_of_triangles
343        vertex_coordinates = num.zeros((3*M, 2), num.float)
344
345        for i in range(M):
346            for j in range(3):
347                k = self.triangles[i,j] # Index of vertex j in triangle i
348                vertex_coordinates[3*i+j,:] = self.nodes[k]
349
350        return vertex_coordinates
351
352
353    def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
354        """Return edge midpoint coordinates for all triangles or from particular triangle.
355
356        Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
357        where the jth midpoint of the ith triangle is located in row 3*i+j and
358        M the number of triangles in the mesh.
359
360        if triangle_id is specified (an integer) the 3 midpoint coordinates
361        for triangle_id are returned.
362
363        Boolean keyword argument absolute determines whether coordinates
364        are to be made absolute by taking georeference into account
365        Default is False as many parts of ANUGA expects relative coordinates.
366        """
367
368        E = self.edge_midpoint_coordinates
369
370        if triangle_id is None:
371            if absolute is True:
372                if not self.geo_reference.is_absolute():
373                    E = self.geo_reference.get_absolute(E)
374            return E
375        else:
376            i = triangle_id
377            msg = 'triangle_id must be an integer'
378            assert int(i) == i, msg
379            assert 0 <= i < self.number_of_triangles
380
381            i3 = 3*i
382            if absolute is True and not self.geo_reference.is_absolute():
383                offset=num.array([self.geo_reference.get_xllcorner(),
384                                  self.geo_reference.get_yllcorner()], num.float)
385                return num.array([E[i3,:]+offset,
386                                  E[i3+1,:]+offset,
387                                  E[i3+2,:]+offset], num.float)
388            else:
389                return num.array([E[i3,:], E[i3+1,:], E[i3+2,:]], num.float)
390
391
392    def get_edge_midpoint_coordinate(self, i, j, absolute=False):
393        """Return coordinates for edge midpoint j of the i'th triangle.
394        Return value is the numeric array slice [x, y]
395        """
396
397        msg = 'edge midpoint id j must be an integer in [0,1,2]'
398        assert j in [0,1,2], msg
399
400        E = self.get_edge_midpoint_coordinates(triangle_id=i, absolute=absolute)
401        return E[j,:]
402
403
404    def compute_edge_midpoint_coordinates(self):
405        """Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
406        where the jth edge midpoint of the ith triangle is located in row 3*i+j.
407
408        This function is used to precompute this important structure. Use
409        get_edge_midpoint_coordinates to retrieve the points.
410
411        Assumes that vertex_coordinates have been computed
412        """
413
414        M = self.number_of_triangles
415        E = num.zeros((3*M, 2), num.float)
416
417        V = self.vertex_coordinates
418
419        V0 = V[0:3*M:3, :]
420        V1 = V[1:3*M:3, :]
421        V2 = V[2:3*M:3, :]
422
423
424        #print V.shape, V0.shape, V1.shape, V2.shape
425
426        #print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
427        E[0:3*M:3, :] = 0.5*(V1+V2)
428        E[1:3*M:3, :] = 0.5*(V2+V0)
429        E[2:3*M:3, :] = 0.5*(V0+V1)
430
431        return E
432
433
434
435    def get_triangles(self, indices=None):
436        """Get mesh triangles.
437
438        Return Mx3 integer array where M is the number of triangles.
439        Each row corresponds to one triangle and the three entries are
440        indices into the mesh nodes which can be obtained using the method
441        get_nodes()
442
443        Optional argument, indices is the set of triangle ids of interest.
444        """
445
446        M = self.number_of_full_triangles
447
448        if indices is None:
449            return self.triangles
450
451        return num.take(self.triangles, indices, axis=0)
452
453    def get_disconnected_triangles(self):
454        """Get mesh based on nodes obtained from get_vertex_coordinates.
455
456        Return array Mx3 array of integers where each row corresponds to
457        a triangle. A triangle is a triplet of indices into
458        point coordinates obtained from get_vertex_coordinates and each
459        index appears only once
460
461        This provides a mesh where no triangles share nodes
462        (hence the name disconnected triangles) and different
463        nodes may have the same coordinates.
464
465        This version of the mesh is useful for storing meshes with
466        discontinuities at each node and is e.g. used for storing
467        data in sww files.
468
469        The triangles created will have the format
470        [[0,1,2],
471         [3,4,5],
472         [6,7,8],
473         ...
474         [3*M-3 3*M-2 3*M-1]]
475        """
476
477        M = len(self) # Number of triangles
478        K = 3*M       # Total number of unique vertices
479        return num.reshape(num.arange(K, dtype=num.int), (M,3))
480
481    def get_unique_vertices(self, indices=None):
482        """FIXME(Ole): This function needs a docstring"""
483
484        triangles = self.get_triangles(indices=indices)
485        unique_verts = {}
486        for triangle in triangles:
487            unique_verts[triangle] = 0
488            unique_verts[triangle] = 0
489            unique_verts[triangle] = 0
490        return unique_verts.keys()
491
492    def get_triangles_and_vertices_per_node(self, node=None):
493        """Get triangles associated with given node.
494
495        Return list of triangle_ids, vertex_ids for specified node.
496        If node in None or absent, this information will be returned
497        for all (full) nodes in a list L where L[v] is the triangle
498        list for node v.
499        """
500
501        triangle_list = []
502        if node is not None:
503            # Get index for this node
504            first = num.sum(self.number_of_triangles_per_node[:node])
505
506            # Get number of triangles for this node
507            count = self.number_of_triangles_per_node[node]
508
509            for i in range(count):
510                index = self.vertex_value_indices[first+i]
511
512                volume_id = index / 3
513                vertex_id = index % 3
514
515                triangle_list.append( (volume_id, vertex_id) )
516
517            triangle_list = num.array(triangle_list, num.int)    #array default#
518        else:
519            # Get info for all nodes recursively.
520            # If need be, we can speed this up by
521            # working directly with the inverted triangle structure
522            for i in range(self.number_of_full_nodes):
523                L = self.get_triangles_and_vertices_per_node(node=i)
524                triangle_list.append(L)
525
526        return triangle_list
527
528    def build_inverted_triangle_structure(self):
529        """Build structure listing triangles belonging to each node
530
531        Two arrays are created and store as mesh attributes
532
533        number_of_triangles_per_node: An integer array of length N
534        listing for each node how many triangles use it. N is the number of
535        nodes in mesh.
536
537        vertex_value_indices: An array of length M listing indices into
538        triangles ordered by node number. The (triangle_id, vertex_id)
539        pairs are obtained from each index as (index/3, index%3) or each
540        index can be used directly into a flattened triangles array. This
541        is for example the case in the quantity.c where this structure is
542        used to average vertex values efficiently.
543
544        Example:
545        a = [0.0, 0.0] # node 0
546        b = [0.0, 2.0] # node 1
547        c = [2.0, 0.0] # node 2
548        d = [0.0, 4.0] # node 3
549        e = [2.0, 2.0] # node 4
550        f = [4.0, 0.0] # node 5
551        nodes = array([a, b, c, d, e, f])
552
553        #                    bac,     bce,     ecf,     dbe
554        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
555
556        For this structure:
557        number_of_triangles_per_node = [1 3 3 1 3 1]
558        which means that node a has 1 triangle associated with it, node b
559        has 3, node has 3 and so on.
560
561        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
562        which reflects the fact that
563        node 0 is used by triangle 0, vertex 1 (index = 1)
564        node 1 is used by triangle 0, vertex 0 (index = 0)
565                   and by triangle 1, vertex 0 (index = 3)
566                   and by triangle 3, vertex 1 (index = 10)
567        node 2 is used by triangle 0, vertex 2 (index = 2)
568                   and by triangle 1, vertex 1 (index = 4)
569                   and by triangle 2, vertex 1 (index = 7)
570        node 3 is used by triangle 3, vertex 0 (index = 9)
571        node 4 is used by triangle 1, vertex 2 (index = 5)
572                   and by triangle 2, vertex 0 (index = 6)
573                   and by triangle 3, vertex 2 (index = 11)
574        node 5 is used by triangle 2, vertex 2 (index = 8)
575
576        Preconditions:
577          self.nodes and self.triangles are defined
578
579        Postcondition:
580          self.number_of_triangles_per_node is built
581          self.vertex_value_indices is built
582        """
583
584        # Count number of triangles per node
585        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes,
586                                                 num.int)       #array default#
587        for volume_id, triangle in enumerate(self.get_triangles()):
588            for vertex_id in triangle:
589                number_of_triangles_per_node[vertex_id] += 1
590
591        # Allocate space for inverted structure
592        number_of_entries = num.sum(number_of_triangles_per_node)
593        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
594
595        # Register (triangle, vertex) indices for each node
596        vertexlist = [None] * self.number_of_full_nodes
597        for volume_id in range(self.number_of_full_triangles):
598            a = self.triangles[volume_id, 0]
599            b = self.triangles[volume_id, 1]
600            c = self.triangles[volume_id, 2]
601
602            for vertex_id, node_id in enumerate([a, b, c]):
603                if vertexlist[node_id] is None:
604                    vertexlist[node_id] = []
605                vertexlist[node_id].append((volume_id, vertex_id))
606
607        # Build inverted triangle index array
608        k = 0
609        for vertices in vertexlist:
610            if vertices is not None:
611                for volume_id, vertex_id in vertices:
612                    vertex_value_indices[k] = 3*volume_id + vertex_id
613                    k += 1
614
615        # Save structure
616        self.number_of_triangles_per_node = number_of_triangles_per_node
617        self.vertex_value_indices = vertex_value_indices
618
619    def get_extent(self, absolute=False):
620        """Return min and max of all x and y coordinates
621
622        Boolean keyword argument absolute determines whether coordinates
623        are to be made absolute by taking georeference into account
624        """
625
626        C = self.get_vertex_coordinates(absolute=absolute)
627        X = C[:,0:6:2].copy()
628        Y = C[:,1:6:2].copy()
629
630        xmin = num.min(X)
631        xmax = num.max(X)
632        ymin = num.min(Y)
633        ymax = num.max(Y)
634
635        return xmin, xmax, ymin, ymax
636
637    def get_areas(self):
638        """Get areas of all individual triangles."""
639
640        return self.areas
641
642    def get_area(self):