source: inundation/pyvolution/mesh.py @ 2709

Last change on this file since 2709 was 2709, checked in by ole, 18 years ago

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

File size: 25.9 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
8from math import pi, sqrt
9       
10
11class Mesh(General_mesh):
12    """Collection of triangular elements (purely geometric)
13
14    A triangular element is defined in terms of three vertex ids,
15    ordered counter clock-wise,
16    each corresponding to a given coordinate set.
17    Vertices from different elements can point to the same
18    coordinate set.
19
20    Coordinate sets are implemented as an N x 2 Numeric array containing
21    x and y coordinates.
22
23
24    To instantiate:
25       Mesh(coordinates, triangles)
26
27    where
28
29      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
30      floats representing all x, y coordinates in the mesh.
31
32      triangles is either a list of 3-tuples or an Nx3 Numeric array of
33      integers representing indices of all vertices in the mesh.
34      Each vertex is identified by its index i in [0, M-1].
35
36
37    Example:
38        a = [0.0, 0.0]
39        b = [0.0, 2.0]
40        c = [2.0,0.0]
41        e = [2.0, 2.0]
42
43        points = [a, b, c, e]
44        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
45        mesh = Mesh(points, triangles)
46
47        #creates two triangles: bac and bce
48
49
50    Mesh takes the optional third argument boundary which is a
51    dictionary mapping from (element_id, edge_id) to boundary tag.
52    The default value is None which will assign the default_boundary_tag
53    as specified in config.py to all boundary edges.
54    """
55
56    #FIXME: Maybe rename coordinates to points (as in a poly file)
57    #But keep 'vertex_coordinates'
58
59    #FIXME: Put in check for angles less than a set minimum
60
61
62    def __init__(self, coordinates, triangles, boundary = None,
63                 tagged_elements = None, geo_reference = None,
64                 use_inscribed_circle = False):
65        """
66        Build triangles from x,y coordinates (sequence of 2-tuples or
67        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
68        or Nx3 Numeric array of non-negative integers).
69        """
70
71
72
73        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
74
75        General_mesh.__init__(self, coordinates, triangles, geo_reference)
76
77        N = self.number_of_elements
78
79        self.use_inscribed_circle = use_inscribed_circle
80
81        #Allocate space for geometric quantities
82        self.centroid_coordinates = zeros((N, 2), Float)
83
84        self.radii = zeros(N, Float)
85
86        self.neighbours = zeros((N, 3), Int)
87        self.neighbour_edges = zeros((N, 3), Int)
88        self.number_of_boundaries = zeros(N, Int)
89        self.surrogate_neighbours = zeros((N, 3), Int)
90
91        #Get x,y coordinates for all triangles and store
92        V = self.vertex_coordinates
93
94        #Initialise each triangle
95        for i in range(N):
96            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
97
98            x0 = V[i, 0]; y0 = V[i, 1]
99            x1 = V[i, 2]; y1 = V[i, 3]
100            x2 = V[i, 4]; y2 = V[i, 5]
101
102            #Compute centroid
103            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
104            self.centroid_coordinates[i] = centroid
105
106
107            if self.use_inscribed_circle == False:
108                #OLD code. Computed radii may exceed that of an
109                #inscribed circle
110
111                #Midpoints
112                m0 = array([(x1 + x2)/2, (y1 + y2)/2])
113                m1 = array([(x0 + x2)/2, (y0 + y2)/2])
114                m2 = array([(x1 + x0)/2, (y1 + y0)/2])
115
116                #The radius is the distance from the centroid of
117                #a triangle to the midpoint of the side of the triangle
118                #closest to the centroid
119                d0 = sqrt(sum( (centroid-m0)**2 ))
120                d1 = sqrt(sum( (centroid-m1)**2 ))
121                d2 = sqrt(sum( (centroid-m2)**2 ))
122
123                self.radii[i] = min(d0, d1, d2)
124
125            else:
126                #NEW code added by Peter Row. True radius
127                #of inscribed circle is computed
128
129                a = sqrt((x0-x1)**2+(y0-y1)**2)
130                b = sqrt((x1-x2)**2+(y1-y2)**2)
131                c = sqrt((x2-x0)**2+(y2-y0)**2)
132
133                self.radii[i]=2.0*self.areas[i]/(a+b+c)
134
135
136            #Initialise Neighbours (-1 means that it is a boundary neighbour)
137            self.neighbours[i, :] = [-1, -1, -1]
138
139            #Initialise edge ids of neighbours
140            #In case of boundaries this slot is not used
141            self.neighbour_edges[i, :] = [-1, -1, -1]
142
143
144        #Build neighbour structure
145        self.build_neighbour_structure()
146
147        #Build surrogate neighbour structure
148        self.build_surrogate_neighbour_structure()
149
150        #Build boundary dictionary mapping (id, edge) to symbolic tags
151        self.build_boundary_dictionary(boundary)
152
153        #Build tagged element  dictionary mapping (tag) to array of elements
154        self.build_tagged_elements_dictionary(tagged_elements)
155
156        #Update boundary indices FIXME: OBSOLETE
157        #self.build_boundary_structure()
158
159        #FIXME check integrity?
160
161    def __repr__(self):
162        return 'Mesh: %d triangles, %d elements, %d boundary segments'\
163               %(self.coordinates.shape[0], len(self), len(self.boundary))
164
165
166    def set_to_inscribed_circle(self,safety_factor = 1):
167        #FIXME phase out eventually
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, verbose = False):
417        """Return bounding polygon for mesh (counter clockwise)
418
419        Using the mesh boundary, derive a bounding polygon for this mesh.
420        If multiple vertex values are present, the algorithm will select the
421        path that contains the mesh.
422        """
423       
424        from Numeric import allclose, sqrt, array, minimum, maximum
425        from utilities.numerical_tools import angle, ensure_numeric       
426
427
428        # Get mesh extent
429        xmin, xmax, ymin, ymax = self.get_extent()
430        pmin = ensure_numeric([xmin, ymin])
431        pmax = ensure_numeric([xmax, ymax])       
432
433
434        # Assemble dictionary of boundary segments and choose starting point
435        segments = {}
436        inverse_segments = {}
437        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
438        for i, edge_id in self.boundary.keys():
439            # Find vertex ids for boundary segment
440            if edge_id == 0: a = 1; b = 2
441            if edge_id == 1: a = 2; b = 0
442            if edge_id == 2: a = 0; b = 1
443
444            A = self.get_vertex_coordinate(i, a) # Start
445            B = self.get_vertex_coordinate(i, b) # End
446
447            # Take the point closest to pmin as starting point
448            # Note: Could be arbitrary, but nice to have
449            # a unique way of selecting
450            dist_A = sqrt(sum((A-pmin)**2))
451            dist_B = sqrt(sum((B-pmin)**2))
452
453            #Find lower leftmost point
454            if dist_A < mindist:
455                mindist = dist_A
456                p0 = A
457            if dist_B < mindist:
458                mindist = dist_B
459                p0 = B
460
461
462            # Sanity check
463            if p0 is None:
464                raise Exception('Impossible')
465
466
467            # Register potential paths from A to B
468            if not segments.has_key(tuple(A)):
469                segments[tuple(A)] = [] # Empty list for candidate points               
470               
471            segments[tuple(A)].append(B)               
472
473
474
475
476        #Start with smallest point and follow boundary (counter clock wise)
477        polygon = [p0]      # Storage for final boundary polygon
478        point_registry = {} # Keep track of storage to avoid multiple runs around boundary
479                            # This will only be the case if there are more than one candidate
480                            # FIXME (Ole): Perhaps we can do away with polygon and use
481                            # only point_registry to save space.
482
483        point_registry[tuple(p0)] = 0                           
484                           
485        #while len(polygon) < len(self.boundary):
486        while len(point_registry) < len(self.boundary):
487
488            candidate_list = segments[tuple(p0)]
489            if len(candidate_list) > 1:
490                # Multiple points detected
491                # Take the candidate that is furthest to the clockwise direction,
492                # as that will follow the boundary.
493
494
495                if verbose:
496                    print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list)
497
498
499                # Choose vector against which all angles will be measured
500                if len(polygon) > 1:   
501                    v_prev = p0 - polygon[-2] # Vector that leads to p0
502                else:
503                    # FIXME (Ole): What do we do if the first point has multiple candidates?
504                    # Being the lower left corner, perhaps we can use the
505                    # vector [1, 0], but I really don't know if this is completely
506                    # watertight.
507                    # Another option might be v_prev = [1.0, 0.0]
508                    v_prev = [1.0, 0.0]
509                   
510
511                   
512                # Choose candidate with minimum angle   
513                minimum_angle = 2*pi
514                for pc in candidate_list:
515                    vc = pc-p0  # Candidate vector
516                   
517                    # Angle between each candidate and the previous vector in [-pi, pi]
518                    ac = angle(vc, v_prev)
519                    if ac > pi: ac = pi-ac
520
521                    # take the minimal angle corresponding to the rightmost vector
522                    if ac < minimum_angle:
523                        minimum_angle = ac
524                        p1 = pc             # Best candidate
525                       
526
527                if verbose is True:
528                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
529               
530            else:
531                p1 = candidate_list[0]
532
533            if point_registry.has_key(tuple(p1)):
534                # We have completed the boundary polygon - yeehaa
535                break 
536            else:
537                point_registry[tuple(p1)] = len(point_registry)
538           
539            polygon.append(p1)
540            p0 = p1
541
542
543        return polygon
544
545
546    def check_integrity(self):
547        """Check that triangles are internally consistent e.g.
548        that area corresponds to edgelengths, that vertices
549        are arranged in a counter-clockwise order, etc etc
550        Neighbour structure will be checked by class Mesh
551        """
552
553        from config import epsilon
554        from utilities.numerical_tools import anglediff
555
556        N = self.number_of_elements
557        #Get x,y coordinates for all vertices for all triangles
558        V = self.get_vertex_coordinates()
559        #Check each triangle
560        for i in range(N):
561           
562            x0 = V[i, 0]; y0 = V[i, 1]
563            x1 = V[i, 2]; y1 = V[i, 3]
564            x2 = V[i, 4]; y2 = V[i, 5]
565
566            #Check that area hasn't been compromised
567            area = self.areas[i]
568            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
569            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
570                  %(x0,y0,x1,y1,x2,y2)
571            assert abs((area - ref)/area) < epsilon, msg
572
573            #Check that points are arranged in counter clock-wise order
574            v0 = [x1-x0, y1-y0]
575            v1 = [x2-x1, y2-y1]
576            v2 = [x0-x2, y0-y2]
577
578            a0 = anglediff(v1, v0)
579            a1 = anglediff(v2, v1)
580            a2 = anglediff(v0, v2)
581
582            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
583            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
584            assert a0 < pi and a1 < pi and a2 < pi, msg
585
586            #Check that normals are orthogonal to edge vectors
587            #Note that normal[k] lies opposite vertex k
588
589            normal0 = self.normals[i, 0:2]
590            normal1 = self.normals[i, 2:4]
591            normal2 = self.normals[i, 4:6]
592
593            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
594
595                #Normalise
596                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
597                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
598
599                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v #Inner product
600               
601                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
602                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
603                msg += ' Inner product is %e.' %x               
604                assert x < epsilon, msg
605
606
607        #Check that all vertices have been registered
608        for v_id, v in enumerate(self.vertexlist):
609
610            msg = 'Vertex %s does not belong to an element.'
611            #assert v is not None, msg
612            if v is None:
613                print msg%v_id
614
615        #Check integrity of neighbour structure
616        for i in range(N):
617#            print i
618            for v in self.triangles[i, :]:
619                #Check that all vertices have been registered
620                assert self.vertexlist[v] is not None
621
622                #Check that this triangle is listed with at least one vertex
623                assert (i, 0) in self.vertexlist[v] or\
624                       (i, 1) in self.vertexlist[v] or\
625                       (i, 2) in self.vertexlist[v]
626
627
628
629            #Check neighbour structure
630            for k, neighbour_id in enumerate(self.neighbours[i,:]):
631
632                #Assert that my neighbour's neighbour is me
633                #Boundaries need not fulfill this
634                if neighbour_id >= 0:
635                    edge = self.neighbour_edges[i, k]
636                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
637                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
638                    assert self.neighbours[neighbour_id, edge] == i ,msg
639
640
641
642        #Check that all boundaries have
643        # unique, consecutive, negative indices
644
645        #L = len(self.boundary)
646        #for i in range(L):
647        #    id, edge = self.boundary_segments[i]
648        #    assert self.neighbours[id, edge] == -i-1
649
650
651        #NOTE: This assert doesn't hold true if there are internal boundaries
652        #FIXME: Look into this further.
653        #FIXME (Ole): In pyvolution mark 3 this is OK again
654        #NOTE: No longer works because neighbour structure is modified by
655        #      domain set_boundary.
656        #for id, edge in self.boundary:
657        #    assert self.neighbours[id,edge] < 0
658        #
659        #NOTE (Ole): I reckon this was resolved late 2004?
660        #
661        #See domain.set_boundary
662
663
664    def get_centroid_coordinates(self):
665        """Return all centroid coordinates.
666        Return all centroid coordinates for all triangles as an Nx2 array
667        (ordered as x0, y0 for each triangle)
668        """
669        return self.centroid_coordinates
670
671
672
673
674    def statistics(self):
675        """Output statistics about mesh
676        """
677
678        from Numeric import arange
679        from utilities.numerical_tools import histogram
680
681        vertex_coordinates = self.vertex_coordinates
682        areas = self.areas
683        x = vertex_coordinates[:,0]
684        y = vertex_coordinates[:,1]
685
686
687        #Setup 10 bins for area histogram
688        m = max(areas)
689        bins = arange(0., m, m/10)
690        hist = histogram(areas, bins)
691
692        str =  '------------------------------------------------\n'
693        str += 'Mesh statistics:\n'
694        str += '  Number of triangles = %d\n' %self.number_of_elements
695        str += '  Extent:\n'
696        str += '    x in [%f, %f]\n' %(min(x), max(x))
697        str += '    y in [%f, %f]\n' %(min(y), max(y))
698        str += '  Areas:\n'
699        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
700        str += '    number of distinct areas: %d\n' %(len(areas))       
701        str += '    Histogram:\n'
702
703        hi = bins[0]
704        for i, count in enumerate(hist):
705            lo = hi
706            if i+1 < len(bins):
707                #Open upper interval               
708                hi = bins[i+1]
709                str += '      [%f, %f[: %d\n' %(lo, hi, count)               
710            else:
711                #Closed upper interval
712                hi = m
713                str += '      [%f, %f]: %d\n' %(lo, hi, count)
714
715        N = len(areas)
716        if N > 10:
717            str += '    Percentiles (10%):\n'
718            areas = areas.tolist()
719            areas.sort()
720
721            k = 0
722            lower = min(areas)
723            for i, a in enumerate(areas):       
724                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
725                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
726                    lower = a
727                    k = i
728                   
729            str += '      %d triangles in [%f, %f]\n'\
730                   %(N-k, lower, max(areas))                   
731               
732                     
733        str += 'Boundary:\n'
734        str += '  Number of boundary segments == %d\n' %(len(self.boundary))
735        str += '  Boundary tags == %s\n' %self.get_boundary_tags() 
736        str += '------------------------------------------------\n'
737       
738
739        return str
740
Note: See TracBrowser for help on using the repository browser.