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

Last change on this file since 3945 was 3945, checked in by ole, 17 years ago

One large step towards major cleanup. This has mainly to do with
the way vertex coordinates are handled internally.

File size: 13.0 KB
Line
1
2from Numeric import concatenate, reshape, take, allclose
3from Numeric import array, zeros, Int, Float, sqrt, sum
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'
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[0]
98        self.number_of_nodes = self.nodes.shape[0]
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[0], 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[1], -n0[0],
205                                  n1[1], -n1[0],
206                                  n2[1], -n2[0]]
207
208            # Edgelengths
209            self.edgelengths[i, :] = [l0, l1, l2]
210
211
212        # Build vertex list
213        if verbose: print 'Building vertex list'
214        self.build_vertexlist()
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[0], 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 node coordinates ordered in an Nx2 array.
243        This is the same format they were provided in the constructor
244        i.e. without any duplication.
245
246        Boolean keyword argument absolute determines whether coordinates
247        are to be made absolute by taking georeference into account
248        Default is False as many parts of ANUGA expects relative coordinates.
249        (To see which, switch to default absolute=True and run tests).
250        """
251
252        N = self.number_of_full_nodes
253        V = self.nodes[:N,:]
254        if absolute is True:
255            if not self.geo_reference.is_absolute():
256                V = self.geo_reference.get_absolute(V)
257
258        return V
259
260
261
262    def get_vertex_coordinates(self, unique=False, absolute=False):
263        """Return all vertex coordinates.
264        Return all vertex coordinates for all triangles as a 3*N x 2 array
265        where the jth vertex of the ith triangle is located in row 3*i+j.
266
267        Boolean keyword unique will cause the points to be returned as
268        they were provided in the constructor i.e. without any duplication
269        in an N x 2 array.
270
271        """
272
273        if unique is True:
274            return self.get_nodes(absolute)
275
276
277        V = self.vertex_coordinates
278        if absolute is True:
279            if not self.geo_reference.is_absolute():
280                V = self.geo_reference.get_absolute(V)
281
282        return V
283
284
285
286    def get_vertex_coordinate(self, i, j, absolute=False):
287        """Return coordinates for vertex j of the i'th triangle.
288        Return value is the numeric array slice [x, y]
289        """
290
291        V = self.get_vertex_coordinates(absolute=absolute)
292        return V[3*i+j, :]
293
294
295    def compute_vertex_coordinates(self):
296        """Return all vertex coordinates for all triangles as a 3*N x 2 array
297        where the jth vertex of the ith triangle is located in row 3*i+j.
298
299        This function is used to precompute this important structure. Use
300        get_vertex coordinates to retrieve the points.
301        """
302
303        N = self.number_of_triangles
304        vertex_coordinates = zeros((3*N, 2), Float)
305
306        for i in range(N):
307            for j in range(3):
308                k = self.triangles[i,j]  #Index of vertex 0
309                v_k = self.nodes[k]
310
311                vertex_coordinates[3*i+j,:] = v_k
312
313
314        return vertex_coordinates
315
316
317    def get_vertices(self, indices=None):
318        """Get connectivity
319        indices is the set of element ids of interest
320        """
321
322        N = self.number_of_full_triangles
323
324        if indices is None:
325            #indices = range(len(self))  #len(self)=number of elements
326            indices = range(N)
327
328        return  take(self.triangles, indices)
329
330
331    #FIXME - merge these two (get_vertices and get_triangles)
332    def get_triangles(self, obj=False):
333        """Get connetivity
334        Return triangles (triplets of indices into point coordinates)
335
336        If obj is True return structure commensurate with replicated
337        points, allowing for discontinuities
338        (FIXME: Need good name for this concept)
339        """
340
341        if obj is True:
342            m = len(self)  #Number of triangles
343            M = 3*m        #Total number of unique vertices
344            T = reshape(array(range(M)).astype(Int), (m,3))
345        else:
346            T = self.triangles
347
348        return T
349
350
351
352    def get_unique_vertices(self,  indices=None):
353        triangles = self.get_vertices(indices=indices)
354        unique_verts = {}
355        for triangle in triangles:
356           unique_verts[triangle[0]] = 0
357           unique_verts[triangle[1]] = 0
358           unique_verts[triangle[2]] = 0
359        return unique_verts.keys()
360
361
362    def build_vertexlist(self):
363        """Build vertexlist index by vertex ids and for each entry (point id)
364        build a list of (triangles, vertex_id) pairs that use the point
365        as vertex.
366
367        Preconditions:
368          self.nodes and self.triangles are defined
369
370        Postcondition:
371          self.vertexlist is built
372        """
373
374        vertexlist = [None]*len(self.nodes)
375        for i in range(self.number_of_triangles):
376
377            a = self.triangles[i, 0]
378            b = self.triangles[i, 1]
379            c = self.triangles[i, 2]
380
381            #Register the vertices v as lists of
382            #(triangle_id, vertex_id) tuples associated with them
383            #This is used for averaging multiple vertex values.
384            for vertex_id, v in enumerate([a,b,c]):
385                if vertexlist[v] is None:
386                    vertexlist[v] = []
387
388                vertexlist[v].append( (i, vertex_id) )
389
390        self.vertexlist = vertexlist
391
392
393    def get_extent(self, absolute=False):
394        """Return min and max of all x and y coordinates
395
396        Boolean keyword argument absolute determines whether coordinates
397        are to be made absolute by taking georeference into account
398        """
399
400
401
402        C = self.get_vertex_coordinates(absolute=absolute)
403        X = C[:,0:6:2].copy()
404        Y = C[:,1:6:2].copy()
405
406        xmin = min(X.flat)
407        xmax = max(X.flat)
408        ymin = min(Y.flat)
409        ymax = max(Y.flat)
410
411        return xmin, xmax, ymin, ymax
412
413
414    def get_area(self):