source: inundation/pyvolution/mesh.py @ 2982

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

bugger!

File size: 26.5 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,
63                 boundary=None,
64                 tagged_elements=None,
65                 geo_reference=None,
66                 use_inscribed_circle=False,
67                 verbose=False):
68        """
69        Build triangles from x,y coordinates (sequence of 2-tuples or
70        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
71        or Nx3 Numeric array of non-negative integers).
72        """
73
74
75
76        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
77
78        General_mesh.__init__(self, coordinates, triangles,
79                              geo_reference, verbose=verbose)
80
81        if verbose: print 'Initialising mesh'         
82
83        N = self.number_of_elements
84
85        self.use_inscribed_circle = use_inscribed_circle
86
87        #Allocate space for geometric quantities
88        self.centroid_coordinates = zeros((N, 2), Float)
89
90        self.radii = zeros(N, Float)
91
92        self.neighbours = zeros((N, 3), Int)
93        self.neighbour_edges = zeros((N, 3), Int)
94        self.number_of_boundaries = zeros(N, Int)
95        self.surrogate_neighbours = zeros((N, 3), Int)
96
97        #Get x,y coordinates for all triangles and store
98        V = self.vertex_coordinates
99
100        #Initialise each triangle
101        if verbose: print 'Mesh: Computing centroids and radii'       
102        for i in range(N):
103            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
104
105            x0 = V[i, 0]; y0 = V[i, 1]
106            x1 = V[i, 2]; y1 = V[i, 3]
107            x2 = V[i, 4]; y2 = V[i, 5]
108
109            #Compute centroid
110            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
111            self.centroid_coordinates[i] = centroid
112
113
114            if self.use_inscribed_circle == False:
115                #OLD code. Computed radii may exceed that of an
116                #inscribed circle
117
118                #Midpoints
119                m0 = array([(x1 + x2)/2, (y1 + y2)/2])
120                m1 = array([(x0 + x2)/2, (y0 + y2)/2])
121                m2 = array([(x1 + x0)/2, (y1 + y0)/2])
122
123                #The radius is the distance from the centroid of
124                #a triangle to the midpoint of the side of the triangle
125                #closest to the centroid
126                d0 = sqrt(sum( (centroid-m0)**2 ))
127                d1 = sqrt(sum( (centroid-m1)**2 ))
128                d2 = sqrt(sum( (centroid-m2)**2 ))
129
130                self.radii[i] = min(d0, d1, d2)
131
132            else:
133                #NEW code added by Peter Row. True radius
134                #of inscribed circle is computed
135
136                a = sqrt((x0-x1)**2+(y0-y1)**2)
137                b = sqrt((x1-x2)**2+(y1-y2)**2)
138                c = sqrt((x2-x0)**2+(y2-y0)**2)
139
140                self.radii[i]=2.0*self.areas[i]/(a+b+c)
141
142
143            #Initialise Neighbours (-1 means that it is a boundary neighbour)
144            self.neighbours[i, :] = [-1, -1, -1]
145
146            #Initialise edge ids of neighbours
147            #In case of boundaries this slot is not used
148            self.neighbour_edges[i, :] = [-1, -1, -1]
149
150
151        #Build neighbour structure
152        if verbose: print 'Mesh: Building neigbour structure'               
153        self.build_neighbour_structure()
154
155        #Build surrogate neighbour structure
156        if verbose: print 'Mesh: Building surrogate neigbour structure'
157        self.build_surrogate_neighbour_structure()
158
159        #Build boundary dictionary mapping (id, edge) to symbolic tags
160        if verbose: print 'Mesh: Building boundary dictionary'
161        self.build_boundary_dictionary(boundary)
162
163        #Build tagged element  dictionary mapping (tag) to array of elements
164        if verbose: print 'Mesh: Building tagged elements dictionary'       
165        self.build_tagged_elements_dictionary(tagged_elements)
166
167        #Update boundary indices FIXME: OBSOLETE
168        #self.build_boundary_structure()
169
170        #FIXME check integrity?
171        if verbose: print 'Mesh: Done'               
172
173    def __repr__(self):
174        return 'Mesh: %d triangles, %d elements, %d boundary segments'\
175               %(self.coordinates.shape[0], len(self), len(self.boundary))
176
177
178    def set_to_inscribed_circle(self,safety_factor = 1):
179        #FIXME phase out eventually
180        N = self.number_of_elements
181        V = self.vertex_coordinates
182
183        #initialising min and max ratio
184        i=0
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        ratio = old_rad/self.radii[i]
193        max_ratio = ratio
194        min_ratio = ratio
195
196        for i in range(N):
197            old_rad = self.radii[i]
198            x0 = V[i, 0]; y0 = V[i, 1]
199            x1 = V[i, 2]; y1 = V[i, 3]
200            x2 = V[i, 4]; y2 = V[i, 5]
201            a = sqrt((x0-x1)**2+(y0-y1)**2)
202            b = sqrt((x1-x2)**2+(y1-y2)**2)
203            c = sqrt((x2-x0)**2+(y2-y0)**2)
204            self.radii[i]=self.areas[i]/(2*(a+b+c))*safety_factor
205            ratio = old_rad/self.radii[i]
206            if ratio >= max_ratio: max_ratio = ratio
207            if ratio <= min_ratio: min_ratio = ratio
208        return max_ratio,min_ratio
209
210    def build_neighbour_structure(self):
211        """Update all registered triangles to point to their neighbours.
212
213        Also, keep a tally of the number of boundaries for each triangle
214
215        Postconditions:
216          neighbours and neighbour_edges is populated
217          number_of_boundaries integer array is defined.
218        """
219
220        #Step 1:
221        #Build dictionary mapping from segments (2-tuple of points)
222        #to left hand side edge (facing neighbouring triangle)
223
224        N = self.number_of_elements
225        neighbourdict = {}
226        for i in range(N):
227
228            #Register all segments as keys mapping to current triangle
229            #and segment id
230            a = self.triangles[i, 0]
231            b = self.triangles[i, 1]
232            c = self.triangles[i, 2]
233            if neighbourdict.has_key((a,b)):
234                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
235                    raise msg
236            if neighbourdict.has_key((b,c)):
237                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
238                    raise msg
239            if neighbourdict.has_key((c,a)):
240                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
241                    raise msg
242
243            neighbourdict[a,b] = (i, 2) #(id, edge)
244            neighbourdict[b,c] = (i, 0) #(id, edge)
245            neighbourdict[c,a] = (i, 1) #(id, edge)
246
247
248        #Step 2:
249        #Go through triangles again, but this time
250        #reverse direction of segments and lookup neighbours.
251        for i in range(N):
252            a = self.triangles[i, 0]
253            b = self.triangles[i, 1]
254            c = self.triangles[i, 2]
255
256            self.number_of_boundaries[i] = 3
257            if neighbourdict.has_key((b,a)):
258                self.neighbours[i, 2] = neighbourdict[b,a][0]
259                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
260                self.number_of_boundaries[i] -= 1
261
262            if neighbourdict.has_key((c,b)):
263                self.neighbours[i, 0] = neighbourdict[c,b][0]
264                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
265                self.number_of_boundaries[i] -= 1
266
267            if neighbourdict.has_key((a,c)):
268                self.neighbours[i, 1] = neighbourdict[a,c][0]
269                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
270                self.number_of_boundaries[i] -= 1
271
272
273    def build_surrogate_neighbour_structure(self):
274        """Build structure where each triangle edge points to its neighbours
275        if they exist. Otherwise point to the triangle itself.
276
277        The surrogate neighbour structure is useful for computing gradients
278        based on centroid values of neighbours.
279
280        Precondition: Neighbour structure is defined
281        Postcondition:
282          Surrogate neighbour structure is defined:
283          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
284          triangles.
285
286        """
287
288        N = self.number_of_elements
289        for i in range(N):
290            #Find all neighbouring volumes that are not boundaries
291            for k in range(3):
292                if self.neighbours[i, k] < 0:
293                    self.surrogate_neighbours[i, k] = i #Point this triangle
294                else:
295                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
296
297
298
299    def build_boundary_dictionary(self, boundary = None):
300        """Build or check the dictionary of boundary tags.
301         self.boundary is a dictionary of tags,
302         keyed by volume id and edge:
303         { (id, edge): tag, ... }
304
305         Postconditions:
306            self.boundary is defined.
307        """
308
309        from config import default_boundary_tag
310
311        if boundary is None:
312            boundary = {}
313            for vol_id in range(self.number_of_elements):
314                for edge_id in range(0, 3):
315                    if self.neighbours[vol_id, edge_id] < 0:
316                        boundary[(vol_id, edge_id)] = default_boundary_tag
317        else:
318            #Check that all keys in given boundary exist
319            for vol_id, edge_id in boundary.keys():
320                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
321                a, b = self.neighbours.shape
322                assert vol_id < a and edge_id < b, msg
323
324                #FIXME: This assert violates internal boundaries (delete it)
325                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
326                #assert self.neighbours[vol_id, edge_id] < 0, msg
327
328            #Check that all boundary segments are assigned a tag
329            for vol_id in range(self.number_of_elements):
330                for edge_id in range(0, 3):
331                    if self.neighbours[vol_id, edge_id] < 0:
332                        if not boundary.has_key( (vol_id, edge_id) ):
333                            msg = 'WARNING: Given boundary does not contain '
334                            msg += 'tags for edge (%d, %d). '\
335                                   %(vol_id, edge_id)
336                            msg += 'Assigning default tag (%s).'\
337                                   %default_boundary_tag
338
339                            #FIXME: Print only as per verbosity
340                            #print msg
341
342                            #FIXME: Make this situation an error in the future
343                            #and make another function which will
344                            #enable default boundary-tags where
345                            #tags a not specified
346                            boundary[ (vol_id, edge_id) ] =\
347                                      default_boundary_tag
348
349
350
351        self.boundary = boundary
352
353
354    def build_tagged_elements_dictionary(self, tagged_elements = None):
355        """Build the dictionary of element tags.
356         self.tagged_elements is a dictionary of element arrays,
357         keyed by tag:
358         { (tag): [e1, e2, e3..] }
359
360         Postconditions:
361            self.element_tag is defined
362        """
363        from Numeric import array, Int
364
365        if tagged_elements is None:
366            tagged_elements = {}
367        else:
368            #Check that all keys in given boundary exist
369            for tag in tagged_elements.keys():
370                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
371
372                msg = 'Not all elements exist. '
373                assert max(tagged_elements[tag]) < self.number_of_elements, msg
374        #print "tagged_elements", tagged_elements
375        self.tagged_elements = tagged_elements
376
377    def build_boundary_structure(self):
378        """Traverse boundary and
379        enumerate neighbour indices from -1 and
380        counting down.
381
382        Precondition:
383            self.boundary is defined.
384        Post condition:
385            neighbour array has unique negative indices for boundary
386            boundary_segments array imposes an ordering on segments
387            (not otherwise available from the dictionary)
388
389        Note: If a segment is listed in the boundary dictionary
390        it *will* become a boundary - even if there is a neighbouring triangle.
391        This would be the case for internal boundaries
392        """
393
394        #FIXME: Now Obsolete - maybe use some comments from here in
395        #domain.set_boundary
396
397        if self.boundary is None:
398            msg = 'Boundary dictionary must be defined before '
399            msg += 'building boundary structure'
400            raise msg
401
402
403        self.boundary_segments = self.boundary.keys()
404        self.boundary_segments.sort()
405
406        index = -1
407        for id, edge in self.boundary_segments:
408
409            #FIXME: One would detect internal boundaries as follows
410            #if self.neighbours[id, edge] > -1:
411            #    print 'Internal boundary'
412
413            self.neighbours[id, edge] = index
414            index -= 1
415
416
417    def get_boundary_tags(self):
418        """Return list of available boundary tags
419        """
420
421        tags = {}
422        for v in self.boundary.values():
423            tags[v] = 1
424
425        return tags.keys()
426
427
428    def get_boundary_polygon(self, verbose = False):
429        """Return bounding polygon for mesh (counter clockwise)
430
431        Using the mesh boundary, derive a bounding polygon for this mesh.
432        If multiple vertex values are present, the algorithm will select the
433        path that contains the mesh.
434        """
435       
436        from Numeric import allclose, sqrt, array, minimum, maximum
437        from utilities.numerical_tools import angle, ensure_numeric       
438
439
440        # Get mesh extent
441        xmin, xmax, ymin, ymax = self.get_extent()
442        pmin = ensure_numeric([xmin, ymin])
443        pmax = ensure_numeric([xmax, ymax])       
444
445
446        # Assemble dictionary of boundary segments and choose starting point
447        segments = {}
448        inverse_segments = {}
449        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
450        for i, edge_id in self.boundary.keys():
451            # Find vertex ids for boundary segment
452            if edge_id == 0: a = 1; b = 2
453            if edge_id == 1: a = 2; b = 0
454            if edge_id == 2: a = 0; b = 1
455
456            A = self.get_vertex_coordinate(i, a) # Start
457            B = self.get_vertex_coordinate(i, b) # End
458
459            # Take the point closest to pmin as starting point
460            # Note: Could be arbitrary, but nice to have
461            # a unique way of selecting
462            dist_A = sqrt(sum((A-pmin)**2))
463            dist_B = sqrt(sum((B-pmin)**2))
464
465            #Find lower leftmost point
466            if dist_A < mindist:
467                mindist = dist_A
468                p0 = A
469            if dist_B < mindist:
470                mindist = dist_B
471                p0 = B
472
473
474            # Sanity check
475            if p0 is None:
476                raise Exception('Impossible')
477
478
479            # Register potential paths from A to B
480            if not segments.has_key(tuple(A)):
481                segments[tuple(A)] = [] # Empty list for candidate points               
482               
483            segments[tuple(A)].append(B)               
484
485
486
487
488        #Start with smallest point and follow boundary (counter clock wise)
489        polygon = [p0]      # Storage for final boundary polygon
490        point_registry = {} # Keep track of storage to avoid multiple runs around boundary
491                            # This will only be the case if there are more than one candidate
492                            # FIXME (Ole): Perhaps we can do away with polygon and use
493                            # only point_registry to save space.
494
495        point_registry[tuple(p0)] = 0                           
496                           
497        #while len(polygon) < len(self.boundary):
498        while len(point_registry) < len(self.boundary):
499
500            candidate_list = segments[tuple(p0)]
501            if len(candidate_list) > 1:
502                # Multiple points detected
503                # Take the candidate that is furthest to the clockwise direction,
504                # as that will follow the boundary.
505
506
507                if verbose:
508                    print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list)
509
510
511                # Choose vector against which all angles will be measured
512                if len(polygon) > 1:   
513                    v_prev = p0 - polygon[-2] # Vector that leads to p0
514                else:
515                    # FIXME (Ole): What do we do if the first point has multiple candidates?
516                    # Being the lower left corner, perhaps we can use the
517                    # vector [1, 0], but I really don't know if this is completely
518                    # watertight.
519                    # Another option might be v_prev = [1.0, 0.0]
520                    v_prev = [1.0, 0.0]
521                   
522
523                   
524                # Choose candidate with minimum angle   
525                minimum_angle = 2*pi
526                for pc in candidate_list:
527                    vc = pc-p0  # Candidate vector
528                   
529                    # Angle between each candidate and the previous vector in [-pi, pi]
530                    ac = angle(vc, v_prev)
531                    if ac > pi: ac = pi-ac
532
533                    # take the minimal angle corresponding to the rightmost vector
534                    if ac < minimum_angle:
535                        minimum_angle = ac
536                        p1 = pc             # Best candidate
537                       
538
539                if verbose is True:
540                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
541               
542            else:
543                p1 = candidate_list[0]
544
545            if point_registry.has_key(tuple(p1)):
546                # We have completed the boundary polygon - yeehaa
547                break 
548            else:
549                point_registry[tuple(p1)] = len(point_registry)
550           
551            polygon.append(p1)
552            p0 = p1
553
554
555        return polygon
556
557
558    def check_integrity(self):
559        """Check that triangles are internally consistent e.g.
560        that area corresponds to edgelengths, that vertices
561        are arranged in a counter-clockwise order, etc etc
562        Neighbour structure will be checked by class Mesh
563        """
564
565        from config import epsilon
566        from utilities.numerical_tools import anglediff
567
568        N = self.number_of_elements
569        #Get x,y coordinates for all vertices for all triangles
570        V = self.get_vertex_coordinates()
571        #Check each triangle
572        for i in range(N):
573           
574            x0 = V[i, 0]; y0 = V[i, 1]
575            x1 = V[i, 2]; y1 = V[i, 3]
576            x2 = V[i, 4]; y2 = V[i, 5]
577
578            #Check that area hasn't been compromised
579            area = self.areas[i]
580            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
581            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
582                  %(x0,y0,x1,y1,x2,y2)
583            assert abs((area - ref)/area) < epsilon, msg
584
585            #Check that points are arranged in counter clock-wise order
586            v0 = [x1-x0, y1-y0]
587            v1 = [x2-x1, y2-y1]
588            v2 = [x0-x2, y0-y2]
589
590            a0 = anglediff(v1, v0)
591            a1 = anglediff(v2, v1)
592            a2 = anglediff(v0, v2)
593
594            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
595            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
596            assert a0 < pi and a1 < pi and a2 < pi, msg
597
598            #Check that normals are orthogonal to edge vectors
599            #Note that normal[k] lies opposite vertex k
600
601            normal0 = self.normals[i, 0:2]
602            normal1 = self.normals[i, 2:4]
603            normal2 = self.normals[i, 4:6]
604
605            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
606
607                #Normalise
608                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
609                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
610
611                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v #Inner product
612               
613                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
614                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
615                msg += ' Inner product is %e.' %x               
616                assert x < epsilon, msg
617
618
619        #Check that all vertices have been registered
620        for v_id, v in enumerate(self.vertexlist):
621
622            msg = 'Vertex %s does not belong to an element.'
623            #assert v is not None, msg
624            if v is None:
625                print msg%v_id
626
627        #Check integrity of neighbour structure
628        for i in range(N):
629#            print i
630            for v in self.triangles[i, :]:
631                #Check that all vertices have been registered
632                assert self.vertexlist[v] is not None
633
634                #Check that this triangle is listed with at least one vertex
635                assert (i, 0) in self.vertexlist[v] or\
636                       (i, 1) in self.vertexlist[v] or\
637                       (i, 2) in self.vertexlist[v]
638
639
640
641            #Check neighbour structure
642            for k, neighbour_id in enumerate(self.neighbours[i,:]):
643
644                #Assert that my neighbour's neighbour is me
645                #Boundaries need not fulfill this
646                if neighbour_id >= 0:
647                    edge = self.neighbour_edges[i, k]
648                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
649                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
650                    assert self.neighbours[neighbour_id, edge] == i ,msg
651
652
653
654        #Check that all boundaries have
655        # unique, consecutive, negative indices
656
657        #L = len(self.boundary)
658        #for i in range(L):
659        #    id, edge = self.boundary_segments[i]
660        #    assert self.neighbours[id, edge] == -i-1
661
662
663        #NOTE: This assert doesn't hold true if there are internal boundaries
664        #FIXME: Look into this further.
665        #FIXME (Ole): In pyvolution mark 3 this is OK again
666        #NOTE: No longer works because neighbour structure is modified by
667        #      domain set_boundary.
668        #for id, edge in self.boundary:
669        #    assert self.neighbours[id,edge] < 0
670        #
671        #NOTE (Ole): I reckon this was resolved late 2004?
672        #
673        #See domain.set_boundary
674
675
676    def get_centroid_coordinates(self):
677        """Return all centroid coordinates.
678        Return all centroid coordinates for all triangles as an Nx2 array
679        (ordered as x0, y0 for each triangle)
680        """
681        return self.centroid_coordinates
682
683
684
685
686    def statistics(self):
687        """Output statistics about mesh
688        """
689
690        from Numeric import arange
691        from utilities.numerical_tools import histogram, create_bins
692
693        vertex_coordinates = self.vertex_coordinates
694        areas = self.areas
695        x = vertex_coordinates[:,0]
696        y = vertex_coordinates[:,1]
697
698
699        #Setup 10 bins for area histogram
700        bins = create_bins(areas, 10)
701        #m = max(areas)
702        #bins = arange(0., m, m/10)
703        hist = histogram(areas, bins)
704
705        str =  '------------------------------------------------\n'
706        str += 'Mesh statistics:\n'
707        str += '  Number of triangles = %d\n' %self.number_of_elements
708        str += '  Extent:\n'
709        str += '    x in [%f, %f]\n' %(min(x), max(x))
710        str += '    y in [%f, %f]\n' %(min(y), max(y))
711        str += '  Areas:\n'
712        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
713        str += '    number of distinct areas: %d\n' %(len(areas))       
714        str += '    Histogram:\n'
715
716        hi = bins[0]
717        for i, count in enumerate(hist):
718            lo = hi
719            if i+1 < len(bins):
720                #Open upper interval               
721                hi = bins[i+1]
722                str += '      [%f, %f[: %d\n' %(lo, hi, count)               
723            else:
724                #Closed upper interval
725                hi = max(areas)
726                str += '      [%f, %f]: %d\n' %(lo, hi, count)
727
728        N = len(areas)
729        if N > 10:
730            str += '    Percentiles (10%):\n'
731            areas = areas.tolist()
732            areas.sort()
733
734            k = 0
735            lower = min(areas)
736            for i, a in enumerate(areas):       
737                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
738                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
739                    lower = a
740                    k = i
741                   
742            str += '      %d triangles in [%f, %f]\n'\
743                   %(N-k, lower, max(areas))                   
744               
745                     
746        str += 'Boundary:\n'
747        str += '  Number of boundary segments == %d\n' %(len(self.boundary))
748        str += '  Boundary tags == %s\n' %self.get_boundary_tags() 
749        str += '------------------------------------------------\n'
750       
751
752        return str
753
Note: See TracBrowser for help on using the repository browser.