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

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

Parallel domains now store only full triangles in sww files.
Still need to remove ghost nodes.

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 triangular elements (purely geometric)
9
10    A triangular element is defined in terms of three vertex ids,
11    ordered counter clock-wise,
12    each corresponding to a given coordinate set.
13    Vertices from different elements can point to the same
14    coordinate set.
15
16    Coordinate sets are implemented as an N x 2 Numeric array containing
17    x and y coordinates.
18
19
20    To instantiate:
21       Mesh(coordinates, triangles)
22
23    where
24
25      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
26      floats representing all x, y coordinates in the mesh.
27
28      triangles is either a list of 3-tuples or an Nx3 Numeric array of
29      integers representing indices of all vertices in the mesh.
30      Each vertex is identified by its index i in [0, M-1].
31
32
33    Example:
34        a = [0.0, 0.0]
35        b = [0.0, 2.0]
36        c = [2.0,0.0]
37        e = [2.0, 2.0]
38
39        points = [a, b, c, e]
40        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
41        mesh = Mesh(points, triangles)
42
43        #creates two triangles: bac and bce
44
45    Other:
46
47      In addition mesh computes an Nx6 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
52        This is a cut down version of mesh from mesh.py
53    """
54
55    #FIXME: It would be a good idea to use geospatial data as an alternative
56    #input
57    def __init__(self, coordinates, triangles,
58                 number_of_full_nodes=None,
59                 number_of_full_triangles=None,
60                 geo_reference=None,
61                 verbose=False):
62        """
63        Build triangles from x,y coordinates (sequence of 2-tuples or
64        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
65        or Nx3 Numeric array of non-negative integers).
66
67        origin is a 3-tuple consisting of UTM zone, easting and northing.
68        If specified coordinates are assumed to be relative to this origin.
69
70
71        number_of_full_nodes and number_of_full_triangles relate to
72        parallelism when each mesh has an extra layer of ghost points and
73        ghost triangles attached to the end of the two arrays.
74        In this case it is usefull to specify the number of real (called full)
75        nodes and triangles. If omitted they will default to all.
76        """
77
78        if verbose: print 'General_mesh: Building basic mesh structure'
79
80        self.triangles = array(triangles,Int)
81        self.coordinates = array(coordinates,Float)
82
83        # Register number of elements and nodes
84        self.number_of_triangles = N = self.triangles.shape[0]
85        self.number_of_nodes = self.coordinates.shape[0]
86
87
88
89        if number_of_full_nodes is None:
90            self.number_of_full_nodes=self.number_of_nodes
91        else:
92            assert int(number_of_full_nodes)
93            self.number_of_full_nodes=number_of_full_nodes
94
95
96        if number_of_full_triangles is None:
97            self.number_of_full_triangles=self.number_of_triangles
98        else:
99            assert int(number_of_full_triangles)
100            self.number_of_full_triangles=number_of_full_triangles
101
102
103        #print self.number_of_full_nodes, self.number_of_nodes
104        #print self.number_of_full_triangles, self.number_of_triangles
105
106
107
108        # FIXME: this stores a geo_reference, but when coords are returned
109        # This geo_ref is not taken into account!
110        if geo_reference is None:
111            self.geo_reference = Geo_reference() #Use defaults
112        else:
113            self.geo_reference = geo_reference
114
115        # Input checks
116        msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples. '
117        msg += 'The supplied array has the shape: %s'\
118               %str(self.triangles.shape)
119        assert len(self.triangles.shape) == 2, msg
120
121        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
122        msg += 'The supplied array has the shape: %s'\
123               %str(self.coordinates.shape)
124        assert len(self.coordinates.shape) == 2, msg
125
126        msg = 'Vertex indices reference non-existing coordinate sets'
127        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
128
129
130        # FIXME: Maybe move to statistics?
131        # Or use with get_extent
132        xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
133                      max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
134
135        self.xy_extent = array(xy_extent, Float)
136
137
138        # Allocate space for geometric quantities
139        self.normals = zeros((N, 6), Float)
140        self.areas = zeros(N, Float)
141        self.edgelengths = zeros((N, 3), Float)
142
143        # Get x,y coordinates for all triangles and store
144        self.vertex_coordinates = V = self.compute_vertex_coordinates()
145
146
147        # Initialise each triangle
148        if verbose:
149            print 'General_mesh: Computing areas, normals and edgelenghts'
150
151        for i in range(N):
152            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
153
154
155            x0 = V[i, 0]; y0 = V[i, 1]
156            x1 = V[i, 2]; y1 = V[i, 3]
157            x2 = V[i, 4]; y2 = V[i, 5]
158
159            # Area
160            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
161
162            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
163            msg += ' is degenerate:  area == %f' %self.areas[i]
164            assert self.areas[i] > 0.0, msg
165
166
167            # Normals
168            # The normal vectors
169            #   - point outward from each edge
170            #   - are orthogonal to the edge
171            #   - have unit length
172            #   - Are enumerated according to the opposite corner:
173            #     (First normal is associated with the edge opposite
174            #     the first vertex, etc)
175            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
176
177            n0 = array([x2 - x1, y2 - y1])
178            l0 = sqrt(sum(n0**2))
179
180            n1 = array([x0 - x2, y0 - y2])
181            l1 = sqrt(sum(n1**2))
182
183            n2 = array([x1 - x0, y1 - y0])
184            l2 = sqrt(sum(n2**2))
185
186            # Normalise
187            n0 /= l0
188            n1 /= l1
189            n2 /= l2
190
191            # Compute and store
192            self.normals[i, :] = [n0[1], -n0[0],
193                                  n1[1], -n1[0],
194                                  n2[1], -n2[0]]
195
196            # Edgelengths
197            self.edgelengths[i, :] = [l0, l1, l2]
198
199
200        # Build vertex list
201        if verbose: print 'Building vertex list'
202        self.build_vertexlist()
203
204
205
206    def __len__(self):
207        return self.number_of_triangles
208
209
210    def __repr__(self):
211        return 'Mesh: %d vertices, %d triangles'\
212               %(self.coordinates.shape[0], len(self))
213
214    def get_normals(self):
215        """Return all normal vectors.
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        return self.normals
220
221
222    def get_normal(self, i, j):
223        """Return normal vector j of the i'th triangle.
224        Return value is the numeric array slice [x, y]
225        """
226        return self.normals[i, 2*j:2*j+2]
227
228
229
230    def get_vertex_coordinates(self, unique=False, obj=False, absolute=False):
231        """Return all vertex coordinates.
232        Return all vertex coordinates for all triangles as an Nx6 array
233        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
234
235        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
236        FIXME, we might make that the default.
237        FIXME Maybe use keyword: continuous = False for this condition?
238        FIXME - Maybe use something referring to unique vertices?
239
240        Boolean keyword argument absolute determines whether coordinates
241        are to be made absolute by taking georeference into account
242        Default is False as many parts of ANUGA expects relative coordinates.
243        (To see which, switch to default absolute=True and run tests).
244        """
245
246        N = self.number_of_full_nodes
247
248        if unique is True:
249            V = self.coordinates[:N,:]
250            if absolute is True:
251                if not self.geo_reference.is_absolute():
252                    V = self.geo_reference.get_absolute(V)
253
254            return V
255
256
257
258        V = self.vertex_coordinates
259        if absolute is True:
260            if not self.geo_reference.is_absolute():
261
262                V0 = self.geo_reference.get_absolute(V[:N,0:2])
263                V1 = self.geo_reference.get_absolute(V[:N,2:4])
264                V2 = self.geo_reference.get_absolute(V[:N,4:6])
265
266                # This does double the memory need
267                V = concatenate( (V0, V1, V2), axis=1 )
268
269
270        if obj is True:
271            N = V.shape[0]
272            return reshape(V, (3*N, 2))
273        else:
274            return V
275
276
277    def get_vertex_coordinate(self, i, j, absolute=False):
278        """Return coordinates for vertex j of the i'th triangle.
279        Return value is the numeric array slice [x, y]
280        """
281
282        V = self.get_vertex_coordinates(absolute=absolute)
283        return V[i, 2*j:2*j+2]
284
285        ##return self.vertex_coordinates[i, 2*j:2*j+2]
286
287
288    def compute_vertex_coordinates(self):
289        """Return vertex coordinates for all triangles as an Nx6 array
290        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
291        """
292
293        #FIXME (Ole): Perhaps they should be ordered as in obj files??
294        #See quantity.get_vertex_values
295        #FIXME (Ole) - oh yes they should
296
297        N = self.number_of_triangles
298        vertex_coordinates = zeros((N, 6), Float)
299
300        for i in range(N):
301            for j in range(3):
302                k = self.triangles[i,j]  #Index of vertex 0
303                v_k = self.coordinates[k]
304                vertex_coordinates[i, 2*j+0] = v_k[0]
305                vertex_coordinates[i, 2*j+1] = v_k[1]
306
307        return vertex_coordinates
308
309    def get_vertices(self, indices=None):
310        """Get connectivity
311        indices is the set of element ids of interest
312        """
313
314        N = self.number_of_full_triangles
315
316        if indices is None:
317            #indices = range(len(self))  #len(self)=number of elements
318            indices = range(N)
319
320        return  take(self.triangles, indices)
321
322
323    #FIXME - merge these two (get_vertices and get_triangles)
324    def get_triangles(self, obj=False):
325        """Get connetivity
326        Return triangles (triplets of indices into point coordinates)
327
328        If obj is True return structure commensurate with replicated
329        points, allowing for discontinuities
330        (FIXME: Need good name for this concept)
331        """
332
333        if obj is True:
334            m = len(self)  #Number of triangles
335            M = 3*m        #Total number of unique vertices
336            T = reshape(array(range(M)).astype(Int), (m,3))
337        else:
338            T = self.triangles
339
340        return T
341
342
343
344    def get_unique_vertices(self,  indices=None):
345        triangles = self.get_vertices(indices=indices)
346        unique_verts = {}
347        for triangle in triangles:
348           unique_verts[triangle[0]] = 0
349           unique_verts[triangle[1]] = 0
350           unique_verts[triangle[2]] = 0
351        return unique_verts.keys()
352
353
354    def build_vertexlist(self):
355        """Build vertexlist index by vertex ids and for each entry (point id)
356        build a list of (triangles, vertex_id) pairs that use the point
357        as vertex.
358
359        Preconditions:
360          self.coordinates and self.triangles are defined
361
362        Postcondition:
363          self.vertexlist is built
364        """
365
366        vertexlist = [None]*len(self.coordinates)
367        for i in range(self.number_of_triangles):
368
369            a = self.triangles[i, 0]
370            b = self.triangles[i, 1]
371            c = self.triangles[i, 2]
372
373            #Register the vertices v as lists of
374            #(triangle_id, vertex_id) tuples associated with them
375            #This is used for averaging multiple vertex values.
376            for vertex_id, v in enumerate([a,b,c]):
377                if vertexlist[v] is None:
378                    vertexlist[v] = []
379
380                vertexlist[v].append( (i, vertex_id) )
381
382        self.vertexlist = vertexlist
383
384
385    def get_extent(self, absolute=False):
386        """Return min and max of all x and y coordinates
387
388        Boolean keyword argument absolute determines whether coordinates
389        are to be made absolute by taking georeference into account
390        """
391
392
393
394        C = self.get_vertex_coordinates(absolute=absolute)
395        X = C[:,0:6:2].copy()
396        Y = C[:,1:6:2].copy()
397
398        xmin = min(X.flat)
399        xmax = max(X.flat)
400        ymin = min(Y.flat)
401        ymax = max(Y.flat)
402
403        return xmin, xmax, ymin, ymax
404
405
406    def get_area(self):
407        """Return total area of mesh
408        """
409
410        return sum(self.areas)
411
412
Note: See TracBrowser for help on using the repository browser.