source: inundation/ga/storm_surge/pyvolution/mesh.py @ 648

Last change on this file since 648 was 648, checked in by ole, 20 years ago

Refactoring boundary structure to allow None Boundary object for internal boundaries.

File size: 16.3 KB
Line 
1"""Classes implementing general 2D geometrical mesh of triangles.
2
3   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
4   Geoscience Australia, 2004   
5"""
6
7from general_mesh import General_mesh
8
9class Mesh(General_mesh):
10    """Collection of triangular elements (purely geometric)
11
12    A triangular element is defined in terms of three vertex ids,
13    ordered counter clock-wise,
14    each corresponding to a given coordinate set.
15    Vertices from different elements can point to the same
16    coordinate set.
17
18    Coordinate sets are implemented as an N x 2 Numeric array containing
19    x and y coordinates.
20   
21
22    To instantiate:
23       Mesh(coordinates, triangles)
24
25    where
26
27      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
28      floats representing all x, y coordinates in the mesh.
29
30      triangles is either a list of 3-tuples or an Nx3 Numeric array of
31      integers representing indices of all vertices in the mesh.
32      Each vertex is identified by its index i in [0, M-1].
33
34
35    Example:
36        a = [0.0, 0.0]
37        b = [0.0, 2.0]
38        c = [2.0,0.0]
39        e = [2.0, 2.0]
40       
41        points = [a, b, c, e]
42        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
43        mesh = Mesh(points, triangles)
44   
45        #creates two triangles: bac and bce
46
47
48    Mesh takes the optional third argument boundary which is a
49    dictionary mapping from (element_id, edge_id) to boundary tag.
50    The default value is None which will assign the defualt_boundary_tag
51    as specified in config.py to all boundary edges.
52    """
53
54    #FIXME: Maybe rename coordinates to points (as in a poly file)
55    #But keep 'vertex_coordinates'
56
57    #FIXME: Put in check for angles less than a set minimum
58   
59    def __init__(self, coordinates, triangles, boundary = None,
60                 tagged_elements = None):
61        """
62        Build triangles from x,y coordinates (sequence of 2-tuples or
63        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
64        or Nx3 Numeric array of non-negative integers).
65        """
66
67        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
68
69        General_mesh.__init__(self, coordinates, triangles)
70
71        N = self.number_of_elements
72       
73        #Allocate space for geometric quantities
74        self.centroid_coordinates = zeros((N, 2), Float)
75
76        self.radii = zeros(N, Float)
77
78        self.neighbours = zeros((N, 3), Int)
79        self.neighbour_edges = zeros((N, 3), Int)
80        self.number_of_boundaries = zeros(N, Int) 
81        self.surrogate_neighbours = zeros((N, 3), Int)       
82       
83        #Get x,y coordinates for all triangles and store
84        V = self.vertex_coordinates
85       
86        #Initialise each triangle
87        for i in range(N):
88            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
89           
90            x0 = V[i, 0]; y0 = V[i, 1]
91            x1 = V[i, 2]; y1 = V[i, 3]
92            x2 = V[i, 4]; y2 = V[i, 5]           
93
94            #Compute centroid
95            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
96            self.centroid_coordinates[i] = centroid
97
98            #Midpoints       
99            m0 = array([(x1 + x2)/2, (y1 + y2)/2])
100            m1 = array([(x0 + x2)/2, (y0 + y2)/2])
101            m2 = array([(x1 + x0)/2, (y1 + y0)/2])
102
103            #The radius is the distance from the centroid of
104            #a triangle to the midpoint of the side of the triangle
105            #closest to the centroid
106            d0 = sqrt(sum( (centroid-m0)**2 )) 
107            d1 = sqrt(sum( (centroid-m1)**2 ))       
108            d2 = sqrt(sum( (centroid-m2)**2 ))
109       
110            self.radii[i] = min(d0, d1, d2)               
111
112            #Initialise Neighbours (-1 means that it is a boundary neighbour)
113            self.neighbours[i, :] = [-1, -1, -1]
114
115            #Initialise edge ids of neighbours
116            #In case of boundaries this slot is not used
117            self.neighbour_edges[i, :] = [-1, -1, -1]
118
119
120        #Build neighbour structure   
121        self.build_neighbour_structure()
122
123        #Build surrogate neighbour structure   
124        self.build_surrogate_neighbour_structure()       
125
126        #Build boundary dictionary mapping (id, edge) to symbolic tags       
127        self.build_boundary_dictionary(boundary)
128
129        #Build tagged element  dictionary mapping (tag) to array of elements
130        self.build_tagged_elements_dictionary(tagged_elements)
131       
132        #Update boundary indices
133        #self.build_boundary_structure()       
134
135       
136    def __repr__(self):
137        return 'Mesh: %d triangles, %d elements, %d boundary segments'\
138               %(self.coordinates.shape[0], len(self), len(self.boundary))
139
140           
141    def build_neighbour_structure(self):
142        """Update all registered triangles to point to their neighbours.
143       
144        Also, keep a tally of the number of boundaries for each triangle
145
146        Postconditions:
147          neighbours and neighbour_edges is populated
148          number_of_boundaries integer array is defined.
149        """
150       
151        #Step 1:
152        #Build dictionary mapping from segments (2-tuple of points)
153        #to left hand side edge (facing neighbouring triangle)
154
155        N = self.number_of_elements       
156        neighbourdict = {}
157        for i in range(N):
158
159            #Register all segments as keys mapping to current triangle
160            #and segment id
161            a = self.triangles[i, 0]
162            b = self.triangles[i, 1]
163            c = self.triangles[i, 2]           
164            neighbourdict[a,b] = (i, 2) #(id, edge)
165            neighbourdict[b,c] = (i, 0) #(id, edge)
166            neighbourdict[c,a] = (i, 1) #(id, edge)
167
168
169        #Step 2:
170        #Go through triangles again, but this time
171        #reverse direction of segments and lookup neighbours.
172        for i in range(N):
173            a = self.triangles[i, 0]
174            b = self.triangles[i, 1]
175            c = self.triangles[i, 2]
176
177            self.number_of_boundaries[i] = 3
178            if neighbourdict.has_key((b,a)):
179                self.neighbours[i, 2] = neighbourdict[b,a][0]
180                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
181                self.number_of_boundaries[i] -= 1
182               
183            if neighbourdict.has_key((c,b)):
184                self.neighbours[i, 0] = neighbourdict[c,b][0]
185                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
186                self.number_of_boundaries[i] -= 1               
187               
188            if neighbourdict.has_key((a,c)):
189                self.neighbours[i, 1] = neighbourdict[a,c][0]
190                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
191                self.number_of_boundaries[i] -= 1             
192
193
194    def build_surrogate_neighbour_structure(self):
195        """Build structure where each triangle edge points to its neighbours
196        if they exist. Otherwise point to the triangle itself.
197
198        The surrogate neighbour structure is useful for computing gradients
199        based on centroid values of neighbours.
200
201        Precondition: Neighbour structure is defined
202        Postcondition:
203          Surrogate neighbour structure is defined:
204          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
205          triangles.
206
207        """
208
209        N = self.number_of_elements
210        for i in range(N):
211            #Find all neighbouring volumes that are not boundaries
212            for k in range(3):
213                if self.neighbours[i, k] < 0:
214                    self.surrogate_neighbours[i, k] = i #Point this triangle
215                else:
216                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
217
218
219
220    def build_boundary_dictionary(self, boundary = None):
221        """Build or check the dictionary of boundary tags.
222         self.boundary is a dictionary of tags,
223         keyed by volume id and edge:
224         { (id, edge): tag, ... }
225         
226         Postconditions:
227            self.boundary is defined.
228        """
229
230        from config import default_boundary_tag
231       
232        if boundary is None:
233            boundary = {}
234            for vol_id in range(self.number_of_elements):
235                for edge_id in range(0, 3):
236                    if self.neighbours[vol_id, edge_id] < 0:
237                        boundary[(vol_id, edge_id)] = default_boundary_tag
238        else:
239            #Check that all keys in given boundary exist
240            for vol_id, edge_id in boundary.keys():
241                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
242                a, b = self.neighbours.shape
243                assert vol_id < a and edge_id < b, msg
244
245                #FIXME: This assert violates internal boundaries (delete it)
246                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
247                #assert self.neighbours[vol_id, edge_id] < 0, msg
248
249            #Check that all boundary segments are assigned a tag
250            for vol_id in range(self.number_of_elements):
251                for edge_id in range(0, 3):
252                    if self.neighbours[vol_id, edge_id] < 0:
253                        if not boundary.has_key( (vol_id, edge_id) ):
254                            msg = 'WARNING: Given boundary does not contain '
255                            msg += 'tags for edge (%d, %d). '\
256                                   %(vol_id, edge_id)
257                            msg += 'Assigning default tag (%s).'\
258                                   %default_boundary_tag
259
260                            #FIXME: Print only as per verbosity
261                            #print msg
262
263                            #FIXME: Make this situation an error in the future
264                            #and make another function which will
265                            #enable default boundary-tags where
266                            #tags a not specified
267                            boundary[ (vol_id, edge_id) ] =\
268                                      default_boundary_tag
269           
270
271           
272        self.boundary = boundary
273           
274
275    def build_tagged_elements_dictionary(self, tagged_elements = None):
276        """Build the dictionary of element tags.
277         self.tagged_elements is a dictionary of element arrays,
278         keyed by tag:
279         { (tag): [e1, e2, e3..] }
280         
281         Postconditions:
282            self.element_tag is defined
283        """
284        from Numeric import array, Int
285       
286        if tagged_elements is None:
287            tagged_elements = {}
288        else:
289            #Check that all keys in given boundary exist
290            for tag in tagged_elements.keys():
291                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
292               
293                msg = 'Not all elements exist. '
294                assert max(tagged_elements[tag]) < self.number_of_elements, msg
295        #print "tagged_elements", tagged_elements               
296        self.tagged_elements = tagged_elements
297           
298    def build_boundary_structure(self):
299        """Traverse boundary and
300        enumerate neighbour indices from -1 and
301        counting down.
302
303        Precondition:
304            self.boundary is defined.
305        Post condition:
306            neighbour array has unique negative indices for boundary
307            boundary_segments array imposes an ordering on segments
308            (not otherwise available from the dictionary)
309
310        Note: If a segment is listed in the boundary dictionary
311        it *will* become a boundary - even if there is a neighbouring triangle.
312        This would be the case for internal boundaries
313        """
314
315        #FIXME: Maybe obsolete
316
317        if self.boundary is None:
318            msg = 'Boundary dictionary must be defined before '
319            msg += 'building boundary structure'
320            raise msg
321
322       
323        self.boundary_segments = self.boundary.keys()
324        self.boundary_segments.sort()
325
326        index = -1       
327        for id, edge in self.boundary_segments:
328
329            #FIXME: One would detect internal boundaries as follows
330            #if self.neighbours[id, edge] > -1:
331            #    print 'Internal boundary'
332           
333            self.neighbours[id, edge] = index
334            index -= 1
335
336   
337    def get_boundary_tags(self):
338        """Return list of available boundary tags
339        """
340       
341        tags = {}
342        for v in self.boundary.values():
343            tags[v] = 1
344
345        return tags.keys()   
346           
347   
348
349    def check_integrity(self):
350        """Check that triangles are internally consistent e.g.
351        that area corresponds to edgelengths, that vertices
352        are arranged in a counter-clockwise order, etc etc
353        Neighbour structure will be checked by class Mesh
354        """
355       
356        from config import epsilon
357        from math import pi
358        from util import anglediff
359
360        N = self.number_of_elements
361       
362        #Get x,y coordinates for all vertices for all triangles
363        V = self.get_vertex_coordinates()
364
365        #Check each triangle
366        for i in range(N):
367            x0 = V[i, 0]; y0 = V[i, 1]
368            x1 = V[i, 2]; y1 = V[i, 3]
369            x2 = V[i, 4]; y2 = V[i, 5]           
370
371            #Check that area hasn't been compromised
372            area = self.areas[i]
373            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
374            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
375                  %(x0,y0,x1,y1,x2,y2)
376            assert abs((area - ref)/area) < epsilon, msg
377       
378            #Check that points are arranged in counter clock-wise order
379            v0 = [x1-x0, y1-y0]
380            v1 = [x2-x1, y2-y1]
381            v2 = [x0-x2, y0-y2]
382
383            a0 = anglediff(v1, v0)
384            a1 = anglediff(v2, v1)
385            a2 = anglediff(v0, v2)       
386
387            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
388            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
389            assert a0 < pi and a1 < pi and a2 < pi, msg
390
391            #Check that normals are orthogonal to edge vectors
392            #Note that normal[k] lies opposite vertex k
393
394            normal0 = self.normals[i, 0:2]
395            normal1 = self.normals[i, 2:4]
396            normal2 = self.normals[i, 4:6]           
397
398            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
399                assert u[0]*v[0] + u[1]*v[1] < epsilon
400
401
402        #Check that all vertices have been registered
403        for v in self.vertexlist:
404            msg = 'Some points do not belong to an element.\n'
405            msg += 'Make sure all points appear as element vertices!'
406            if v is None:
407                print 'WARNING (mesh.py): %s' %msg
408                break
409
410            #Should this warning be an exception?
411            #assert v is not None, msg
412
413               
414        #Check integrity of neighbour structure
415        for i in range(N):
416            for v in self.triangles[i, :]:
417                #Check that all vertices have been registered
418                assert self.vertexlist[v] is not None
419
420                #Check that this triangle is listed with at least one vertex
421                assert (i, 0) in self.vertexlist[v] or\
422                       (i, 1) in self.vertexlist[v] or\
423                       (i, 2) in self.vertexlist[v]
424               
425               
426
427            #Check neighbour structure
428            for k, neighbour_id in enumerate(self.neighbours[i,:]):
429
430                #Assert that my neighbour's neighbour is me
431                #Boundaries need not fulfill this
432                if neighbour_id >= 0:
433                    edge = self.neighbour_edges[i, k]
434                    assert self.neighbours[neighbour_id, edge] == i
435
436
437
438        #Check that all boundaries have
439        # unique, consecutive, negative indices                   
440
441        #L = len(self.boundary)
442        #for i in range(L):
443        #    id, edge = self.boundary_segments[i]
444        #    assert self.neighbours[id, edge] == -i-1
445
446
447        #NOTE: This assert doesn't hold true if there are internal boundaries
448        #FIXME: Look into this further.
449        #FIXME (Ole): In pyvolution mark 3 this is OK again
450        #NOTE: No longer works because neighbour structure is modified by
451        #      domain set_boundary.
452        #for id, edge in self.boundary:
453        #    assert self.neighbours[id,edge] < 0
454
455
456       
457           
458    def get_centroid_coordinates(self):
459        """Return all centroid coordinates.
460        Return all centroid coordinates for all triangles as an Nx2 array
461        (ordered as x0, y0 for each triangle)       
462        """
463        return self.centroid_coordinates
Note: See TracBrowser for help on using the repository browser.