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

Last change on this file since 8200 was 8200, checked in by steve, 13 years ago

Pulled the number_of_full_nodes and triangle out of generic_mesh up into generic_domain. These "numbers" are now calculated from ghost_recv and full_send data

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