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

Last change on this file since 8219 was 8219, checked in by steve, 11 years ago

Rudy created a mesh which produced a degenerate surrogate triangle. Changed throwing an error to defaulting to a first order reconstruction

File size: 43.1 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 = -((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
678            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
679            msg += 'Wrong area: %f  %f'\
680                  %(area, ref)
681            assert abs((area - ref)/area) < epsilon, msg
682
683            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
684            msg += ' is degenerate:  area == %f' % self.areas[i]
685            assert area > 0.0, msg
686
687            # Check that points are arranged in counter clock-wise order
688            v0 = [x1-x0, y1-y0]
689            v1 = [x2-x1, y2-y1]
690            v2 = [x0-x2, y0-y2]
691            a0 = anglediff(v1, v0)
692            a1 = anglediff(v2, v1)
693            a2 = anglediff(v0, v2)
694
695            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
696            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
697            assert a0 < pi and a1 < pi and a2 < pi, msg
698
699            # Check that normals are orthogonal to edge vectors
700            # Note that normal[k] lies opposite vertex k
701
702            normal0 = self.normals[i, 0:2]
703            normal1 = self.normals[i, 2:4]
704            normal2 = self.normals[i, 4:6]
705
706            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
707
708                # Normalise
709                l_u = num.sqrt(u[0]*u[0] + u[1]*u[1])
710                l_v = num.sqrt(v[0]*v[0] + v[1]*v[1])
711
712                msg = 'Normal vector in triangle %d does not have unit length' %i
713                assert num.allclose(l_v, 1), msg
714
715                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
716
717                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
718                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
719                msg += ' Inner product is %e.' %x
720                assert x < epsilon, msg
721
722
723
724
725        # Check neighbour structure
726        for i in range(N):
727            # For each triangle
728
729            for k, neighbour_id in enumerate(self.neighbours[i,:]):
730
731                #Assert that my neighbour's neighbour is me
732                #Boundaries need not fulfill this
733                if neighbour_id >= 0:
734                    edge = self.neighbour_edges[i, k]
735                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
736                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
737                    assert self.neighbours[neighbour_id, edge] == i ,msg
738
739
740
741        #Check that all boundaries have
742        # unique, consecutive, negative indices
743
744        #L = len(self.boundary)
745        #for i in range(L):
746        #    id, edge = self.boundary_segments[i]
747        #    assert self.neighbours[id, edge] == -i-1
748
749
750        #NOTE: This assert doesn't hold true if there are internal boundaries
751        #FIXME: Look into this further.
752        #FIXME (Ole): In pyvolution mark 3 this is OK again
753        #NOTE: No longer works because neighbour structure is modified by
754        #      domain set_boundary.
755        #for id, edge in self.boundary:
756        #    assert self.neighbours[id,edge] < 0
757        #
758        #NOTE (Ole): I reckon this was resolved late 2004?
759        #
760        #See domain.set_boundary
761
762
763
764        #Check integrity of inverted triangle structure
765
766        V = self.vertex_value_indices[:] #Take a copy
767        V = num.sort(V)
768        assert num.allclose(V, range(3*N))
769
770        assert num.sum(self.number_of_triangles_per_node) ==\
771                       len(self.vertex_value_indices)
772
773        # Check number of triangles per node
774        count = [0]*self.number_of_nodes
775        for triangle in self.triangles:
776            for i in triangle:
777                count[i] += 1
778
779        assert num.allclose(count, self.number_of_triangles_per_node)
780
781
782        # Check integrity of vertex_value_indices
783        current_node = 0
784        k = 0 # Track triangles touching on node
785        for index in self.vertex_value_indices:
786
787            if self.number_of_triangles_per_node[current_node] == 0:
788                # Node is lone - i.e. not part of the mesh
789                continue
790
791            k += 1
792
793            volume_id = index / 3
794            vertex_id = index % 3
795
796            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
797                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
798            assert self.triangles[volume_id, vertex_id] == current_node, msg
799
800            if self.number_of_triangles_per_node[current_node] == k:
801                # Move on to next node
802                k = 0
803                current_node += 1
804
805
806    def get_lone_vertices(self):
807        """Return a list of vertices that are not connected to any triangles.
808
809        """
810        return self.lone_vertices
811
812    def get_centroid_coordinates(self, absolute=False):
813        """Return all centroid coordinates.
814        Return all centroid coordinates for all triangles as an Nx2 array
815        (ordered as x0, y0 for each triangle)
816
817        Boolean keyword argument absolute determines whether coordinates
818        are to be made absolute by taking georeference into account
819        Default is False as many parts of ANUGA expects relative coordinates.
820        """
821
822        V = self.centroid_coordinates
823        if absolute is True:
824            if not self.geo_reference.is_absolute():
825                V = self.geo_reference.get_absolute(V)
826
827        return V
828
829
830    def get_radii(self):
831        """Return all radii.
832        Return radius of inscribed cirle for all triangles
833        """
834        return self.radii
835
836
837
838    def statistics(self):
839        """Output statistics about mesh
840        """
841
842        from anuga.utilities.numerical_tools import histogram, create_bins
843
844        vertex_coordinates = self.vertex_coordinates # Relative coordinates
845        areas = self.areas
846        x = vertex_coordinates[:,0]
847        y = vertex_coordinates[:,1]
848
849
850        #Setup 10 bins for area histogram
851        bins = create_bins(areas, 10)
852        #m = max(areas)
853        #bins = arange(0., m, m/10)
854        hist = histogram(areas, bins)
855
856        str =  '------------------------------------------------\n'
857        str += 'Mesh statistics:\n'
858        str += '  Number of triangles = %d\n' %len(self)
859        str += '  Extent [m]:\n'
860        str += '    x in [%f, %f]\n' %(min(x), max(x))
861        str += '    y in [%f, %f]\n' %(min(y), max(y))
862        str += '  Areas [m^2]:\n'
863        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
864        str += '    number of distinct areas: %d\n' %(len(areas))
865        str += '    Histogram:\n'
866
867        hi = bins[0]
868        for i, count in enumerate(hist):
869            lo = hi
870            if i+1 < len(bins):
871                #Open upper interval
872                hi = bins[i+1]
873                str += '      [%f, %f[: %d\n' %(lo, hi, count)
874            else:
875                #Closed upper interval
876                hi = max(areas)
877                str += '      [%f, %f]: %d\n' %(lo, hi, count)
878
879        N = len(areas)
880        if N > 10:
881            str += '    Percentiles (10%):\n'
882            areas = areas.tolist()
883            areas.sort()
884
885            k = 0
886            lower = min(areas)
887            for i, a in enumerate(areas):
888                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas
889                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
890                    lower = a
891                    k = i
892
893            str += '      %d triangles in [%f, %f]\n'\
894                   %(N-k, lower, max(areas))
895
896
897        str += '  Boundary:\n'
898        str += '    Number of boundary segments == %d\n' %(len(self.boundary))
899        str += '    Boundary tags == %s\n' %self.get_boundary_tags()
900        str += '------------------------------------------------\n'
901
902
903        return str
904
905
906    def get_triangle_containing_point(self, point):
907        """Return triangle id for triangle containing specified point (x,y)
908
909        If point isn't within mesh, raise exception
910
911        """
912
913        # FIXME(Ole): This function is currently brute force
914        # because I needed it for diagnostics.
915        # We should make it fast - probably based on the
916        # quad tree structure.
917        from anuga.geometry.polygon import is_outside_polygon,\
918             is_inside_polygon
919
920        polygon = self.get_boundary_polygon()
921
922        if is_outside_polygon(point, polygon):
923            msg = 'Point %s is outside mesh' %str(point)
924            raise Exception(msg)
925
926
927        V = self.get_vertex_coordinates(absolute=True)
928
929        # FIXME: Horrible brute force
930        for i, triangle in enumerate(self.triangles):
931            poly = V[3*i:3*i+3]
932
933            if is_inside_polygon(point, poly, closed=True):
934                return i
935
936        msg = 'Point %s not found within a triangle' %str(point)
937        raise Exception(msg)
938
939
940
941
942    def get_intersecting_segments(self, polyline,
943                                  use_cache=False,
944                                  verbose=False):
945        """Find edges intersected by polyline
946
947        Input:
948            polyline - list of points forming a segmented line
949            use_cache
950            verbose
951
952        Output:
953            list of instances of class Triangle_intersection
954
955        The polyline may break inside any triangle causing multiple
956        segments per triangle - consequently the same triangle may
957        appear in several entries.
958
959        If a polyline segment coincides with a triangle edge,
960        the the entire shared segment will be used.
961        Onle one of the triangles thus intersected will be used and that
962        is the first one encountered.
963
964        Intersections with single vertices are ignored.
965
966        Resulting segments are unsorted
967        """
968
969        V = self.get_vertex_coordinates()
970        N = len(self)
971
972        # Adjust polyline to mesh spatial origin
973        polyline = self.geo_reference.get_relative(polyline)
974
975        if use_cache is True:
976            segments = cache(get_intersecting_segments,
977                             (V, N, polyline),
978                             {'verbose': verbose},
979                             verbose=verbose)
980        else:
981            segments = get_intersecting_segments(V, N, polyline,
982                                                 verbose=verbose)
983
984
985        return segments
986
987
988
989    def get_triangle_neighbours(self, tri_id):
990        """ Given a triangle id, Return an array of the
991        3 neighbour triangle id's.
992
993        Negative returned triangle id's represent a boundary as a neighbour.
994
995        If the given triangle id is bad, return an empty list.
996        """
997
998        try:
999            return self.neighbours[tri_id,:]
1000        except IndexError:
1001            return []
1002
1003
1004    def get_interpolation_object(self):
1005        """Get object I that will allow linear interpolation using this mesh
1006
1007        This is a time consuming process but it needs only to be
1008        once for the mesh.
1009
1010        Interpolation can then be done using
1011
1012        result = I.interpolate_block(vertex_values, interpolation_points)
1013
1014        where vertex values have been obtained from a quantity using
1015        vertex_values, triangles = self.get_vertex_values()
1016        """
1017
1018        if hasattr(self, 'interpolation_object'):
1019            I = self.interpolation_object
1020        else:
1021            from anuga.fit_interpolate.interpolate import Interpolate
1022
1023            # Get discontinuous mesh - this will match internal
1024            # representation of vertex values
1025            triangles = self.get_disconnected_triangles()
1026            vertex_coordinates = self.get_vertex_coordinates()
1027
1028            I = Interpolate(vertex_coordinates, triangles)
1029            self.interpolation_object = I
1030
1031        return I
1032
1033
1034class Triangle_intersection:
1035    """Store information about line segments intersecting a triangle
1036
1037    Attributes are
1038
1039        segment: Line segment intersecting triangle [[x0,y0], [x1, y1]]
1040        normal: [a,b] right hand normal to segment
1041        length: Length of intersecting segment
1042        triangle_id: id (in mesh) of triangle being intersected
1043
1044    """
1045
1046
1047    def __init__(self,
1048                 segment=None,
1049                 normal=None,
1050                 length=None,
1051                 triangle_id=None):
1052        self.segment = segment
1053        self.normal = normal
1054        self.length = length
1055        self.triangle_id = triangle_id
1056
1057
1058    def __repr__(self):
1059        s = 'Triangle_intersection('
1060        s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\
1061             %(self.segment,
1062               self.normal,
1063               self.length,
1064               self.triangle_id)
1065
1066        return s
1067
1068
1069
1070def _get_intersecting_segments(V, N, line,
1071                               verbose=False):
1072    """Find edges intersected by line
1073
1074    Input:
1075        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1076        N: Number of triangles in mesh
1077        line - list of two points forming a segmented line
1078        verbose
1079    Output:
1080        list of instances of class Triangle_intersection
1081
1082    This method is used by the public method
1083    get_intersecting_segments(self, polyline) which also contains
1084    more documentation.
1085    """
1086
1087    from anuga.geometry.polygon import intersection
1088    from anuga.geometry.polygon import is_inside_polygon
1089
1090    msg = 'Line segment must contain exactly two points'
1091    assert len(line) == 2, msg
1092
1093    # Origin of intersecting line to be used for
1094    # establishing direction
1095    xi0 = line[0][0]
1096    eta0 = line[0][1]
1097
1098
1099    # Check intersection with edge segments for all triangles
1100    # FIXME (Ole): This should be implemented in C
1101    triangle_intersections={} # Keep track of segments already done
1102    for i in range(N):
1103        # Get nodes and edge segments for each triangle
1104        x0, y0 = V[3*i, :]
1105        x1, y1 = V[3*i+1, :]
1106        x2, y2 = V[3*i+2, :]
1107
1108
1109        edge_segments = [[[x0,y0], [x1, y1]],
1110                          [[x1,y1], [x2, y2]],
1111                          [[x2,y2], [x0, y0]]]
1112
1113        # Find segments that are intersected by line
1114
1115        intersections = {} # Use dictionary to record points only once
1116        for edge in edge_segments:
1117
1118            status, value = intersection(line, edge)
1119            #if value is not None: log.critical('Triangle %d, status=%s, '
1120            #                                   'value=%s'
1121            #                                   % (i, str(status), str(value)))
1122
1123            if status == 1:
1124                # Normal intersection of one edge or vertex
1125                intersections[tuple(value)] = i
1126
1127                # Exclude singular intersections with vertices
1128                #if not(allclose(value, edge[0]) or\
1129                #       allclose(value, edge[1])):
1130                #    intersections.append(value)
1131
1132            if status == 2:
1133                # Edge is sharing a segment with line
1134
1135                # This is usually covered by the two
1136                # vertices that would have been picked up
1137                # under status == 1.
1138                # However, if coinciding line stops partway
1139                # along this edge, it will be recorded here.
1140                intersections[tuple(value[0,:])] = i
1141                intersections[tuple(value[1,:])] = i
1142
1143
1144        if len(intersections) == 1:
1145            # Check if either line end point lies fully within this triangle
1146            # If this is the case accept that as one end of the intersecting
1147            # segment
1148
1149            poly = V[3*i:3*i+3]
1150            if is_inside_polygon(line[1], poly, closed=False):
1151                intersections[tuple(line[1])] = i
1152            elif is_inside_polygon(line[0], poly, closed=False):
1153                intersections[tuple(line[0])] = i
1154            else:
1155                # Ignore situations where one vertex is touch, for instance
1156                continue
1157
1158
1159        msg = 'There can be only two or no intersections'
1160        assert len(intersections) in [0,2], msg
1161
1162
1163        if len(intersections) == 2:
1164
1165            # Calculate attributes for this segment
1166
1167
1168            # End points of intersecting segment
1169            points = intersections.keys()
1170            x0, y0 = points[0]
1171            x1, y1 = points[1]
1172
1173
1174            # Determine which end point is closer to the origin of the line
1175            # This is necessary for determining the direction of
1176            # the line and the normals
1177
1178            # Distances from line origin to the two intersections
1179            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
1180            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
1181            d0 = num.sqrt(num.sum(z0**2))
1182            d1 = num.sqrt(num.sum(z1**2))
1183
1184            if d1 < d0:
1185                # Swap
1186                xi, eta = x0, y0
1187                x0, y0 = x1, y1
1188                x1, y1 = xi, eta
1189
1190            # (x0,y0) is now the origin of the intersecting segment
1191
1192
1193            # Normal direction:
1194            # Right hand side relative to line direction
1195            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
1196            length = num.sqrt(num.sum(vector**2))      # Segment length
1197            normal = num.array([vector[1], -vector[0]], num.float)/length
1198
1199
1200            segment = ((x0,y0), (x1, y1))
1201            T = Triangle_intersection(segment=segment,
1202                                      normal=normal,
1203                                      length=length,
1204                                      triangle_id=i)
1205
1206
1207            # Add segment unless it was done earlier
1208            if not triangle_intersections.has_key(segment):
1209                triangle_intersections[segment] = T
1210
1211
1212    # Return segments as a list
1213    return triangle_intersections.values()
1214
1215
1216def get_intersecting_segments(V, N, polyline,
1217                              verbose=False):
1218    """Internal function to find edges intersected by Polyline
1219
1220    Input:
1221        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
1222        N: Number of triangles in mesh
1223        polyline - list of points forming a segmented line
1224        verbose
1225    Output:
1226        list of instances of class Triangle_intersection
1227
1228    This method is used by the public method
1229    get_intersecting_segments(self, polyline) which also contains
1230    more documentation.
1231    """
1232
1233    msg = 'Polyline must contain at least two points'
1234    assert len(polyline) >= 2, msg
1235
1236
1237    # For all segments in polyline
1238    triangle_intersections = []
1239    for i, point0 in enumerate(polyline[:-1]):
1240
1241        point1 = polyline[i+1]
1242        if verbose:
1243            log.critical('Extracting mesh intersections from line:')
1244            log.critical('(%.2f, %.2f) - (%.2f, %.2f)'
1245                         % (point0[0], point0[1], point1[0], point1[1]))
1246
1247        line = [point0, point1]
1248        triangle_intersections += _get_intersecting_segments(V, N, line,
1249                                                             verbose=verbose)
1250
1251
1252    msg = 'No segments found'
1253    assert len(triangle_intersections) > 0, msg
1254
1255
1256    return triangle_intersections
1257
1258
1259
1260
1261
1262def segment_midpoints(segments):
1263    """Calculate midpoints of all segments
1264
1265    Inputs:
1266       segments: List of instances of class Segment
1267
1268    Ouputs:
1269       midpoints: List of points
1270    """
1271
1272    midpoints = []
1273    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
1274    for segment in segments:
1275        assert isinstance(segment, Triangle_intersection), msg
1276
1277        midpoint = num.sum(num.array(segment.segment, num.float), axis=0)/2
1278        midpoints.append(midpoint)
1279
1280    return midpoints
1281
1282
1283
Note: See TracBrowser for help on using the repository browser.