source: anuga_core/source/anuga/pyvolution/neighbour_mesh.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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