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

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

Initial commit of numpy changes. Still a long way to go.

File size: 42.4 KB
Line
1"""Classes implementing general 2D triangular mesh with neighbour structure.
2
3   This structure is purely geometrical. Anything relating to quantities
4   or timestepping is implemented in subclass domain.py.
5
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia, 2004
8"""
9
10from general_mesh import General_mesh
11from 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
470        # Get mesh extent
471        xmin, xmax, ymin, ymax = self.get_extent(absolute=True)
472        pmin = ensure_numeric([xmin, ymin])
473        pmax = ensure_numeric([xmax, ymax])
474
475
476        # Assemble dictionary of boundary segments and choose starting point
477        segments = {}
478        inverse_segments = {}
479        p0 = None
480        mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh
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
491            # Take the point closest to pmin as starting point
492            # Note: Could be arbitrary, but nice to have
493            # a unique way of selecting
494            dist_A = num.sqrt(num.sum((A-pmin)**2))
495            dist_B = num.sqrt(num.sum((B-pmin)**2))
496
497            # Find lower leftmost point
498            if dist_A < mindist:
499                mindist = dist_A
500                p0 = A
501            if dist_B < mindist:
502                mindist = dist_B
503                p0 = B
504
505
506            # Sanity check
507            if p0 is None:
508                msg = 'Impossible: p0 is None!?'
509                raise Exception, msg
510
511            # Register potential paths from A to B
512            if not segments.has_key(tuple(A)):
513                segments[tuple(A)] = [] # Empty list for candidate points
514
515            segments[tuple(A)].append(B)
516
517
519        polygon = [list(p0)]# Storage for final boundary polygon
520        point_registry = {} # Keep track of storage to avoid multiple runs
521                            # around boundary. This will only be the case if
522                            # there are more than one candidate.
523                            # FIXME (Ole): Perhaps we can do away with polygon
524                            # and use only point_registry to save space.
525
526        point_registry[tuple(p0)] = 0
527
528        while len(point_registry) < len(self.boundary):
529
530            candidate_list = segments[tuple(p0)]
531            if len(candidate_list) > 1:
532                # Multiple points detected (this will be the case for meshes
533                # with duplicate points as those used for discontinuous
534                # triangles with vertices stored uniquely).
535                # Take the candidate that is furthest to the clockwise
536                # direction, as that will follow the boundary.
537                #
538                # This will also be the case for pathological triangles
539                # that have no neighbours.
540
541                if verbose:
542                    print 'Point %s has multiple candidates: %s'\
543                          %(str(p0), candidate_list)
544
545                # Check that previous are not in candidate list
546                #for p in candidate_list:
547                #    assert not allclose(p0, p)
548
549                # Choose vector against which all angles will be measured
550                if len(polygon) > 1:
551                    v_prev = p0 - polygon[-2] # Vector that leads to p0
552                                              # from previous point
553                else:
554                    # FIXME (Ole): What do we do if the first point has
555                    # multiple candidates?
556                    # Being the lower left corner, perhaps we can use the
557                    # vector [1, 0], but I really don't know if this is
558                    # completely watertight.
559                    v_prev = [1.0, 0.0]
560
561
562                # Choose candidate with minimum angle
563                minimum_angle = 2*pi
564                for pc in candidate_list:
565
566                    vc = pc-p0  # Candidate vector (from p0 to candidate pt)
567
568                    # Angle between each candidate and the previous vector
569                    # in [-pi, pi]
570                    ac = angle(vc, v_prev)
571                    if ac > pi:
572                        # Give preference to angles on the right hand side
573                        # of v_prev
574                        # print 'pc = %s, changing angle from %f to %f'\
575                        # %(pc, ac*180/pi, (ac-2*pi)*180/pi)
576                        ac = ac-2*pi
577
578                    # Take the minimal angle corresponding to the
579                    # rightmost vector
580                    if ac < minimum_angle:
581                        minimum_angle = ac
582                        p1 = pc             # Best candidate
583
584
585                if verbose is True:
586                    print '  Best candidate %s, angle %f'\
587                          %(p1, minimum_angle*180/pi)
588
589            else:
590                p1 = candidate_list
591
592
593            if point_registry.has_key(tuple(p1)):
594                # We have reached a point already visited.
595
596                if num.allclose(p1, polygon):
597                    # If it is the initial point, the polygon is complete.
598
599                    if verbose is True:
600                        print '  Stop criterion fulfilled at point %s' %p1
601                        print polygon
602
603                    # We have completed the boundary polygon - yeehaa
604                    break
605                else:
606                    # The point already visited is not the initial point
607                    # This would be a pathological triangle, but the
608                    # algorithm must be able to deal with this
609                    pass
610
611            else:
612                # We are still finding new points on the boundary
613                point_registry[tuple(p1)] = len(point_registry)
614
615            polygon.append(list(p1)) # De-numeric each point :-)
616            p0 = p1
617
618
619        return polygon
620
621
622    def check_integrity(self):
623        """Check that triangles are internally consistent e.g.
624        that area corresponds to edgelengths, that vertices
625        are arranged in a counter-clockwise order, etc etc
626        Neighbour structure will be checked by class Mesh
627        """
628
629        from anuga.config import epsilon
630        from anuga.utilities.numerical_tools import anglediff
631
632        N = len(self)
633
634        # Get x,y coordinates for all vertices for all triangles
635        V = self.get_vertex_coordinates()
636
637        # Check each triangle
638        for i in range(N):
639
640            x0, y0 = V[3*i, :]
641            x1, y1 = V[3*i+1, :]
642            x2, y2 = V[3*i+2, :]
643
644            # Check that area hasn't been compromised
645            area = self.areas[i]
646            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
647            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
648                  %(x0,y0,x1,y1,x2,y2)
649            assert abs((area - ref)/area) < epsilon, msg
650
651            # Check that points are arranged in counter clock-wise order
652            v0 = [x1-x0, y1-y0]
653            v1 = [x2-x1, y2-y1]
654            v2 = [x0-x2, y0-y2]
655            a0 = anglediff(v1, v0)
656            a1 = anglediff(v2, v1)
657            a2 = anglediff(v0, v2)
658
659            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
660            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
661            assert a0 < pi and a1 < pi and a2 < pi, msg
662
663            # Check that normals are orthogonal to edge vectors
664            # Note that normal[k] lies opposite vertex k
665
666            normal0 = self.normals[i, 0:2]
667            normal1 = self.normals[i, 2:4]
668            normal2 = self.normals[i, 4:6]
669
670            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
671
672                # Normalise
673                l_u = num.sqrt(u*u + u*u)
674                l_v = num.sqrt(v*v + v*v)
675
676                msg = 'Normal vector in triangle %d does not have unit length' %i
677                assert num.allclose(l_v, 1), msg
678
679                x = (u*v + u*v)/l_u # Inner product
680
681                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
682                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
683                msg += ' Inner product is %e.' %x
684                assert x < epsilon, msg
685
686
687
688
689        #Check neighbour structure
690        for i in range(N):
691            # For each triangle
692
693            for k, neighbour_id in enumerate(self.neighbours[i,:]):
694
695                #Assert that my neighbour's neighbour is me
696                #Boundaries need not fulfill this
697                if neighbour_id >= 0:
698                    edge = self.neighbour_edges[i, k]
699                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
700                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
701                    assert self.neighbours[neighbour_id, edge] == i ,msg
702
703
704
705        #Check that all boundaries have
706        # unique, consecutive, negative indices
707
708        #L = len(self.boundary)
709        #for i in range(L):
710        #    id, edge = self.boundary_segments[i]
711        #    assert self.neighbours[id, edge] == -i-1
712
713
714        #NOTE: This assert doesn't hold true if there are internal boundaries
715        #FIXME: Look into this further.
716        #FIXME (Ole): In pyvolution mark 3 this is OK again
717        #NOTE: No longer works because neighbour structure is modified by
718        #      domain set_boundary.
719        #for id, edge in self.boundary:
720        #    assert self.neighbours[id,edge] < 0
721        #
722        #NOTE (Ole): I reckon this was resolved late 2004?
723        #
724        #See domain.set_boundary
725
726
727
728        #Check integrity of inverted triangle structure
729
730        V = self.vertex_value_indices[:] #Take a copy
731        V = num.sort(V)
732        assert num.allclose(V, range(3*N))
733
734        assert num.sum(self.number_of_triangles_per_node) ==\
735                       len(self.vertex_value_indices)
736
737        # Check number of triangles per node
738        count = *self.number_of_nodes
739        for triangle in self.triangles:
740            for i in triangle:
741                count[i] += 1
742
743        assert num.allclose(count, self.number_of_triangles_per_node)
744
745
746        # Check integrity of vertex_value_indices
747        current_node = 0
748        k = 0 # Track triangles touching on node
749        for index in self.vertex_value_indices:
750
751            if self.number_of_triangles_per_node[current_node] == 0:
752                # Node is lone - i.e. not part of the mesh
753                continue
754
755            k += 1
756
757            volume_id = index / 3
758            vertex_id = index % 3
759
760            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
761                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
762            assert self.triangles[volume_id, vertex_id] == current_node, msg
763
764            if self.number_of_triangles_per_node[current_node] == k:
765                # Move on to next node
766                k = 0
767                current_node += 1
768
769
770    def get_lone_vertices(self):
771        """Return a list of vertices that are not connected to any triangles.
772
773        """
774        return self.lone_vertices
775
776    def get_centroid_coordinates(self, absolute=False):
777        """Return all centroid coordinates.
778        Return all centroid coordinates for all triangles as an Nx2 array
779        (ordered as x0, y0 for each triangle)
780
781        Boolean keyword argument absolute determines whether coordinates
782        are to be made absolute by taking georeference into account
783        Default is False as many parts of ANUGA expects relative coordinates.
784        """
785
786        V = self.centroid_coordinates
787        if absolute is True:
788            if not self.geo_reference.is_absolute():
789                V = self.geo_reference.get_absolute(V)
790
791        return V
792
793
796        Return radius of inscribed cirle for all triangles
797        """
799
800
801
802    def statistics(self):
804        """
805
806        from anuga.utilities.numerical_tools import histogram, create_bins
807
808        vertex_coordinates = self.vertex_coordinates # Relative coordinates
809        areas = self.areas
810        x = vertex_coordinates[:,0]
811        y = vertex_coordinates[:,1]
812
813
814        #Setup 10 bins for area histogram
815        bins = create_bins(areas, 10)
816        #m = max(areas)
817        #bins = arange(0., m, m/10)
818        hist = histogram(areas, bins)
819
820        str =  '------------------------------------------------\n'
821        str += 'Mesh statistics:\n'
822        str += '  Number of triangles = %d\n' %len(self)
823        str += '  Extent [m]:\n'
824        str += '    x in [%f, %f]\n' %(min(x), max(x))
825        str += '    y in [%f, %f]\n' %(min(y), max(y))
826        str += '  Areas [m^2]:\n'
827        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
828        str += '    number of distinct areas: %d\n' %(len(areas))
829        str += '    Histogram:\n'
830
831        hi = bins
832        for i, count in enumerate(hist):
833            lo = hi
834            if i+1 < len(bins):
835                #Open upper interval
836                hi = bins[i+1]
837                str += '      [%f, %f[: %d\n' %(lo, hi, count)
838            else:
839                #Closed upper interval
840                hi = max(areas)
841                str += '      [%f, %f]: %d\n' %(lo, hi, count)
842
843        N = len(areas)
844        if N > 10:
845            str += '    Percentiles (10%):\n'
846            areas = areas.tolist()
847            areas.sort()
848
849            k = 0
850            lower = min(areas)
851            for i, a in enumerate(areas):
852                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas
853                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
854                    lower = a
855                    k = i
856
857            str += '      %d triangles in [%f, %f]\n'\
858                   %(N-k, lower, max(areas))
859
860
861        str += '  Boundary:\n'
862        str += '    Number of boundary segments == %d\n' %(len(self.boundary))
863        str += '    Boundary tags == %s\n' %self.get_boundary_tags()
864        str += '------------------------------------------------\n'
865
866
867        return str
868
869
870    def get_triangle_containing_point(self, point):
871        """Return triangle id for triangle containing specifiend point (x,y)
872
873        If point isn't within mesh, raise exception
874
875        """
876
877        # FIXME(Ole): This function is currently brute force
878        # because I needed it for diagnostics.
879        # We should make it fast - probably based on the
881        from anuga.utilities.polygon import is_outside_polygon,\
882             is_inside_polygon
883
884        polygon = self.get_boundary_polygon()
885        #print polygon, point
886
887        if is_outside_polygon(point, polygon):
888            msg = 'Point %s is outside mesh' %str(point)
889            raise Exception, msg
890
891
892        V = self.get_vertex_coordinates(absolute=True)
893
894        # FIXME: Horrible brute force
895        for i, triangle in enumerate(self.triangles):
896            poly = V[3*i:3*i+3]
897            #print i, poly
898
899            if is_inside_polygon(point, poly, closed=True):
900                return i
901
902        return
903
904
905
906    def get_intersecting_segments(self, polyline,
907                                  use_cache=False,
908                                  verbose=False):
909        """Find edges intersected by polyline
910
911        Input:
912            polyline - list of points forming a segmented line
913            use_cache
914            verbose
915
916        Output:
917            list of instances of class Triangle_intersection
918
919        The polyline may break inside any triangle causing multiple
920        segments per triangle - consequently the same triangle may
921        appear in several entries.
922
923        If a polyline segment coincides with a triangle edge,
924        the the entire shared segment will be used.
925        Onle one of the triangles thus intersected will be used and that
926        is the first one encountered.
927
928        Intersections with single vertices are ignored.
929
930        Resulting segments are unsorted
931        """
932
933        V = self.get_vertex_coordinates()
934        N = len(self)
935
936        # Adjust polyline to mesh spatial origin
937        polyline = self.geo_reference.get_relative(polyline)
938
939        if use_cache is True:
940            segments = cache(get_intersecting_segments,
941                             (V, N, polyline),
942                             {'verbose': verbose},
943                             verbose=verbose)
944        else:
945            segments = get_intersecting_segments(V, N, polyline,
946                                                 verbose=verbose)
947
948
949        return segments
950
951
952
953    def get_triangle_neighbours(self, tri_id):
954        """ Given a triangle id, Return an array of the
955        3 neighbour triangle id's.
956
957        Negative returned triangle id's represent a boundary as a neighbour.
958
959        If the given triangle id is bad, return an empty list.
960        """
961
962        try:
963            return self.neighbours[tri_id,:]
964        except IndexError:
965            return []
966
967
968    def get_interpolation_object(self):
969        """Get object I that will allow linear interpolation using this mesh
970
971        This is a time consuming process but it needs only to be
972        once for the mesh.
973
974        Interpolation can then be done using
975
976        result = I.interpolate_block(vertex_values, interpolation_points)
977
978        where vertex values have been obtained from a quantity using
979        vertex_values, triangles = self.get_vertex_values()
980        """
981
982        if hasattr(self, 'interpolation_object'):
983            I = self.interpolation_object
984        else:
985            from anuga.fit_interpolate.interpolate import Interpolate
986
987            # Get discontinuous mesh - this will match internal
988            # representation of vertex values
989            triangles = self.get_disconnected_triangles()
990            vertex_coordinates = self.get_vertex_coordinates()
991
992            I = Interpolate(vertex_coordinates, triangles)
993            self.interpolation_object = I
994
995        return I
996
997
998class Triangle_intersection:
999    """Store information about line segments intersecting a triangle
1000
1001    Attributes are
1002
1003        segment: Line segment intersecting triangle [[x0,y0], [x1, y1]]
1004        normal: [a,b] right hand normal to segment
1005        length: Length of intersecting segment
1006        triangle_id: id (in mesh) of triangle being intersected
1007
1008    """
1009
1010
1011    def __init__(self,
1012                 segment=None,
1013                 normal=None,
1014                 length=None,
1015                 triangle_id=None):
1016        self.segment = segment
1017        self.normal = normal
1018        self.length = length
1019        self.triangle_id = triangle_id
1020
1021
1022    def __repr__(self):
1023        s = 'Triangle_intersection('
1024        s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\
1025             %(self.segment,
1026               self.normal,
1027               self.length,
1028               self.triangle_id)
1029
1030        return s
1031
1032
1033
1034def _get_intersecting_segments(V, N, line,
1035                               verbose=False):
1036    """Find edges intersected by line
1037
1038    Input:
1039        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1040        N: Number of triangles in mesh
1041        line - list of two points forming a segmented line
1042        verbose
1043    Output:
1044        list of instances of class Triangle_intersection
1045
1046    This method is used by the public method
1047    get_intersecting_segments(self, polyline) which also contains
1048    more documentation.
1049    """
1050
1051    from anuga.utilities.polygon import intersection
1052    from anuga.utilities.polygon import is_inside_polygon
1053
1054    msg = 'Line segment must contain exactly two points'
1055    assert len(line) == 2, msg
1056
1057    # Origin of intersecting line to be used for
1058    # establishing direction
1059    xi0 = line
1060    eta0 = line
1061
1062
1063    # Check intersection with edge segments for all triangles
1064    # FIXME (Ole): This should be implemented in C
1065    triangle_intersections={} # Keep track of segments already done
1066    for i in range(N):
1067        # Get nodes and edge segments for each triangle
1068        x0, y0 = V[3*i, :]
1069        x1, y1 = V[3*i+1, :]
1070        x2, y2 = V[3*i+2, :]
1071
1072
1073        edge_segments = [[[x0,y0], [x1, y1]],
1074                          [[x1,y1], [x2, y2]],
1075                          [[x2,y2], [x0, y0]]]
1076
1077        # Find segments that are intersected by line
1078
1079        intersections = {} # Use dictionary to record points only once
1080        for edge in edge_segments:
1081
1082            status, value = intersection(line, edge)
1083            #if value is not None: print 'Triangle %d, Got' %i, status, value
1084
1085            if status == 1:
1086                # Normal intersection of one edge or vertex
1087                intersections[tuple(value)] = i
1088
1089                # Exclude singular intersections with vertices
1090                #if not(allclose(value, edge) or\
1091                #       allclose(value, edge)):
1092                #    intersections.append(value)
1093
1094            if status == 2:
1095                # Edge is sharing a segment with line
1096
1097                # This is usually covered by the two
1098                # vertices that would have been picked up
1099                # under status == 1.
1100                # However, if coinciding line stops partway
1101                # along this edge, it will be recorded here.
1102                intersections[tuple(value[0,:])] = i
1103                intersections[tuple(value[1,:])] = i
1104
1105
1106        if len(intersections) == 1:
1107            # Check if either line end point lies fully within this triangle
1108            # If this is the case accept that as one end of the intersecting
1109            # segment
1110
1111            poly = V[3*i:3*i+3]
1112            if is_inside_polygon(line, poly, closed=False):
1113                intersections[tuple(line)] = i
1114            elif is_inside_polygon(line, poly, closed=False):
1115                intersections[tuple(line)] = i
1116            else:
1117                # Ignore situations where one vertex is touch, for instance
1118                continue
1119
1120
1121        msg = 'There can be only two or no intersections'
1122        assert len(intersections) in [0,2], msg
1123
1124
1125        if len(intersections) == 2:
1126
1127            # Calculate attributes for this segment
1128
1129
1130            # End points of intersecting segment
1131            points = intersections.keys()
1132            x0, y0 = points
1133            x1, y1 = points
1134
1135
1136            # Determine which end point is closer to the origin of the line
1137            # This is necessary for determining the direction of
1138            # the line and the normals
1139
1140            # Distances from line origin to the two intersections
1141            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
1142            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
1143            d0 = num.sqrt(num.sum(z0**2))
1144            d1 = num.sqrt(num.sum(z1**2))
1145
1146            if d1 < d0:
1147                # Swap
1148                xi, eta = x0, y0
1149                x0, y0 = x1, y1
1150                x1, y1 = xi, eta
1151
1152            # (x0,y0) is now the origin of the intersecting segment
1153
1154
1155            # Normal direction:
1156            # Right hand side relative to line direction
1157            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
1158            length = num.sqrt(num.sum(vector**2))      # Segment length
1159            normal = num.array([vector, -vector], num.float)/length
1160
1161
1162            segment = ((x0,y0), (x1, y1))
1163            T = Triangle_intersection(segment=segment,
1164                                      normal=normal,
1165                                      length=length,
1166                                      triangle_id=i)
1167
1168
1169            # Add segment unless it was done earlier
1170            if not triangle_intersections.has_key(segment):
1171                triangle_intersections[segment] = T
1172
1173
1174    # Return segments as a list
1175    return triangle_intersections.values()
1176
1177
1178def get_intersecting_segments(V, N, polyline,
1179                              verbose=False):
1180    """Internal function to find edges intersected by Polyline
1181
1182    Input:
1183        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1184        N: Number of triangles in mesh
1185        polyline - list of points forming a segmented line
1186        verbose
1187    Output:
1188        list of instances of class Triangle_intersection
1189
1190    This method is used by the public method
1191    get_intersecting_segments(self, polyline) which also contains
1192    more documentation.
1193    """
1194
1195    msg = 'Polyline must contain at least two points'
1196    assert len(polyline) >= 2, msg
1197
1198
1199    # For all segments in polyline
1200    triangle_intersections = []
1201    for i, point0 in enumerate(polyline[:-1]):
1202
1203        point1 = polyline[i+1]
1204        if verbose:
1205            print 'Extracting mesh intersections from line:',
1206            print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0, point0,
1207                                                  point1, point1)
1208
1209        line = [point0, point1]
1210        triangle_intersections += _get_intersecting_segments(V, N, line,
1211                                                             verbose=verbose)
1212
1213
1214    msg = 'No segments found'
1215    assert len(triangle_intersections) > 0, msg
1216
1217
1218    return triangle_intersections
1219
1220
1221
1222
1223
1224def segment_midpoints(segments):
1225    """Calculate midpoints of all segments
1226
1227    Inputs:
1228       segments: List of instances of class Segment
1229
1230    Ouputs:
1231       midpoints: List of points
1232    """
1233
1234    midpoints = []
1235    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
1236    for segment in segments:
1237        assert isinstance(segment, Triangle_intersection), msg
1238
1239        midpoint = num.sum(num.array(segment.segment, num.float))/2
1240        midpoints.append(midpoint)
1241
1242    return midpoints
1243
1244
1245
Note: See TracBrowser for help on using the repository browser.