source: branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py @ 6883

Last change on this file since 6883 was 6843, checked in by rwilson, 16 years ago

Removed debug prints.

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