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

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

Modified General_mesh.get_node(absolute=True) to return a copy if modifying the object node data.

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