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

Last change on this file since 7317 was 7317, checked in by rwilson, 15 years ago

Replaced 'print' statements with log.critical() calls.

File size: 41.7 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.utilities.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        return
887
888
889
890    def get_intersecting_segments(self, polyline,
891                                  use_cache=False,
892                                  verbose=False):
893        """Find edges intersected by polyline
894
895        Input:
896            polyline - list of points forming a segmented line
897            use_cache
898            verbose
899
900        Output:
901            list of instances of class Triangle_intersection
902
903        The polyline may break inside any triangle causing multiple
904        segments per triangle - consequently the same triangle may
905        appear in several entries.
906
907        If a polyline segment coincides with a triangle edge,
908        the the entire shared segment will be used.
909        Onle one of the triangles thus intersected will be used and that
910        is the first one encountered.
911
912        Intersections with single vertices are ignored.
913
914        Resulting segments are unsorted
915        """
916
917        V = self.get_vertex_coordinates()
918        N = len(self)
919
920        # Adjust polyline to mesh spatial origin
921        polyline = self.geo_reference.get_relative(polyline)
922
923        if use_cache is True:
924            segments = cache(get_intersecting_segments,
925                             (V, N, polyline),
926                             {'verbose': verbose},
927                             verbose=verbose)
928        else:
929            segments = get_intersecting_segments(V, N, polyline,
930                                                 verbose=verbose)
931
932
933        return segments
934
935
936
937    def get_triangle_neighbours(self, tri_id):
938        """ Given a triangle id, Return an array of the
939        3 neighbour triangle id's.
940
941        Negative returned triangle id's represent a boundary as a neighbour.
942
943        If the given triangle id is bad, return an empty list.
944        """
945
946        try:
947            return self.neighbours[tri_id,:]
948        except IndexError:
949            return []
950
951
952    def get_interpolation_object(self):
953        """Get object I that will allow linear interpolation using this mesh
954
955        This is a time consuming process but it needs only to be
956        once for the mesh.
957
958        Interpolation can then be done using
959
960        result = I.interpolate_block(vertex_values, interpolation_points)
961
962        where vertex values have been obtained from a quantity using
963        vertex_values, triangles = self.get_vertex_values()
964        """
965
966        if hasattr(self, 'interpolation_object'):
967            I = self.interpolation_object
968        else:
969            from anuga.fit_interpolate.interpolate import Interpolate
970
971            # Get discontinuous mesh - this will match internal
972            # representation of vertex values
973            triangles = self.get_disconnected_triangles()
974            vertex_coordinates = self.get_vertex_coordinates()
975
976            I = Interpolate(vertex_coordinates, triangles)
977            self.interpolation_object = I
978
979        return I
980
981
982class Triangle_intersection:
983    """Store information about line segments intersecting a triangle
984
985    Attributes are
986
987        segment: Line segment intersecting triangle [[x0,y0], [x1, y1]]
988        normal: [a,b] right hand normal to segment
989        length: Length of intersecting segment
990        triangle_id: id (in mesh) of triangle being intersected
991
992    """
993
994
995    def __init__(self,
996                 segment=None,
997                 normal=None,
998                 length=None,
999                 triangle_id=None):
1000        self.segment = segment
1001        self.normal = normal
1002        self.length = length
1003        self.triangle_id = triangle_id
1004
1005
1006    def __repr__(self):
1007        s = 'Triangle_intersection('
1008        s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\
1009             %(self.segment,
1010               self.normal,
1011               self.length,
1012               self.triangle_id)
1013
1014        return s
1015
1016
1017
1018def _get_intersecting_segments(V, N, line,
1019                               verbose=False):
1020    """Find edges intersected by line
1021
1022    Input:
1023        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1024        N: Number of triangles in mesh
1025        line - list of two points forming a segmented line
1026        verbose
1027    Output:
1028        list of instances of class Triangle_intersection
1029
1030    This method is used by the public method
1031    get_intersecting_segments(self, polyline) which also contains
1032    more documentation.
1033    """
1034
1035    from anuga.utilities.polygon import intersection
1036    from anuga.utilities.polygon import is_inside_polygon
1037
1038    msg = 'Line segment must contain exactly two points'
1039    assert len(line) == 2, msg
1040
1041    # Origin of intersecting line to be used for
1042    # establishing direction
1043    xi0 = line[0][0]
1044    eta0 = line[0][1]
1045
1046
1047    # Check intersection with edge segments for all triangles
1048    # FIXME (Ole): This should be implemented in C
1049    triangle_intersections={} # Keep track of segments already done
1050    for i in range(N):
1051        # Get nodes and edge segments for each triangle
1052        x0, y0 = V[3*i, :]
1053        x1, y1 = V[3*i+1, :]
1054        x2, y2 = V[3*i+2, :]
1055
1056
1057        edge_segments = [[[x0,y0], [x1, y1]],
1058                          [[x1,y1], [x2, y2]],
1059                          [[x2,y2], [x0, y0]]]
1060
1061        # Find segments that are intersected by line
1062
1063        intersections = {} # Use dictionary to record points only once
1064        for edge in edge_segments:
1065
1066            status, value = intersection(line, edge)
1067            #if value is not None: log.critical('Triangle %d, status=%s, '
1068            #                                   'value=%s'
1069            #                                   % (i, str(status), str(value)))
1070
1071            if status == 1:
1072                # Normal intersection of one edge or vertex
1073                intersections[tuple(value)] = i
1074
1075                # Exclude singular intersections with vertices
1076                #if not(allclose(value, edge[0]) or\
1077                #       allclose(value, edge[1])):
1078                #    intersections.append(value)
1079
1080            if status == 2:
1081                # Edge is sharing a segment with line
1082
1083                # This is usually covered by the two
1084                # vertices that would have been picked up
1085                # under status == 1.
1086                # However, if coinciding line stops partway
1087                # along this edge, it will be recorded here.
1088                intersections[tuple(value[0,:])] = i
1089                intersections[tuple(value[1,:])] = i
1090
1091
1092        if len(intersections) == 1:
1093            # Check if either line end point lies fully within this triangle
1094            # If this is the case accept that as one end of the intersecting
1095            # segment
1096
1097            poly = V[3*i:3*i+3]
1098            if is_inside_polygon(line[1], poly, closed=False):
1099                intersections[tuple(line[1])] = i
1100            elif is_inside_polygon(line[0], poly, closed=False):
1101                intersections[tuple(line[0])] = i
1102            else:
1103                # Ignore situations where one vertex is touch, for instance
1104                continue
1105
1106
1107        msg = 'There can be only two or no intersections'
1108        assert len(intersections) in [0,2], msg
1109
1110
1111        if len(intersections) == 2:
1112
1113            # Calculate attributes for this segment
1114
1115
1116            # End points of intersecting segment
1117            points = intersections.keys()
1118            x0, y0 = points[0]
1119            x1, y1 = points[1]
1120
1121
1122            # Determine which end point is closer to the origin of the line
1123            # This is necessary for determining the direction of
1124            # the line and the normals
1125
1126            # Distances from line origin to the two intersections
1127            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
1128            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
1129            d0 = num.sqrt(num.sum(z0**2))
1130            d1 = num.sqrt(num.sum(z1**2))
1131
1132            if d1 < d0:
1133                # Swap
1134                xi, eta = x0, y0
1135                x0, y0 = x1, y1
1136                x1, y1 = xi, eta
1137
1138            # (x0,y0) is now the origin of the intersecting segment
1139
1140
1141            # Normal direction:
1142            # Right hand side relative to line direction
1143            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
1144            length = num.sqrt(num.sum(vector**2))      # Segment length
1145            normal = num.array([vector[1], -vector[0]], num.float)/length
1146
1147
1148            segment = ((x0,y0), (x1, y1))
1149            T = Triangle_intersection(segment=segment,
1150                                      normal=normal,
1151                                      length=length,
1152                                      triangle_id=i)
1153
1154
1155            # Add segment unless it was done earlier
1156            if not triangle_intersections.has_key(segment):
1157                triangle_intersections[segment] = T
1158
1159
1160    # Return segments as a list
1161    return triangle_intersections.values()
1162
1163
1164def get_intersecting_segments(V, N, polyline,
1165                              verbose=False):
1166    """Internal function to find edges intersected by Polyline
1167
1168    Input:
1169        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1170        N: Number of triangles in mesh
1171        polyline - list of points forming a segmented line
1172        verbose
1173    Output:
1174        list of instances of class Triangle_intersection
1175
1176    This method is used by the public method
1177    get_intersecting_segments(self, polyline) which also contains
1178    more documentation.
1179    """
1180
1181    msg = 'Polyline must contain at least two points'
1182    assert len(polyline) >= 2, msg
1183
1184
1185    # For all segments in polyline
1186    triangle_intersections = []
1187    for i, point0 in enumerate(polyline[:-1]):
1188
1189        point1 = polyline[i+1]
1190        if verbose:
1191            log.critical('Extracting mesh intersections from line:')
1192            log.critical('(%.2f, %.2f) - (%.2f, %.2f)'
1193                         % (point0[0], point0[1], point1[0], point1[1]))
1194
1195        line = [point0, point1]
1196        triangle_intersections += _get_intersecting_segments(V, N, line,
1197                                                             verbose=verbose)
1198
1199
1200    msg = 'No segments found'
1201    assert len(triangle_intersections) > 0, msg
1202
1203
1204    return triangle_intersections
1205
1206
1207
1208
1209
1210def segment_midpoints(segments):
1211    """Calculate midpoints of all segments
1212
1213    Inputs:
1214       segments: List of instances of class Segment
1215
1216    Ouputs:
1217       midpoints: List of points
1218    """
1219
1220    midpoints = []
1221    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
1222    for segment in segments:
1223        assert isinstance(segment, Triangle_intersection), msg
1224
1225        midpoint = num.sum(num.array(segment.segment, num.float), axis=0)/2
1226        midpoints.append(midpoint)
1227
1228    return midpoints
1229
1230
1231
Note: See TracBrowser for help on using the repository browser.