source: inundation/pyvolution/mesh.py @ 2689

Last change on this file since 2689 was 2683, checked in by nick, 18 years ago

removed print statements from least_squares.py and mesh.py

File size: 23.4 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 utilities.numerical_tools import anglediff
498
499        N = self.number_of_elements
500        #Get x,y coordinates for all vertices for all triangles
501        V = self.get_vertex_coordinates()
502        #Check each triangle
503        for i in range(N):
504           
505            x0 = V[i, 0]; y0 = V[i, 1]
506            x1 = V[i, 2]; y1 = V[i, 3]
507            x2 = V[i, 4]; y2 = V[i, 5]
508
509            #Check that area hasn't been compromised
510            area = self.areas[i]
511            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
512            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
513                  %(x0,y0,x1,y1,x2,y2)
514            assert abs((area - ref)/area) < epsilon, msg
515
516            #Check that points are arranged in counter clock-wise order
517            v0 = [x1-x0, y1-y0]
518            v1 = [x2-x1, y2-y1]
519            v2 = [x0-x2, y0-y2]
520
521            a0 = anglediff(v1, v0)
522            a1 = anglediff(v2, v1)
523            a2 = anglediff(v0, v2)
524
525            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
526            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
527            assert a0 < pi and a1 < pi and a2 < pi, msg
528
529            #Check that normals are orthogonal to edge vectors
530            #Note that normal[k] lies opposite vertex k
531
532            normal0 = self.normals[i, 0:2]
533            normal1 = self.normals[i, 2:4]
534            normal2 = self.normals[i, 4:6]
535
536            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
537
538                #Normalise
539                from math import sqrt
540                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
541                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
542
543                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v #Inner product
544               
545                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
546                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
547                msg += ' Inner product is %e.' %x               
548                assert x < epsilon, msg
549
550
551        #Check that all vertices have been registered
552        for v_id, v in enumerate(self.vertexlist):
553
554            msg = 'Vertex %s does not belong to an element.'
555            #assert v is not None, msg
556            if v is None:
557                print msg%v_id
558
559        #Check integrity of neighbour structure
560        for i in range(N):
561#            print i
562            for v in self.triangles[i, :]:
563                #Check that all vertices have been registered
564                assert self.vertexlist[v] is not None
565
566                #Check that this triangle is listed with at least one vertex
567                assert (i, 0) in self.vertexlist[v] or\
568                       (i, 1) in self.vertexlist[v] or\
569                       (i, 2) in self.vertexlist[v]
570
571
572
573            #Check neighbour structure
574            for k, neighbour_id in enumerate(self.neighbours[i,:]):
575
576                #Assert that my neighbour's neighbour is me
577                #Boundaries need not fulfill this
578                if neighbour_id >= 0:
579                    edge = self.neighbour_edges[i, k]
580                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
581                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
582                    assert self.neighbours[neighbour_id, edge] == i ,msg
583
584
585
586        #Check that all boundaries have
587        # unique, consecutive, negative indices
588
589        #L = len(self.boundary)
590        #for i in range(L):
591        #    id, edge = self.boundary_segments[i]
592        #    assert self.neighbours[id, edge] == -i-1
593
594
595        #NOTE: This assert doesn't hold true if there are internal boundaries
596        #FIXME: Look into this further.
597        #FIXME (Ole): In pyvolution mark 3 this is OK again
598        #NOTE: No longer works because neighbour structure is modified by
599        #      domain set_boundary.
600        #for id, edge in self.boundary:
601        #    assert self.neighbours[id,edge] < 0
602        #
603        #NOTE (Ole): I reckon this was resolved late 2004?
604        #
605        #See domain.set_boundary
606
607
608    def get_centroid_coordinates(self):
609        """Return all centroid coordinates.
610        Return all centroid coordinates for all triangles as an Nx2 array
611        (ordered as x0, y0 for each triangle)
612        """
613        return self.centroid_coordinates
614
615
616
617
618    def statistics(self):
619        """Output statistics about mesh
620        """
621
622        from Numeric import arange
623        from utilities.numerical_tools import histogram
624
625        vertex_coordinates = self.vertex_coordinates
626        areas = self.areas
627        x = vertex_coordinates[:,0]
628        y = vertex_coordinates[:,1]
629
630
631        #Setup 10 bins for area histogram
632        m = max(areas)
633        bins = arange(0., m, m/10)
634        hist = histogram(areas, bins)
635
636        str =  '------------------------------------------------\n'
637        str += 'Mesh statistics:\n'
638        str += '  Number of triangles = %d\n' %self.number_of_elements
639        str += '  Extent:\n'
640        str += '    x in [%f, %f]\n' %(min(x), max(x))
641        str += '    y in [%f, %f]\n' %(min(y), max(y))
642        str += '  Areas:\n'
643        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
644        str += '    number of distinct areas: %d\n' %(len(areas))       
645        str += '    Histogram:\n'
646
647        hi = bins[0]
648        for i, count in enumerate(hist):
649            lo = hi
650            if i+1 < len(bins):
651                #Open upper interval               
652                hi = bins[i+1]
653                str += '      [%f, %f[: %d\n' %(lo, hi, count)               
654            else:
655                #Closed upper interval
656                hi = m
657                str += '      [%f, %f]: %d\n' %(lo, hi, count)
658
659        N = len(areas)
660        if N > 10:
661            str += '    Percentiles (10%):\n'
662            areas = areas.tolist()
663            areas.sort()
664
665            k = 0
666            lower = min(areas)
667            for i, a in enumerate(areas):       
668                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
669                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
670                    lower = a
671                    k = i
672                   
673            str += '      %d triangles in [%f, %f]\n'\
674                   %(N-k, lower, max(areas))                   
675               
676                     
677        str += 'Boundary:\n'
678        str += '  Number of boundary segments == %d\n' %(len(self.boundary))
679        str += '  Boundary tags == %s\n' %self.get_boundary_tags() 
680        str += '------------------------------------------------\n'
681       
682
683        return str
684
Note: See TracBrowser for help on using the repository browser.