# source:branches/numpy/anuga/abstract_2d_finite_volumes/neighbour_mesh.py@6428

Last change on this file since 6428 was 6428, checked in by rwilson, 14 years ago

numpy changes

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