# source:anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/neighbour_mesh.py@5899

Last change on this file since 5899 was 5899, checked in by rwilson, 15 years ago

Initial NumPy? changes (again!).

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