# source:anuga_core/source/anuga/pyvolution/general_mesh.py@3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 11.2 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 anuga.pyvolution 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                 geo_reference=None,
59                 verbose=False):
60        """
61        Build triangles from x,y coordinates (sequence of 2-tuples or
62        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
63        or Nx3 Numeric array of non-negative integers).
64
65        origin is a 3-tuple consisting of UTM zone, easting and northing.
66        If specified coordinates are assumed to be relative to this origin.
67        """
68
69        if verbose: print 'General_mesh: Building basic mesh structure'
70
71        self.triangles = array(triangles,Int)
72        self.coordinates = array(coordinates,Float)
73
74
75        # FIXME: this stores a geo_reference, but when coords are returned
76        # This geo_ref is not taken into account!
77        if geo_reference is None:
78            self.geo_reference = Geo_reference() #Use defaults
79        else:
80            self.geo_reference = geo_reference
81
82        #Input checks
83        msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples'
84        assert len(self.triangles.shape) == 2, msg
85
86        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
87        assert len(self.coordinates.shape) == 2, msg
88
89        msg = 'Vertex indices reference non-existing coordinate sets'
90        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
91
92
93        #Register number of elements (N)
94        self.number_of_elements = N = self.triangles.shape[0]
95
96        # FIXME: Maybe move to statistics?
97        # Or use with get_extent
98        xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
99                      max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
100
101        self.xy_extent = array(xy_extent, Float)
102
103
104        #Allocate space for geometric quantities
105        self.normals = zeros((N, 6), Float)
106        self.areas = zeros(N, Float)
107        self.edgelengths = zeros((N, 3), Float)
108
109        #Get x,y coordinates for all triangles and store
110        self.vertex_coordinates = V = self.compute_vertex_coordinates()
111
112
113        #Initialise each triangle
114        if verbose:
115            print 'General_mesh: Computing areas, normals and edgelenghts'
116
117        for i in range(N):
118            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
119
120
121            x0 = V[i, 0]; y0 = V[i, 1]
122            x1 = V[i, 2]; y1 = V[i, 3]
123            x2 = V[i, 4]; y2 = V[i, 5]
124
125            #Area
126            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
127
128            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
129            msg += ' is degenerate:  area == %f' %self.areas[i]
130            assert self.areas[i] > 0.0, msg
131
132
133            #Normals
134            #The normal vectors
135            # - point outward from each edge
136            # - are orthogonal to the edge
137            # - have unit length
138            # - Are enumerated according to the opposite corner:
139            #   (First normal is associated with the edge opposite
140            #    the first vertex, etc)
141            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
142
143            n0 = array([x2 - x1, y2 - y1])
144            l0 = sqrt(sum(n0**2))
145
146            n1 = array([x0 - x2, y0 - y2])
147            l1 = sqrt(sum(n1**2))
148
149            n2 = array([x1 - x0, y1 - y0])
150            l2 = sqrt(sum(n2**2))
151
152            #Normalise
153            n0 /= l0
154            n1 /= l1
155            n2 /= l2
156
157            #Compute and store
158            self.normals[i, :] = [n0[1], -n0[0],
159                                  n1[1], -n1[0],
160                                  n2[1], -n2[0]]
161
162            #Edgelengths
163            self.edgelengths[i, :] = [l0, l1, l2]
164
165
166        #Build vertex list
167        if verbose: print 'Building vertex list'
168        self.build_vertexlist()
169
170
171
172    def __len__(self):
173        return self.number_of_elements
174
175    def __repr__(self):
176        return 'Mesh: %d triangles, %d elements'\
177               %(self.coordinates.shape[0], len(self))
178
179    def get_normals(self):
180        """Return all normal vectors.
181        Return normal vectors for all triangles as an Nx6 array
182        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
183        """
184        return self.normals
185
186
187    def get_normal(self, i, j):
188        """Return normal vector j of the i'th triangle.
189        Return value is the numeric array slice [x, y]
190        """
191        return self.normals[i, 2*j:2*j+2]
192
193
194
195    def get_vertex_coordinates(self, obj=False, absolute=False):
196        """Return all vertex coordinates.
197        Return all vertex coordinates for all triangles as an Nx6 array
198        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
199
200        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
201        FIXME, we might make that the default.
202        FIXME Maybe use keyword: continuous = False for this condition?
203        FIXME - Maybe use something referring to unique vertices?
204
205        Boolean keyword argument absolute determines whether coordinates
206        are to be made absolute by taking georeference into account
207        Default is False as many parts of ANUGA expects relative coordinates.
208        (To see which, switch to default absolute=True and run tests).
209        """
210
211        V = self.vertex_coordinates
212        if absolute is True:
213            if not self.geo_reference.is_absolute():
214
215                V0 = self.geo_reference.get_absolute(V[:,0:2])
216                V1 = self.geo_reference.get_absolute(V[:,2:4])
217                V2 = self.geo_reference.get_absolute(V[:,4:6])
218
219                # This does double the memory need
220                V = concatenate( (V0, V1, V2), axis=1 )
221
222
223        if obj is True:
224
225            N = V.shape[0]
226            return reshape(V, (3*N, 2))
227        else:
228            return V
229
230
231    def get_vertex_coordinate(self, i, j, absolute=False):
232        """Return coordinates for vertex j of the i'th triangle.
233        Return value is the numeric array slice [x, y]
234        """
235
236        V = self.get_vertex_coordinates(absolute=absolute)
237        return V[i, 2*j:2*j+2]
238
239        ##return self.vertex_coordinates[i, 2*j:2*j+2]
240
241
242    def compute_vertex_coordinates(self):
243        """Return vertex coordinates for all triangles as an Nx6 array
244        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
245        """
246
247        #FIXME (Ole): Perhaps they should be ordered as in obj files??
248        #See quantity.get_vertex_values
249        #FIXME (Ole) - oh yes they should
250
251        N = self.number_of_elements
252        vertex_coordinates = zeros((N, 6), Float)
253
254        for i in range(N):
255            for j in range(3):
256                k = self.triangles[i,j]  #Index of vertex 0
257                v_k = self.coordinates[k]
258                vertex_coordinates[i, 2*j+0] = v_k[0]
259                vertex_coordinates[i, 2*j+1] = v_k[1]
260
261        return vertex_coordinates
262
263    def get_vertices(self, indices=None):
264        """Get connectivity
265        indices is the set of element ids of interest
266        """
267
268        if (indices ==  None):
269            indices = range(len(self))  #len(self)=number of elements
270
271        return  take(self.triangles, indices)
272
273    #FIXME - merge these two (get_vertices and get_triangles)
274    def get_triangles(self, obj=False):
275        """Get connetivity
276        Return triangles (triplets of indices into point coordinates)
277
278        If obj is True return structure commensurate with replicated
279        points, allowing for discontinuities
280        (FIXME: Need good name for this concept)
281        """
282
283        if obj is True:
284            m = len(self)  #Number of triangles
285            M = 3*m        #Total number of unique vertices
286            T = reshape(array(range(M)).astype(Int), (m,3))
287        else:
288            T = self.triangles
289
290        return T
291
292
293
294    def get_unique_vertices(self,  indices=None):
295        triangles = self.get_vertices(indices=indices)
296        unique_verts = {}
297        for triangle in triangles:
298           unique_verts[triangle[0]] = 0
299           unique_verts[triangle[1]] = 0
300           unique_verts[triangle[2]] = 0
301        return unique_verts.keys()
302
303    def build_vertexlist(self):
304        """Build vertexlist index by vertex ids and for each entry (point id)
305        build a list of (triangles, vertex_id) pairs that use the point
306        as vertex.
307
308        Preconditions:
309          self.coordinates and self.triangles are defined
310
311        Postcondition:
312          self.vertexlist is built
313        """
314
315        vertexlist = [None]*len(self.coordinates)
316        for i in range(self.number_of_elements):
317
318            a = self.triangles[i, 0]
319            b = self.triangles[i, 1]
320            c = self.triangles[i, 2]
321
322            #Register the vertices v as lists of
323            #(triangle_id, vertex_id) tuples associated with them
324            #This is used for averaging multiple vertex values.
325            for vertex_id, v in enumerate([a,b,c]):
326                if vertexlist[v] is None:
327                    vertexlist[v] = []
328
329                vertexlist[v].append( (i, vertex_id) )
330
331        self.vertexlist = vertexlist
332
333
334    def get_extent(self, absolute=False):
335        """Return min and max of all x and y coordinates
336
337        Boolean keyword argument absolute determines whether coordinates
338        are to be made absolute by taking georeference into account
339        """
340
341
342
343        C = self.get_vertex_coordinates(absolute=absolute)
344        X = C[:,0:6:2].copy()
345        Y = C[:,1:6:2].copy()
346
347        xmin = min(X.flat)
348        xmax = max(X.flat)
349        ymin = min(Y.flat)
350        ymax = max(Y.flat)
351
352        return xmin, xmax, ymin, ymax
353
354
355    def get_area(self):