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

Last change on this file since 8280 was 8280, checked in by neweyv, 12 years ago

Additional information written to the log file. This information will assist in creating an algorithm to predict the runtime of ANUGA given various parameters.

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