source: anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py @ 5521

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

Investigated performance profile of get_flow_through_cross_sections.
Most time appear to be spent in the interpolation routines.

File size: 39.3 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
12from Numeric import array, allclose
13       
14
15class Mesh(General_mesh):
16    """Collection of triangular elements (purely geometric)
17
18    A triangular element is defined in terms of three vertex ids,
19    ordered counter clock-wise,
20    each corresponding to a given coordinate set.
21    Vertices from different elements can point to the same
22    coordinate set.
23
24    Coordinate sets are implemented as an N x 2 Numeric array containing
25    x and y coordinates.
26
27
28    To instantiate:
29       Mesh(coordinates, triangles)
30
31    where
32
33      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
34      floats representing all x, y coordinates in the mesh.
35
36      triangles is either a list of 3-tuples or an Nx3 Numeric array of
37      integers representing indices of all vertices in the mesh.
38      Each vertex is identified by its index i in [0, M-1].
39
40
41    Example:
42        a = [0.0, 0.0]
43        b = [0.0, 2.0]
44        c = [2.0,0.0]
45        e = [2.0, 2.0]
46
47        points = [a, b, c, e]
48        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
49        mesh = Mesh(points, triangles)
50
51        #creates two triangles: bac and bce
52
53
54    Mesh takes the optional third argument boundary which is a
55    dictionary mapping from (element_id, edge_id) to boundary tag.
56    The default value is None which will assign the default_boundary_tag
57    as specified in config.py to all boundary edges.
58    """
59
60    #FIXME: Maybe rename coordinates to points (as in a poly file)
61    #But keep 'vertex_coordinates'
62
63    #FIXME: Put in check for angles less than a set minimum
64
65
66    def __init__(self, coordinates, triangles,
67                 boundary=None,
68                 tagged_elements=None,
69                 geo_reference=None,
70                 number_of_full_nodes=None,
71                 number_of_full_triangles=None,
72                 use_inscribed_circle=False,
73                 verbose=False):
74        """
75        Build triangles from x,y coordinates (sequence of 2-tuples or
76        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
77        or Nx3 Numeric array of non-negative integers).
78        """
79
80
81
82        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
83
84        General_mesh.__init__(self, coordinates, triangles,
85                              number_of_full_nodes=\
86                              number_of_full_nodes,
87                              number_of_full_triangles=\
88                              number_of_full_triangles,
89                              geo_reference=geo_reference,
90                              verbose=verbose)
91
92        if verbose: print 'Initialising mesh'         
93
94        N = len(self) #Number_of_triangles
95
96        self.use_inscribed_circle = use_inscribed_circle
97
98        #Allocate space for geometric quantities
99        self.centroid_coordinates = zeros((N, 2), Float)
100
101        self.radii = zeros(N, Float)
102
103        self.neighbours = zeros((N, 3), Int)
104        self.neighbour_edges = zeros((N, 3), Int)
105        self.number_of_boundaries = zeros(N, Int)
106        self.surrogate_neighbours = zeros((N, 3), Int)
107
108        #Get x,y coordinates for all triangles and store
109        V = self.vertex_coordinates # Relative coordinates
110
111        #Initialise each triangle
112        if verbose: print 'Mesh: Computing centroids and radii'       
113        for i in range(N):
114            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
115
116            x0, y0 = V[3*i, :]
117            x1, y1 = V[3*i+1, :]
118            x2, y2 = V[3*i+2, :]                       
119
120            #x0 = V[i, 0]; y0 = V[i, 1]
121            #x1 = V[i, 2]; y1 = V[i, 3]
122            #x2 = V[i, 4]; y2 = V[i, 5]
123
124            #Compute centroid
125            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
126            self.centroid_coordinates[i] = centroid
127
128
129            if self.use_inscribed_circle == False:
130                #OLD code. Computed radii may exceed that of an
131                #inscribed circle
132
133                #Midpoints
134                m0 = array([(x1 + x2)/2, (y1 + y2)/2])
135                m1 = array([(x0 + x2)/2, (y0 + y2)/2])
136                m2 = array([(x1 + x0)/2, (y1 + y0)/2])
137
138                #The radius is the distance from the centroid of
139                #a triangle to the midpoint of the side of the triangle
140                #closest to the centroid
141                d0 = sqrt(sum( (centroid-m0)**2 ))
142                d1 = sqrt(sum( (centroid-m1)**2 ))
143                d2 = sqrt(sum( (centroid-m2)**2 ))
144
145                self.radii[i] = min(d0, d1, d2)
146
147            else:
148                #NEW code added by Peter Row. True radius
149                #of inscribed circle is computed
150
151                a = sqrt((x0-x1)**2+(y0-y1)**2)
152                b = sqrt((x1-x2)**2+(y1-y2)**2)
153                c = sqrt((x2-x0)**2+(y2-y0)**2)
154
155                self.radii[i]=2.0*self.areas[i]/(a+b+c)
156
157
158            #Initialise Neighbours (-1 means that it is a boundary neighbour)
159            self.neighbours[i, :] = [-1, -1, -1]
160
161            #Initialise edge ids of neighbours
162            #In case of boundaries this slot is not used
163            self.neighbour_edges[i, :] = [-1, -1, -1]
164
165
166        #Build neighbour structure
167        if verbose: print 'Mesh: Building neigbour structure'               
168        self.build_neighbour_structure()
169
170        #Build surrogate neighbour structure
171        if verbose: print 'Mesh: Building surrogate neigbour structure'
172        self.build_surrogate_neighbour_structure()
173
174        #Build boundary dictionary mapping (id, edge) to symbolic tags
175        if verbose: print 'Mesh: Building boundary dictionary'
176        self.build_boundary_dictionary(boundary)
177
178        #Build tagged element  dictionary mapping (tag) to array of elements
179        if verbose: print 'Mesh: Building tagged elements dictionary'       
180        self.build_tagged_elements_dictionary(tagged_elements)
181
182        # Build a list of vertices that are not connected to any triangles
183        self.lone_vertices = []
184        #Check that all vertices have been registered
185        for node, count in enumerate(self.number_of_triangles_per_node):   
186            #msg = 'Node %d does not belong to an element.' %node
187            #assert count > 0, msg
188            if count == 0:
189                self.lone_vertices.append(node)
190               
191        #Update boundary indices FIXME: OBSOLETE
192        #self.build_boundary_structure()
193
194        #FIXME check integrity?
195        if verbose: print 'Mesh: Done'               
196
197    def __repr__(self):
198        return General_mesh.__repr__(self) + ', %d boundary segments'\
199               %(len(self.boundary))
200
201
202    def set_to_inscribed_circle(self,safety_factor = 1):
203        #FIXME phase out eventually
204        N = self.number_of_triangles
205        V = self.vertex_coordinates
206
207        #initialising min and max ratio
208        i=0
209        old_rad = self.radii[i]
210        x0 = V[i, 0]; y0 = V[i, 1]
211        x1 = V[i, 2]; y1 = V[i, 3]
212        x2 = V[i, 4]; y2 = V[i, 5]
213        a = sqrt((x0-x1)**2+(y0-y1)**2)
214        b = sqrt((x1-x2)**2+(y1-y2)**2)
215        c = sqrt((x2-x0)**2+(y2-y0)**2)
216        ratio = old_rad/self.radii[i]
217        max_ratio = ratio
218        min_ratio = ratio
219
220        for i in range(N):
221            old_rad = self.radii[i]
222            x0 = V[i, 0]; y0 = V[i, 1]
223            x1 = V[i, 2]; y1 = V[i, 3]
224            x2 = V[i, 4]; y2 = V[i, 5]
225            a = sqrt((x0-x1)**2+(y0-y1)**2)
226            b = sqrt((x1-x2)**2+(y1-y2)**2)
227            c = sqrt((x2-x0)**2+(y2-y0)**2)
228            self.radii[i]=self.areas[i]/(2*(a+b+c))*safety_factor
229            ratio = old_rad/self.radii[i]
230            if ratio >= max_ratio: max_ratio = ratio
231            if ratio <= min_ratio: min_ratio = ratio
232        return max_ratio,min_ratio
233
234
235
236    def build_neighbour_structure(self):
237        """Update all registered triangles to point to their neighbours.
238
239        Also, keep a tally of the number of boundaries for each triangle
240
241        Postconditions:
242          neighbours and neighbour_edges is populated
243          number_of_boundaries integer array is defined.
244        """
245
246        #Step 1:
247        #Build dictionary mapping from segments (2-tuple of points)
248        #to left hand side edge (facing neighbouring triangle)
249
250        N = len(self) #Number_of_triangles
251        neighbourdict = {}
252        for i in range(N):
253
254            #Register all segments as keys mapping to current triangle
255            #and segment id
256            a = self.triangles[i, 0]
257            b = self.triangles[i, 1]
258            c = self.triangles[i, 2]
259            if neighbourdict.has_key((a,b)):
260                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
261                    raise Exception, msg
262            if neighbourdict.has_key((b,c)):
263                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
264                    raise msg
265            if neighbourdict.has_key((c,a)):
266                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
267                    raise msg
268
269            neighbourdict[a,b] = (i, 2) #(id, edge)
270            neighbourdict[b,c] = (i, 0) #(id, edge)
271            neighbourdict[c,a] = (i, 1) #(id, edge)
272
273
274        #Step 2:
275        #Go through triangles again, but this time
276        #reverse direction of segments and lookup neighbours.
277        for i in range(N):
278            a = self.triangles[i, 0]
279            b = self.triangles[i, 1]
280            c = self.triangles[i, 2]
281
282            self.number_of_boundaries[i] = 3
283            if neighbourdict.has_key((b,a)):
284                self.neighbours[i, 2] = neighbourdict[b,a][0]
285                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
286                self.number_of_boundaries[i] -= 1
287
288            if neighbourdict.has_key((c,b)):
289                self.neighbours[i, 0] = neighbourdict[c,b][0]
290                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
291                self.number_of_boundaries[i] -= 1
292
293            if neighbourdict.has_key((a,c)):
294                self.neighbours[i, 1] = neighbourdict[a,c][0]
295                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
296                self.number_of_boundaries[i] -= 1
297
298
299    def build_surrogate_neighbour_structure(self):
300        """Build structure where each triangle edge points to its neighbours
301        if they exist. Otherwise point to the triangle itself.
302
303        The surrogate neighbour structure is useful for computing gradients
304        based on centroid values of neighbours.
305
306        Precondition: Neighbour structure is defined
307        Postcondition:
308          Surrogate neighbour structure is defined:
309          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
310          triangles.
311
312        """
313
314        N = len(self) #Number of triangles
315        for i in range(N):
316            #Find all neighbouring volumes that are not boundaries
317            for k in range(3):
318                if self.neighbours[i, k] < 0:
319                    self.surrogate_neighbours[i, k] = i #Point this triangle
320                else:
321                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
322
323
324
325    def build_boundary_dictionary(self, boundary = None):
326        """Build or check the dictionary of boundary tags.
327         self.boundary is a dictionary of tags,
328         keyed by volume id and edge:
329         { (id, edge): tag, ... }
330
331         Postconditions:
332            self.boundary is defined.
333        """
334
335        from anuga.config import default_boundary_tag
336
337        if boundary is None:
338            boundary = {}
339            for vol_id in range(len(self)):
340                for edge_id in range(0, 3):
341                    if self.neighbours[vol_id, edge_id] < 0:
342                        boundary[(vol_id, edge_id)] = default_boundary_tag
343        else:
344            #Check that all keys in given boundary exist
345            for vol_id, edge_id in boundary.keys():
346                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
347                a, b = self.neighbours.shape
348                assert vol_id < a and edge_id < b, msg
349
350                #FIXME: This assert violates internal boundaries (delete it)
351                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
352                #assert self.neighbours[vol_id, edge_id] < 0, msg
353
354            #Check that all boundary segments are assigned a tag
355            for vol_id in range(len(self)):
356                for edge_id in range(0, 3):
357                    if self.neighbours[vol_id, edge_id] < 0:
358                        if not boundary.has_key( (vol_id, edge_id) ):
359                            msg = 'WARNING: Given boundary does not contain '
360                            msg += 'tags for edge (%d, %d). '\
361                                   %(vol_id, edge_id)
362                            msg += 'Assigning default tag (%s).'\
363                                   %default_boundary_tag
364
365                            #FIXME: Print only as per verbosity
366                            #print msg
367
368                            #FIXME: Make this situation an error in the future
369                            #and make another function which will
370                            #enable default boundary-tags where
371                            #tags a not specified
372                            boundary[ (vol_id, edge_id) ] =\
373                                      default_boundary_tag
374
375
376
377        self.boundary = boundary
378
379
380    def build_tagged_elements_dictionary(self, tagged_elements = None):
381        """Build the dictionary of element tags.
382         self.tagged_elements is a dictionary of element arrays,
383         keyed by tag:
384         { (tag): [e1, e2, e3..] }
385
386         Postconditions:
387            self.element_tag is defined
388        """
389        from Numeric import array, Int
390
391        if tagged_elements is None:
392            tagged_elements = {}
393        else:
394            #Check that all keys in given boundary exist
395            for tag in tagged_elements.keys():
396                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
397
398                msg = 'Not all elements exist. '
399                assert max(tagged_elements[tag]) < len(self), msg
400        #print "tagged_elements", tagged_elements
401        self.tagged_elements = tagged_elements
402
403    def build_boundary_structure(self):
404        """Traverse boundary and
405        enumerate neighbour indices from -1 and
406        counting down.
407
408        Precondition:
409            self.boundary is defined.
410        Post condition:
411            neighbour array has unique negative indices for boundary
412            boundary_segments array imposes an ordering on segments
413            (not otherwise available from the dictionary)
414
415        Note: If a segment is listed in the boundary dictionary
416        it *will* become a boundary - even if there is a neighbouring triangle.
417        This would be the case for internal boundaries
418        """
419
420        #FIXME: Now Obsolete - maybe use some comments from here in
421        #domain.set_boundary
422
423        if self.boundary is None:
424            msg = 'Boundary dictionary must be defined before '
425            msg += 'building boundary structure'
426            raise msg
427
428
429        self.boundary_segments = self.boundary.keys()
430        self.boundary_segments.sort()
431
432        index = -1
433        for id, edge in self.boundary_segments:
434
435            #FIXME: One would detect internal boundaries as follows
436            #if self.neighbours[id, edge] > -1:
437            #    print 'Internal boundary'
438
439            self.neighbours[id, edge] = index
440            index -= 1
441
442
443    def get_boundary_tags(self):
444        """Return list of available boundary tags
445        """
446
447        tags = {}
448        for v in self.boundary.values():
449            tags[v] = 1
450
451        return tags.keys()
452
453
454    def get_boundary_polygon(self, verbose=False):
455        """Return bounding polygon for mesh (counter clockwise)
456
457        Using the mesh boundary, derive a bounding polygon for this mesh.
458        If multiple vertex values are present (vertices stored uniquely),
459        the algorithm will select the path that contains the entire mesh.
460
461        All points are in absolute UTM coordinates
462        """
463       
464        from Numeric import allclose, sqrt, array, minimum, maximum
465        from anuga.utilities.numerical_tools import angle, ensure_numeric     
466
467
468        # Get mesh extent
469        xmin, xmax, ymin, ymax = self.get_extent(absolute=True)
470        pmin = ensure_numeric([xmin, ymin])
471        pmax = ensure_numeric([xmax, ymax])       
472
473
474        # Assemble dictionary of boundary segments and choose starting point
475        segments = {}
476        inverse_segments = {}
477        p0 = None
478        mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh
479        for i, edge_id in self.boundary.keys():
480            # Find vertex ids for boundary segment
481            if edge_id == 0: a = 1; b = 2
482            if edge_id == 1: a = 2; b = 0
483            if edge_id == 2: a = 0; b = 1
484
485            A = self.get_vertex_coordinate(i, a, absolute=True) # Start
486            B = self.get_vertex_coordinate(i, b, absolute=True) # End
487
488
489            # Take the point closest to pmin as starting point
490            # Note: Could be arbitrary, but nice to have
491            # a unique way of selecting
492            dist_A = sqrt(sum((A-pmin)**2))
493            dist_B = sqrt(sum((B-pmin)**2))
494
495            # Find lower leftmost point
496            if dist_A < mindist:
497                mindist = dist_A
498                p0 = A
499            if dist_B < mindist:
500                mindist = dist_B
501                p0 = B
502
503
504            # Sanity check
505            if p0 is None:
506                raise Exception('Impossible')
507
508
509            # Register potential paths from A to B
510            if not segments.has_key(tuple(A)):
511                segments[tuple(A)] = [] # Empty list for candidate points
512
513            segments[tuple(A)].append(B)               
514
515
516        # Start with smallest point and follow boundary (counter clock wise)
517        polygon = [list(p0)]# Storage for final boundary polygon
518        point_registry = {} # Keep track of storage to avoid multiple runs
519                            # around boundary. This will only be the case if
520                            # there are more than one candidate.
521                            # FIXME (Ole): Perhaps we can do away with polygon
522                            # and use only point_registry to save space.
523
524        point_registry[tuple(p0)] = 0                           
525                           
526        while len(point_registry) < len(self.boundary):
527
528            candidate_list = segments[tuple(p0)]
529            if len(candidate_list) > 1:
530                # Multiple points detected (this will be the case for meshes
531                # with duplicate points as those used for discontinuous
532                # triangles with vertices stored uniquely).
533                # Take the candidate that is furthest to the clockwise
534                # direction, as that will follow the boundary.
535                #
536                # This will also be the case for pathological triangles
537                # that have no neighbours.
538
539                if verbose:
540                    print 'Point %s has multiple candidates: %s'\
541                          %(str(p0), candidate_list)
542
543                # Check that previous are not in candidate list
544                #for p in candidate_list:
545                #    assert not allclose(p0, p)
546
547                # Choose vector against which all angles will be measured
548                if len(polygon) > 1:   
549                    v_prev = p0 - polygon[-2] # Vector that leads to p0
550                                              # from previous point
551                else:
552                    # FIXME (Ole): What do we do if the first point has
553                    # multiple candidates?
554                    # Being the lower left corner, perhaps we can use the
555                    # vector [1, 0], but I really don't know if this is
556                    # completely watertight.
557                    v_prev = [1.0, 0.0]
558                   
559
560                # Choose candidate with minimum angle   
561                minimum_angle = 2*pi
562                for pc in candidate_list:
563
564                    vc = pc-p0  # Candidate vector (from p0 to candidate pt)
565                   
566                    # Angle between each candidate and the previous vector
567                    # in [-pi, pi]
568                    ac = angle(vc, v_prev)
569                    if ac > pi:
570                        # Give preference to angles on the right hand side
571                        # of v_prev
572                        # print 'pc = %s, changing angle from %f to %f'\
573                        # %(pc, ac*180/pi, (ac-2*pi)*180/pi)
574                        ac = ac-2*pi
575
576                    # Take the minimal angle corresponding to the
577                    # rightmost vector
578                    if ac < minimum_angle:
579                        minimum_angle = ac
580                        p1 = pc             # Best candidate
581                       
582
583                if verbose is True:
584                    print '  Best candidate %s, angle %f'\
585                          %(p1, minimum_angle*180/pi)
586                   
587            else:
588                p1 = candidate_list[0]
589
590               
591            if point_registry.has_key(tuple(p1)):
592                # We have reached a point already visited.
593               
594                if allclose(p1, polygon[0]):
595                    # If it is the initial point, the polygon is complete.
596                   
597                    if verbose is True:
598                        print '  Stop criterion fulfilled at point %s' %p1
599                        print polygon               
600                       
601                    # We have completed the boundary polygon - yeehaa
602                    break
603                else:   
604                    # The point already visited is not the initial point
605                    # This would be a pathological triangle, but the
606                    # algorithm must be able to deal with this
607                    pass
608   
609            else:
610                # We are still finding new points on the boundary
611                point_registry[tuple(p1)] = len(point_registry)
612           
613            polygon.append(list(p1)) # De-Numeric each point :-)
614            p0 = p1
615
616
617        return polygon
618
619
620    def check_integrity(self):
621        """Check that triangles are internally consistent e.g.
622        that area corresponds to edgelengths, that vertices
623        are arranged in a counter-clockwise order, etc etc
624        Neighbour structure will be checked by class Mesh
625        """
626
627        from anuga.config import epsilon
628        from anuga.utilities.numerical_tools import anglediff
629
630        from Numeric import sort, allclose
631
632        N = len(self)
633
634        # Get x,y coordinates for all vertices for all triangles
635        V = self.get_vertex_coordinates()
636
637        # Check each triangle
638        for i in range(N):
639
640            x0, y0 = V[3*i, :]
641            x1, y1 = V[3*i+1, :]
642            x2, y2 = V[3*i+2, :]
643           
644            # Check that area hasn't been compromised
645            area = self.areas[i]
646            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
647            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
648                  %(x0,y0,x1,y1,x2,y2)
649            assert abs((area - ref)/area) < epsilon, msg
650
651            # Check that points are arranged in counter clock-wise order
652            v0 = [x1-x0, y1-y0]
653            v1 = [x2-x1, y2-y1]
654            v2 = [x0-x2, y0-y2]
655            a0 = anglediff(v1, v0)
656            a1 = anglediff(v2, v1)
657            a2 = anglediff(v0, v2)
658
659            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
660            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
661            assert a0 < pi and a1 < pi and a2 < pi, msg
662
663            # Check that normals are orthogonal to edge vectors
664            # Note that normal[k] lies opposite vertex k
665
666            normal0 = self.normals[i, 0:2]
667            normal1 = self.normals[i, 2:4]
668            normal2 = self.normals[i, 4:6]
669
670            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
671
672                # Normalise
673                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
674                l_v = sqrt(v[0]*v[0] + v[1]*v[1])
675
676                msg = 'Normal vector in triangle %d does not have unit length' %i
677                assert allclose(l_v, 1), msg
678
679                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
680               
681                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
682                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
683                msg += ' Inner product is %e.' %x               
684                assert x < epsilon, msg
685
686
687
688
689        #Check neighbour structure
690        for i in range(N):
691            # For each triangle
692           
693            for k, neighbour_id in enumerate(self.neighbours[i,:]):
694
695                #Assert that my neighbour's neighbour is me
696                #Boundaries need not fulfill this
697                if neighbour_id >= 0:
698                    edge = self.neighbour_edges[i, k]
699                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
700                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
701                    assert self.neighbours[neighbour_id, edge] == i ,msg
702
703
704
705        #Check that all boundaries have
706        # unique, consecutive, negative indices
707
708        #L = len(self.boundary)
709        #for i in range(L):
710        #    id, edge = self.boundary_segments[i]
711        #    assert self.neighbours[id, edge] == -i-1
712
713
714        #NOTE: This assert doesn't hold true if there are internal boundaries
715        #FIXME: Look into this further.
716        #FIXME (Ole): In pyvolution mark 3 this is OK again
717        #NOTE: No longer works because neighbour structure is modified by
718        #      domain set_boundary.
719        #for id, edge in self.boundary:
720        #    assert self.neighbours[id,edge] < 0
721        #
722        #NOTE (Ole): I reckon this was resolved late 2004?
723        #
724        #See domain.set_boundary
725
726
727
728        #Check integrity of inverted triangle structure
729
730        V = self.vertex_value_indices[:] #Take a copy
731        V = sort(V)
732        assert allclose(V, range(3*N))
733
734        assert sum(self.number_of_triangles_per_node) ==\
735               len(self.vertex_value_indices)
736
737        # Check number of triangles per node
738        count = [0]*self.number_of_nodes
739        for triangle in self.triangles:
740            for i in triangle:
741                count[i] += 1
742
743        assert allclose(count, self.number_of_triangles_per_node)
744
745
746        # Check integrity of vertex_value_indices
747        current_node = 0
748        k = 0 # Track triangles touching on node
749        for index in self.vertex_value_indices:
750
751            if self.number_of_triangles_per_node[current_node] == 0:
752                # Node is lone - i.e. not part of the mesh
753                continue
754           
755            k += 1
756           
757            volume_id = index / 3
758            vertex_id = index % 3
759           
760            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
761                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
762            assert self.triangles[volume_id, vertex_id] == current_node, msg
763                       
764            if self.number_of_triangles_per_node[current_node] == k:
765                # Move on to next node
766                k = 0
767                current_node += 1
768
769
770    def get_lone_vertices(self):
771        """Return a list of vertices that are not connected to any triangles.
772
773        """
774        return self.lone_vertices
775
776    def get_centroid_coordinates(self, absolute=False):
777        """Return all centroid coordinates.
778        Return all centroid coordinates for all triangles as an Nx2 array
779        (ordered as x0, y0 for each triangle)
780
781        Boolean keyword argument absolute determines whether coordinates
782        are to be made absolute by taking georeference into account
783        Default is False as many parts of ANUGA expects relative coordinates.
784        """
785
786        V = self.centroid_coordinates
787        if absolute is True:
788            if not self.geo_reference.is_absolute():
789                V = self.geo_reference.get_absolute(V)
790           
791        return V
792
793       
794    def get_radii(self):
795        """Return all radii.
796        Return radius of inscribed cirle for all triangles
797        """
798        return self.radii       
799
800
801
802    def statistics(self):
803        """Output statistics about mesh
804        """
805
806        from Numeric import arange
807        from anuga.utilities.numerical_tools import histogram, create_bins
808
809        vertex_coordinates = self.vertex_coordinates # Relative coordinates
810        areas = self.areas
811        x = vertex_coordinates[:,0]
812        y = vertex_coordinates[:,1]
813
814
815        #Setup 10 bins for area histogram
816        bins = create_bins(areas, 10)
817        #m = max(areas)
818        #bins = arange(0., m, m/10)
819        hist = histogram(areas, bins)
820
821        str =  '------------------------------------------------\n'
822        str += 'Mesh statistics:\n'
823        str += '  Number of triangles = %d\n' %len(self)
824        str += '  Extent [m]:\n'
825        str += '    x in [%f, %f]\n' %(min(x), max(x))
826        str += '    y in [%f, %f]\n' %(min(y), max(y))
827        str += '  Areas [m^2]:\n'
828        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
829        str += '    number of distinct areas: %d\n' %(len(areas))       
830        str += '    Histogram:\n'
831
832        hi = bins[0]
833        for i, count in enumerate(hist):
834            lo = hi
835            if i+1 < len(bins):
836                #Open upper interval               
837                hi = bins[i+1]
838                str += '      [%f, %f[: %d\n' %(lo, hi, count)               
839            else:
840                #Closed upper interval
841                hi = max(areas)
842                str += '      [%f, %f]: %d\n' %(lo, hi, count)
843
844        N = len(areas)
845        if N > 10:
846            str += '    Percentiles (10%):\n'
847            areas = areas.tolist()
848            areas.sort()
849
850            k = 0
851            lower = min(areas)
852            for i, a in enumerate(areas):       
853                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
854                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
855                    lower = a
856                    k = i
857                   
858            str += '      %d triangles in [%f, %f]\n'\
859                   %(N-k, lower, max(areas))                   
860               
861                     
862        str += '  Boundary:\n'
863        str += '    Number of boundary segments == %d\n' %(len(self.boundary))
864        str += '    Boundary tags == %s\n' %self.get_boundary_tags() 
865        str += '------------------------------------------------\n'
866       
867
868        return str
869   
870
871    def get_triangle_containing_point(self, point):
872        """Return triangle id for triangle containing specifiend point (x,y)
873
874        If point isn't within mesh, raise exception
875
876        """
877       
878        # FIXME(Ole): This function is currently brute force
879        # because I needed it for diagnostics.
880        # We should make it fast - probably based on the
881        # quad tree structure.
882        from anuga.utilities.polygon import is_outside_polygon,\
883             is_inside_polygon
884
885        polygon = self.get_boundary_polygon()
886        #print polygon, point
887       
888        if is_outside_polygon(point, polygon):
889            msg = 'Point %s is outside mesh' %str(point)
890            raise Exception, msg
891
892
893        V = self.get_vertex_coordinates(absolute=True)
894
895        # FIXME: Horrible brute force
896        for i, triangle in enumerate(self.triangles):
897            poly = V[3*i:3*i+3]
898            #print i, poly
899
900            if is_inside_polygon(point, poly, closed=True):
901                return i
902               
903        return
904
905
906    def _get_intersecting_segments(self, line,
907                                   verbose=False):
908      """Find edges intersected by line
909
910      Input:
911          line - list of two points forming a segmented line
912          verbose
913      Output:
914          list of instances of class Triangle_intersection
915
916      This method is used by the public method
917      get_intersecting_segments(self, polyline) which also contains
918      more documentation.
919      """
920
921      from anuga.utilities.polygon import intersection
922      from anuga.utilities.polygon import is_inside_polygon
923     
924      msg = 'Line segment must contain exactly two points'
925      assert len(line) == 2, msg
926
927      # Origin of intersecting line to be used for
928      # establishing direction
929      xi0 = line[0][0]
930      eta0 = line[0][1]
931
932     
933      # Check intersection with edge segments for all triangles
934      # FIXME (Ole): This should be implemented in C
935      V = self.get_vertex_coordinates()
936      N = len(self)
937      triangle_intersections={} # Keep track of segments already done
938      for i in range(N):
939          # Get nodes and edge segments for each triangle
940          x0, y0 = V[3*i, :]
941          x1, y1 = V[3*i+1, :]
942          x2, y2 = V[3*i+2, :]
943           
944
945          edge_segments = [[[x0,y0], [x1, y1]],
946                            [[x1,y1], [x2, y2]],
947                            [[x2,y2], [x0, y0]]]
948
949          # Find segments that are intersected by line
950         
951          intersections = {} # Use dictionary to record points only once
952          for edge in edge_segments:
953
954              status, value = intersection(line, edge)
955              #if value is not None: print 'Triangle %d, Got' %i, status, value
956                 
957              if status == 1:
958                  # Normal intersection of one edge or vertex
959                  intersections[tuple(value)] = i                 
960
961                  # Exclude singular intersections with vertices
962                  #if not(allclose(value, edge[0]) or\
963                  #       allclose(value, edge[1])):
964                  #    intersections.append(value)
965
966              if status == 2:
967                  # Edge is sharing a segment with line
968
969                  # This is usually covered by the two
970                  # vertices that would have been picked up
971                  # under status == 1.
972                  # However, if coinciding line stops partway
973                  # along this edge, it will be recorded here.
974                  intersections[tuple(value[0,:])] = i
975                  intersections[tuple(value[1,:])] = i                                   
976
977                 
978          if len(intersections) == 1:
979              # Check if either line end point lies fully within this triangle
980              # If this is the case accept that as one end of the intersecting
981              # segment
982
983              poly = V[3*i:3*i+3]
984              if is_inside_polygon(line[1], poly, closed=False):
985                  intersections[tuple(line[1])] = i
986              elif is_inside_polygon(line[0], poly, closed=False):
987                  intersections[tuple(line[0])] = i         
988              else:
989                  # Ignore situations where one vertex is touch, for instance                   
990                  continue
991
992
993          msg = 'There can be only two or no intersections'
994          assert len(intersections) in [0,2], msg
995
996
997          if len(intersections) == 2:
998
999              # Calculate attributes for this segment
1000
1001
1002              # End points of intersecting segment
1003              points = intersections.keys()
1004              x0, y0 = points[0]
1005              x1, y1 = points[1]
1006
1007
1008              # Determine which end point is closer to the origin of the line
1009              # This is necessary for determining the direction of
1010              # the line and the normals
1011
1012              # Distances from line origin to the two intersections
1013              z0 = array([x0 - xi0, y0 - eta0])
1014              z1 = array([x1 - xi0, y1 - eta0])             
1015              d0 = sqrt(sum(z0**2)) 
1016              d1 = sqrt(sum(z1**2))
1017                 
1018              if d1 < d0:
1019                  # Swap
1020                  xi, eta = x0, y0
1021                  x0, y0 = x1, y1
1022                  x1, y1 = xi, eta
1023
1024              # (x0,y0) is now the origin of the intersecting segment
1025                 
1026
1027              # Normal direction:
1028              # Right hand side relative to line direction
1029              vector = array([x1 - x0, y1 - y0]) # Segment vector
1030              length = sqrt(sum(vector**2))      # Segment length
1031              normal = array([vector[1], -vector[0]])/length
1032
1033
1034              segment = ((x0,y0), (x1, y1))   
1035              T = Triangle_intersection(segment=segment,
1036                                        normal=normal,
1037                                        length=length,
1038                                        triangle_id=i)
1039
1040
1041              # Add segment unless it was done earlier
1042              if not triangle_intersections.has_key(segment):   
1043                  triangle_intersections[segment] = T
1044
1045
1046      # Return segments as a list           
1047      return triangle_intersections.values()
1048
1049
1050
1051    def get_intersecting_segments(self, polyline,
1052                                  verbose=False):
1053      """Find edges intersected by polyline
1054
1055      Input:
1056          polyline - list of points forming a segmented line
1057          verbose
1058
1059      Output:
1060          list of instances of class Triangle_intersection
1061
1062      The polyline may break inside any triangle causing multiple
1063      segments per triangle - consequently the same triangle may
1064      appear in several entries.
1065
1066      If a polyline segment coincides with a triangle edge,
1067      the the entire shared segment will be used.
1068      Onle one of the triangles thus intersected will be used and that
1069      is the first one encoutered.
1070
1071      Intersections with single vertices are ignored.
1072
1073      Resulting segments are unsorted
1074      """
1075
1076      msg = 'Polyline must contain at least two points'
1077      assert len(polyline) >= 2, msg
1078
1079      # For all segments in polyline
1080      triangle_intersections = []
1081      for i, point0 in enumerate(polyline[:-1]):
1082
1083          point1 = polyline[i+1]
1084          if verbose:
1085              print 'Extracting mesh intersections from line:',
1086              print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0[0], point0[1],
1087                                                    point1[0], point1[1])
1088             
1089         
1090          line = [point0, point1]
1091
1092          triangle_intersections += self._get_intersecting_segments(line)
1093
1094
1095      return triangle_intersections
1096 
1097
1098 
1099
1100    def get_triangle_neighbours(self, tri_id):
1101        """ Given a triangle id, Return an array of the
1102        3 neighbour triangle id's.
1103
1104        Negative returned triangle id's represent a boundary as a neighbour.
1105       
1106        If the given triangle id is bad, return an empty list.
1107        """
1108
1109        try:
1110            return self.neighbours[tri_id,:]
1111        except IndexError:
1112            return []
1113       
1114
1115
1116class Triangle_intersection:
1117    """Store information about line segments intersecting a triangle
1118   
1119    Attributes are
1120
1121        segment: Line segment intersecting triangle [[x0,y0], [x1, y1]]
1122        normal: [a,b] right hand normal to segment
1123        length: Length of intersecting segment
1124        triangle_id: id (in mesh) of triangle being intersected
1125
1126    """
1127
1128
1129    def __init__(self,
1130                 segment=None,
1131                 normal=None,
1132                 length=None,
1133                 triangle_id=None):
1134        self.segment = segment             
1135        self.normal = normal
1136        self.length = length
1137        self.triangle_id = triangle_id
1138       
1139
1140    def __repr__(self):
1141        s = 'Triangle_intersection('
1142        s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\
1143             %(self.segment,
1144               self.normal,
1145               self.length,
1146               self.triangle_id)
1147   
1148        return s
Note: See TracBrowser for help on using the repository browser.