# source:anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py@4592

Last change on this file since 4592 was 4592, checked in by ole, 16 years ago

Allowed set_quantity restriction by polygon to observe geo_referencing.

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