source: inundation/pyvolution/mesh.py @ 3062

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

Added FIXMES about ticket:81

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