source: inundation/pyvolution/mesh.py @ 1794

Last change on this file since 1794 was 1670, checked in by ole, 19 years ago

Refactored Interpolation_function and tested

File size: 21.2 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 default_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
60    def __init__(self, coordinates, triangles, boundary = None,
61                 tagged_elements = None, geo_reference = None,
62                 use_inscribed_circle = False):
63        """
64        Build triangles from x,y coordinates (sequence of 2-tuples or
65        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
66        or Nx3 Numeric array of non-negative integers).
67        """
68
69
70
71        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
72
73        General_mesh.__init__(self, coordinates, triangles, geo_reference)
74
75        N = self.number_of_elements
76
77        self.use_inscribed_circle = use_inscribed_circle
78
79        #Allocate space for geometric quantities
80        self.centroid_coordinates = zeros((N, 2), Float)
81
82        self.radii = zeros(N, Float)
83
84        self.neighbours = zeros((N, 3), Int)
85        self.neighbour_edges = zeros((N, 3), Int)
86        self.number_of_boundaries = zeros(N, Int)
87        self.surrogate_neighbours = zeros((N, 3), Int)
88
89        #Get x,y coordinates for all triangles and store
90        V = self.vertex_coordinates
91
92        #Initialise each triangle
93        for i in range(N):
94            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
95
96            x0 = V[i, 0]; y0 = V[i, 1]
97            x1 = V[i, 2]; y1 = V[i, 3]
98            x2 = V[i, 4]; y2 = V[i, 5]
99
100            #Compute centroid
101            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
102            self.centroid_coordinates[i] = centroid
103
104
105            if self.use_inscribed_circle == False:
106                #OLD code. Computed radii may exceed that of an
107                #inscribed circle
108
109                #Midpoints
110                m0 = array([(x1 + x2)/2, (y1 + y2)/2])
111                m1 = array([(x0 + x2)/2, (y0 + y2)/2])
112                m2 = array([(x1 + x0)/2, (y1 + y0)/2])
113
114                #The radius is the distance from the centroid of
115                #a triangle to the midpoint of the side of the triangle
116                #closest to the centroid
117                d0 = sqrt(sum( (centroid-m0)**2 ))
118                d1 = sqrt(sum( (centroid-m1)**2 ))
119                d2 = sqrt(sum( (centroid-m2)**2 ))
120
121                self.radii[i] = min(d0, d1, d2)
122
123            else:
124                #NEW code added by Peter Row. True radius
125                #of inscribed circle is computed
126
127                from math import sqrt
128                a = sqrt((x0-x1)**2+(y0-y1)**2)
129                b = sqrt((x1-x2)**2+(y1-y2)**2)
130                c = sqrt((x2-x0)**2+(y2-y0)**2)
131
132                self.radii[i]=2.0*self.areas[i]/(a+b+c)
133
134
135            #Initialise Neighbours (-1 means that it is a boundary neighbour)
136            self.neighbours[i, :] = [-1, -1, -1]
137
138            #Initialise edge ids of neighbours
139            #In case of boundaries this slot is not used
140            self.neighbour_edges[i, :] = [-1, -1, -1]
141
142
143        #Build neighbour structure
144        self.build_neighbour_structure()
145
146        #Build surrogate neighbour structure
147        self.build_surrogate_neighbour_structure()
148
149        #Build boundary dictionary mapping (id, edge) to symbolic tags
150        self.build_boundary_dictionary(boundary)
151
152        #Build tagged element  dictionary mapping (tag) to array of elements
153        self.build_tagged_elements_dictionary(tagged_elements)
154
155        #Update boundary indices FIXME: OBSOLETE
156        #self.build_boundary_structure()
157
158        #FIXME check integrity?
159
160    def __repr__(self):
161        return 'Mesh: %d triangles, %d elements, %d boundary segments'\
162               %(self.coordinates.shape[0], len(self), len(self.boundary))
163
164
165    def set_to_inscribed_circle(self,safety_factor = 1):
166        #FIXME phase out eventually
167        from math import sqrt
168        N = self.number_of_elements
169        V = self.vertex_coordinates
170
171        #initialising min and max ratio
172        i=0
173        old_rad = self.radii[i]
174        x0 = V[i, 0]; y0 = V[i, 1]
175        x1 = V[i, 2]; y1 = V[i, 3]
176        x2 = V[i, 4]; y2 = V[i, 5]
177        a = sqrt((x0-x1)**2+(y0-y1)**2)
178        b = sqrt((x1-x2)**2+(y1-y2)**2)
179        c = sqrt((x2-x0)**2+(y2-y0)**2)
180        ratio = old_rad/self.radii[i]
181        max_ratio = ratio
182        min_ratio = ratio
183
184        for i in range(N):
185            old_rad = self.radii[i]
186            x0 = V[i, 0]; y0 = V[i, 1]
187            x1 = V[i, 2]; y1 = V[i, 3]
188            x2 = V[i, 4]; y2 = V[i, 5]
189            a = sqrt((x0-x1)**2+(y0-y1)**2)
190            b = sqrt((x1-x2)**2+(y1-y2)**2)
191            c = sqrt((x2-x0)**2+(y2-y0)**2)
192            self.radii[i]=self.areas[i]/(2*(a+b+c))*safety_factor
193            ratio = old_rad/self.radii[i]
194            if ratio >= max_ratio: max_ratio = ratio
195            if ratio <= min_ratio: min_ratio = ratio
196        return max_ratio,min_ratio
197
198    def build_neighbour_structure(self):
199        """Update all registered triangles to point to their neighbours.
200
201        Also, keep a tally of the number of boundaries for each triangle
202
203        Postconditions:
204          neighbours and neighbour_edges is populated
205          number_of_boundaries integer array is defined.
206        """
207
208        #Step 1:
209        #Build dictionary mapping from segments (2-tuple of points)
210        #to left hand side edge (facing neighbouring triangle)
211
212        N = self.number_of_elements
213        neighbourdict = {}
214        for i in range(N):
215
216            #Register all segments as keys mapping to current triangle
217            #and segment id
218            a = self.triangles[i, 0]
219            b = self.triangles[i, 1]
220            c = self.triangles[i, 2]
221            if neighbourdict.has_key((a,b)):
222                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
223                    raise msg
224            if neighbourdict.has_key((b,c)):
225                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
226                    raise msg
227            if neighbourdict.has_key((c,a)):
228                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
229                    raise msg
230
231            neighbourdict[a,b] = (i, 2) #(id, edge)
232            neighbourdict[b,c] = (i, 0) #(id, edge)
233            neighbourdict[c,a] = (i, 1) #(id, edge)
234
235
236        #Step 2:
237        #Go through triangles again, but this time
238        #reverse direction of segments and lookup neighbours.
239        for i in range(N):
240            a = self.triangles[i, 0]
241            b = self.triangles[i, 1]
242            c = self.triangles[i, 2]
243
244            self.number_of_boundaries[i] = 3
245            if neighbourdict.has_key((b,a)):
246                self.neighbours[i, 2] = neighbourdict[b,a][0]
247                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
248                self.number_of_boundaries[i] -= 1
249
250            if neighbourdict.has_key((c,b)):
251                self.neighbours[i, 0] = neighbourdict[c,b][0]
252                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
253                self.number_of_boundaries[i] -= 1
254
255            if neighbourdict.has_key((a,c)):
256                self.neighbours[i, 1] = neighbourdict[a,c][0]
257                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
258                self.number_of_boundaries[i] -= 1
259
260
261    def build_surrogate_neighbour_structure(self):
262        """Build structure where each triangle edge points to its neighbours
263        if they exist. Otherwise point to the triangle itself.
264
265        The surrogate neighbour structure is useful for computing gradients
266        based on centroid values of neighbours.
267
268        Precondition: Neighbour structure is defined
269        Postcondition:
270          Surrogate neighbour structure is defined:
271          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
272          triangles.
273
274        """
275
276        N = self.number_of_elements
277        for i in range(N):
278            #Find all neighbouring volumes that are not boundaries
279            for k in range(3):
280                if self.neighbours[i, k] < 0:
281                    self.surrogate_neighbours[i, k] = i #Point this triangle
282                else:
283                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
284
285
286
287    def build_boundary_dictionary(self, boundary = None):
288        """Build or check the dictionary of boundary tags.
289         self.boundary is a dictionary of tags,
290         keyed by volume id and edge:
291         { (id, edge): tag, ... }
292
293         Postconditions:
294            self.boundary is defined.
295        """
296
297        from config import default_boundary_tag
298
299        if boundary is None:
300            boundary = {}
301            for vol_id in range(self.number_of_elements):
302                for edge_id in range(0, 3):
303                    if self.neighbours[vol_id, edge_id] < 0:
304                        boundary[(vol_id, edge_id)] = default_boundary_tag
305        else:
306            #Check that all keys in given boundary exist
307            for vol_id, edge_id in boundary.keys():
308                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
309                a, b = self.neighbours.shape
310                assert vol_id < a and edge_id < b, msg
311
312                #FIXME: This assert violates internal boundaries (delete it)
313                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
314                #assert self.neighbours[vol_id, edge_id] < 0, msg
315
316            #Check that all boundary segments are assigned a tag
317            for vol_id in range(self.number_of_elements):
318                for edge_id in range(0, 3):
319                    if self.neighbours[vol_id, edge_id] < 0:
320                        if not boundary.has_key( (vol_id, edge_id) ):
321                            msg = 'WARNING: Given boundary does not contain '
322                            msg += 'tags for edge (%d, %d). '\
323                                   %(vol_id, edge_id)
324                            msg += 'Assigning default tag (%s).'\
325                                   %default_boundary_tag
326
327                            #FIXME: Print only as per verbosity
328                            #print msg
329
330                            #FIXME: Make this situation an error in the future
331                            #and make another function which will
332                            #enable default boundary-tags where
333                            #tags a not specified
334                            boundary[ (vol_id, edge_id) ] =\
335                                      default_boundary_tag
336
337
338
339        self.boundary = boundary
340
341
342    def build_tagged_elements_dictionary(self, tagged_elements = None):
343        """Build the dictionary of element tags.
344         self.tagged_elements is a dictionary of element arrays,
345         keyed by tag:
346         { (tag): [e1, e2, e3..] }
347
348         Postconditions:
349            self.element_tag is defined
350        """
351        from Numeric import array, Int
352
353        if tagged_elements is None:
354            tagged_elements = {}
355        else:
356            #Check that all keys in given boundary exist
357            for tag in tagged_elements.keys():
358                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
359
360                msg = 'Not all elements exist. '
361                assert max(tagged_elements[tag]) < self.number_of_elements, msg
362        #print "tagged_elements", tagged_elements
363        self.tagged_elements = tagged_elements
364
365    def build_boundary_structure(self):
366        """Traverse boundary and
367        enumerate neighbour indices from -1 and
368        counting down.
369
370        Precondition:
371            self.boundary is defined.
372        Post condition:
373            neighbour array has unique negative indices for boundary
374            boundary_segments array imposes an ordering on segments
375            (not otherwise available from the dictionary)
376
377        Note: If a segment is listed in the boundary dictionary
378        it *will* become a boundary - even if there is a neighbouring triangle.
379        This would be the case for internal boundaries
380        """
381
382        #FIXME: Now Obsolete - maybe use some comments from here in
383        #domain.set_boundary
384
385        if self.boundary is None:
386            msg = 'Boundary dictionary must be defined before '
387            msg += 'building boundary structure'
388            raise msg
389
390
391        self.boundary_segments = self.boundary.keys()
392        self.boundary_segments.sort()
393
394        index = -1
395        for id, edge in self.boundary_segments:
396
397            #FIXME: One would detect internal boundaries as follows
398            #if self.neighbours[id, edge] > -1:
399            #    print 'Internal boundary'
400
401            self.neighbours[id, edge] = index
402            index -= 1
403
404
405    def get_boundary_tags(self):
406        """Return list of available boundary tags
407        """
408
409        tags = {}
410        for v in self.boundary.values():
411            tags[v] = 1
412
413        return tags.keys()
414
415
416    def get_boundary_polygon(self):
417        """Return bounding polygon as a list of points
418
419        FIXME: If triangles are listed as discontinuous
420        (e.g vertex coordinates listed multiple times),
421        this may not work as expected.
422        """
423        from Numeric import allclose, sqrt, array, minimum, maximum
424
425
426
427        #V = self.get_vertex_coordinates()
428        segments = {}
429
430        #pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
431        #pmax = (max(self.coordinates[:,0]), max(self.coordinates[:,1]))
432
433        #FIXME:Can this be written more compactly, e.g.
434        #using minimum and maximium?
435        pmin = array( [min(self.coordinates[:,0]),
436                       min(self.coordinates[:,1]) ] )
437
438        pmax = array( [max(self.coordinates[:,0]),
439                       max(self.coordinates[:,1]) ] )
440
441        mindist = sqrt(sum( (pmax-pmin)**2 ))
442        for i, edge_id in self.boundary.keys():
443            #Find vertex ids for boundary segment
444            if edge_id == 0: a = 1; b = 2
445            if edge_id == 1: a = 2; b = 0
446            if edge_id == 2: a = 0; b = 1
447
448            A = tuple(self.get_vertex_coordinate(i, a))
449            B = tuple(self.get_vertex_coordinate(i, b))
450
451            #Take the point closest to pmin as starting point
452            #Note: Could be arbitrary, but nice to have
453            #a unique way of selecting
454            dist_A = sqrt(sum( (A-pmin)**2 ))
455            dist_B = sqrt(sum( (B-pmin)**2 ))
456
457            #Find minimal point
458            if dist_A < mindist:
459                mindist = dist_A
460                p0 = A
461            if dist_B < mindist:
462                mindist = dist_B
463                p0 = B
464
465
466            if p0 is None:
467                raise 'Weird'
468                p0 = A #We need a starting point (FIXME)
469
470                print 'A', A
471                print 'B', B
472                print 'pmin', pmin
473                print
474
475            segments[A] = B
476
477
478        #Start with smallest point and follow boundary (counter clock wise)
479        polygon = [p0]
480        while len(polygon) < len(self.boundary):
481            p1 = segments[p0]
482            polygon.append(p1)
483            p0 = p1
484
485        return polygon
486
487
488    def check_integrity(self):
489        """Check that triangles are internally consistent e.g.
490        that area corresponds to edgelengths, that vertices
491        are arranged in a counter-clockwise order, etc etc
492        Neighbour structure will be checked by class Mesh
493        """
494
495        from config import epsilon
496        from math import pi
497        from util import anglediff
498
499        N = self.number_of_elements
500
501        #Get x,y coordinates for all vertices for all triangles
502        V = self.get_vertex_coordinates()
503
504        #Check each triangle
505        for i in range(N):
506            x0 = V[i, 0]; y0 = V[i, 1]
507            x1 = V[i, 2]; y1 = V[i, 3]
508            x2 = V[i, 4]; y2 = V[i, 5]
509
510            #Check that area hasn't been compromised
511            area = self.areas[i]
512            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
513            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
514                  %(x0,y0,x1,y1,x2,y2)
515            assert abs((area - ref)/area) < epsilon, msg
516
517            #Check that points are arranged in counter clock-wise order
518            v0 = [x1-x0, y1-y0]
519            v1 = [x2-x1, y2-y1]
520            v2 = [x0-x2, y0-y2]
521
522            a0 = anglediff(v1, v0)
523            a1 = anglediff(v2, v1)
524            a2 = anglediff(v0, v2)
525
526            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
527            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
528            assert a0 < pi and a1 < pi and a2 < pi, msg
529
530            #Check that normals are orthogonal to edge vectors
531            #Note that normal[k] lies opposite vertex k
532
533            normal0 = self.normals[i, 0:2]
534            normal1 = self.normals[i, 2:4]
535            normal2 = self.normals[i, 4:6]
536
537            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
538                assert u[0]*v[0] + u[1]*v[1] < epsilon
539
540
541        #Check that all vertices have been registered
542        for v_id, v in enumerate(self.vertexlist):
543
544            msg = 'Vertex %s does not belong to an element.'
545            #assert v is not None, msg
546            if v is None:
547                print msg%v_id
548
549        #Check integrity of neighbour structure
550        for i in range(N):
551            for v in self.triangles[i, :]:
552                #Check that all vertices have been registered
553                assert self.vertexlist[v] is not None
554
555                #Check that this triangle is listed with at least one vertex
556                assert (i, 0) in self.vertexlist[v] or\
557                       (i, 1) in self.vertexlist[v] or\
558                       (i, 2) in self.vertexlist[v]
559
560
561
562            #Check neighbour structure
563            for k, neighbour_id in enumerate(self.neighbours[i,:]):
564
565                #Assert that my neighbour's neighbour is me
566                #Boundaries need not fulfill this
567                if neighbour_id >= 0:
568                    edge = self.neighbour_edges[i, k]
569                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
570                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
571                    assert self.neighbours[neighbour_id, edge] == i ,msg
572
573
574
575        #Check that all boundaries have
576        # unique, consecutive, negative indices
577
578        #L = len(self.boundary)
579        #for i in range(L):
580        #    id, edge = self.boundary_segments[i]
581        #    assert self.neighbours[id, edge] == -i-1
582
583
584        #NOTE: This assert doesn't hold true if there are internal boundaries
585        #FIXME: Look into this further.
586        #FIXME (Ole): In pyvolution mark 3 this is OK again
587        #NOTE: No longer works because neighbour structure is modified by
588        #      domain set_boundary.
589        #for id, edge in self.boundary:
590        #    assert self.neighbours[id,edge] < 0
591        #
592        #NOTE (Ole): I reckon this was resolved late 2004?
593        #
594        #See domain.set_boundary
595
596
597    def get_centroid_coordinates(self):
598        """Return all centroid coordinates.
599        Return all centroid coordinates for all triangles as an Nx2 array
600        (ordered as x0, y0 for each triangle)
601        """
602        return self.centroid_coordinates
Note: See TracBrowser for help on using the repository browser.