source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py @ 7968

Last change on this file since 7968 was 7968, checked in by steve, 14 years ago

Commiting changes of refactor of culvert operator

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