source: inundation/pyvolution/mesh.py @ 3072

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

Made boundary polygon return absolute UTM coordinates as per ticket:81

File size: 27.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
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 # Relative 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        All points are in absolute UTM coordinates
436        """
437       
438        from Numeric import allclose, sqrt, array, minimum, maximum
439        from utilities.numerical_tools import angle, ensure_numeric       
440
441
442        # Get mesh extent
443        xmin, xmax, ymin, ymax = self.get_extent(absolute=True)
444        pmin = ensure_numeric([xmin, ymin])
445        pmax = ensure_numeric([xmax, ymax])       
446
447
448        # Assemble dictionary of boundary segments and choose starting point
449        segments = {}
450        inverse_segments = {}
451        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
452        for i, edge_id in self.boundary.keys():
453            # Find vertex ids for boundary segment
454            if edge_id == 0: a = 1; b = 2
455            if edge_id == 1: a = 2; b = 0
456            if edge_id == 2: a = 0; b = 1
457
458            A = self.get_vertex_coordinate(i, a, absolute=True) # Start
459            B = self.get_vertex_coordinate(i, b, absolute=True) # End
460
461
462            # Take the point closest to pmin as starting point
463            # Note: Could be arbitrary, but nice to have
464            # a unique way of selecting
465            dist_A = sqrt(sum((A-pmin)**2))
466            dist_B = sqrt(sum((B-pmin)**2))
467
468            #Find lower leftmost point
469            if dist_A < mindist:
470                mindist = dist_A
471                p0 = A
472            if dist_B < mindist:
473                mindist = dist_B
474                p0 = B
475
476
477            # Sanity check
478            if p0 is None:
479                raise Exception('Impossible')
480
481
482            # Register potential paths from A to B
483            if not segments.has_key(tuple(A)):
484                segments[tuple(A)] = [] # Empty list for candidate points               
485               
486            segments[tuple(A)].append(B)               
487
488
489
490           
491        #Start with smallest point and follow boundary (counter clock wise)
492        polygon = [p0]      # Storage for final boundary polygon
493        point_registry = {} # Keep track of storage to avoid multiple runs around
494                            # boundary. This will only be the case if there are
495                            # more than one candidate.
496                            # FIXME (Ole): Perhaps we can do away with polygon
497                            # and use only point_registry to save space.
498
499        point_registry[tuple(p0)] = 0                           
500                           
501        #while len(polygon) < len(self.boundary):
502        while len(point_registry) < len(self.boundary):
503
504            candidate_list = segments[tuple(p0)]
505            if len(candidate_list) > 1:
506                # Multiple points detected
507                # Take the candidate that is furthest to the clockwise direction,
508                # as that will follow the boundary.
509
510
511                if verbose:
512                    print 'Point %s has multiple candidates: %s'\
513                          %(str(p0), candidate_list)
514
515
516                # Choose vector against which all angles will be measured
517                if len(polygon) > 1:   
518                    v_prev = p0 - polygon[-2] # Vector that leads to p0
519                else:
520                    # FIXME (Ole): What do we do if the first point has multiple
521                    # candidates?
522                    # Being the lower left corner, perhaps we can use the
523                    # vector [1, 0], but I really don't know if this is completely
524                    # watertight.
525                    # Another option might be v_prev = [1.0, 0.0]
526                    v_prev = [1.0, 0.0]
527                   
528
529                   
530                # Choose candidate with minimum angle   
531                minimum_angle = 2*pi
532                for pc in candidate_list:
533                    vc = pc-p0  # Candidate vector
534                   
535                    # Angle between each candidate and the previous vector
536                    # in [-pi, pi]
537                    ac = angle(vc, v_prev)
538                    if ac > pi: ac = pi-ac
539
540                    # take the minimal angle corresponding to the rightmost vector
541                    if ac < minimum_angle:
542                        minimum_angle = ac
543                        p1 = pc             # Best candidate
544                       
545
546                if verbose is True:
547                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
548               
549            else:
550                p1 = candidate_list[0]
551
552            if point_registry.has_key(tuple(p1)):
553                # We have completed the boundary polygon - yeehaa
554                break 
555            else:
556                point_registry[tuple(p1)] = len(point_registry)
557           
558            polygon.append(p1)
559            p0 = p1
560
561
562        return polygon
563
564
565    def check_integrity(self):
566        """Check that triangles are internally consistent e.g.
567        that area corresponds to edgelengths, that vertices
568        are arranged in a counter-clockwise order, etc etc
569        Neighbour structure will be checked by class Mesh
570        """
571
572        from config import epsilon
573        from utilities.numerical_tools import anglediff
574
575        N = self.number_of_elements
576        #Get x,y coordinates for all vertices for all triangles
577        V = self.get_vertex_coordinates()
578        #Check each triangle
579        for i in range(N):
580           
581            x0 = V[i, 0]; y0 = V[i, 1]
582            x1 = V[i, 2]; y1 = V[i, 3]
583            x2 = V[i, 4]; y2 = V[i, 5]
584
585            #Check that area hasn't been compromised
586            area = self.areas[i]
587            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
588            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
589                  %(x0,y0,x1,y1,x2,y2)
590            assert abs((area - ref)/area) < epsilon, msg
591
592            #Check that points are arranged in counter clock-wise order
593            v0 = [x1-x0, y1-y0]
594            v1 = [x2-x1, y2-y1]
595            v2 = [x0-x2, y0-y2]
596
597            a0 = anglediff(v1, v0)
598            a1 = anglediff(v2, v1)
599            a2 = anglediff(v0, v2)
600
601            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
602            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
603            assert a0 < pi and a1 < pi and a2 < pi, msg
604
605            #Check that normals are orthogonal to edge vectors
606            #Note that normal[k] lies opposite vertex k
607
608            normal0 = self.normals[i, 0:2]
609            normal1 = self.normals[i, 2:4]
610            normal2 = self.normals[i, 4:6]
611
612            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
613
614                #Normalise
615                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
616                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
617
618                x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v #Inner product
619               
620                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
621                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
622                msg += ' Inner product is %e.' %x               
623                assert x < epsilon, msg
624
625        self.lone_vertices = []
626        #Check that all vertices have been registered
627        for v_id, v in enumerate(self.vertexlist):
628
629            #msg = 'Vertex %s does not belong to an element.'
630            #assert v is not None, msg
631            if v is None:
632                #print msg%v_id
633                self.lone_vertices.append(v_id)
634
635        #Check integrity of neighbour structure
636        for i in range(N):
637#            print i
638            for v in self.triangles[i, :]:
639                #Check that all vertices have been registered
640                assert self.vertexlist[v] is not None
641
642                #Check that this triangle is listed with at least one vertex
643                assert (i, 0) in self.vertexlist[v] or\
644                       (i, 1) in self.vertexlist[v] or\
645                       (i, 2) in self.vertexlist[v]
646
647
648
649            #Check neighbour structure
650            for k, neighbour_id in enumerate(self.neighbours[i,:]):
651
652                #Assert that my neighbour's neighbour is me
653                #Boundaries need not fulfill this
654                if neighbour_id >= 0:
655                    edge = self.neighbour_edges[i, k]
656                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
657                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
658                    assert self.neighbours[neighbour_id, edge] == i ,msg
659
660
661
662        #Check that all boundaries have
663        # unique, consecutive, negative indices
664
665        #L = len(self.boundary)
666        #for i in range(L):
667        #    id, edge = self.boundary_segments[i]
668        #    assert self.neighbours[id, edge] == -i-1
669
670
671        #NOTE: This assert doesn't hold true if there are internal boundaries
672        #FIXME: Look into this further.
673        #FIXME (Ole): In pyvolution mark 3 this is OK again
674        #NOTE: No longer works because neighbour structure is modified by
675        #      domain set_boundary.
676        #for id, edge in self.boundary:
677        #    assert self.neighbours[id,edge] < 0
678        #
679        #NOTE (Ole): I reckon this was resolved late 2004?
680        #
681        #See domain.set_boundary
682
683    def get_lone_vertices(self):
684        """Return a list of vertices that are not connected to any triangles.
685
686        Precondition
687        FIXME(DSG - DSG) Pull the code out of check integrity that builds this
688                         structure.
689        check_integrity has to have been called.
690        """
691        return self.lone_vertices
692
693    def get_centroid_coordinates(self):
694        """Return all centroid coordinates.
695        Return all centroid coordinates for all triangles as an Nx2 array
696        (ordered as x0, y0 for each triangle)
697        """
698        return self.centroid_coordinates
699
700
701
702
703    def statistics(self):
704        """Output statistics about mesh
705        """
706
707        from Numeric import arange
708        from utilities.numerical_tools import histogram, create_bins
709
710        vertex_coordinates = self.vertex_coordinates # Relative coordinates
711        areas = self.areas
712        x = vertex_coordinates[:,0]
713        y = vertex_coordinates[:,1]
714
715
716        #Setup 10 bins for area histogram
717        bins = create_bins(areas, 10)
718        #m = max(areas)
719        #bins = arange(0., m, m/10)
720        hist = histogram(areas, bins)
721
722        str =  '------------------------------------------------\n'
723        str += 'Mesh statistics:\n'
724        str += '  Number of triangles = %d\n' %self.number_of_elements
725        str += '  Extent:\n'
726        str += '    x in [%f, %f]\n' %(min(x), max(x))
727        str += '    y in [%f, %f]\n' %(min(y), max(y))
728        str += '  Areas:\n'
729        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
730        str += '    number of distinct areas: %d\n' %(len(areas))       
731        str += '    Histogram:\n'
732
733        hi = bins[0]
734        for i, count in enumerate(hist):
735            lo = hi
736            if i+1 < len(bins):
737                #Open upper interval               
738                hi = bins[i+1]
739                str += '      [%f, %f[: %d\n' %(lo, hi, count)               
740            else:
741                #Closed upper interval
742                hi = max(areas)
743                str += '      [%f, %f]: %d\n' %(lo, hi, count)
744
745        N = len(areas)
746        if N > 10:
747            str += '    Percentiles (10%):\n'
748            areas = areas.tolist()
749            areas.sort()
750
751            k = 0
752            lower = min(areas)
753            for i, a in enumerate(areas):       
754                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
755                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
756                    lower = a
757                    k = i
758                   
759            str += '      %d triangles in [%f, %f]\n'\
760                   %(N-k, lower, max(areas))                   
761               
762                     
763        str += 'Boundary:\n'
764        str += '  Number of boundary segments == %d\n' %(len(self.boundary))
765        str += '  Boundary tags == %s\n' %self.get_boundary_tags() 
766        str += '------------------------------------------------\n'
767       
768
769        return str
770
Note: See TracBrowser for help on using the repository browser.