source: anuga_core/source/obsolete_code/obs_mesh.py @ 5571

Last change on this file since 5571 was 4906, checked in by duncan, 16 years ago

removing dead code
File size: 41.3 KB
Line 
1
2SET_COLOUR='red'
3
4class Triangle(MeshObject):
5    """
6    A triangle element, defined by 3 vertices.
7    Attributes based on the Triangle program.
8    """
9
10    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
11        """
12        Vertices, the initial arguments, are listed in counterclockwise order.
13        """
14        self.vertices= [vertex1,vertex2, vertex3 ]
15       
16        if attribute is None:
17            self.attribute =""
18        else:
19            self.attribute = attribute #this is a string
20           
21        if neighbors is None:
22            self.neighbors=[]
23        else:
24            self.neighbors=neighbors
25
26    def replace(self,new_triangle):
27        self = new_triangle
28
29    def longestSideID(self):
30        ax = self.vertices[0].x
31        ay = self.vertices[0].y
32       
33        bx = self.vertices[1].x
34        by = self.vertices[1].y
35       
36        cx = self.vertices[2].x
37        cy = self.vertices[2].y
38
39        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
40        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
41        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
42 
43        len = [lenA,lenB,lenC]
44        return len.index(max(len))
45
46    def rotate(self,offset):
47        """
48        permute the order of the sides of the triangle
49        offset must be 0,1 or 2
50        """
51
52        if offset == 0:
53            pass
54        else:
55            if offset == 1:
56                self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
57                self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
58            if offset == 2:
59                self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
60                self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
61
62    def rotate_longest_side(self):
63        self.rotate(self.longestSideID())
64
65    def getVertices(self):
66        return self.vertices
67   
68    def get_vertices(self):
69        """
70        Return a list of the vertices.  The x and y values will be relative
71        Easting and Northings for the zone of the current geo_ref.
72        """
73        return self.vertices
74   
75    def calcArea(self):
76        ax = self.vertices[0].x
77        ay = self.vertices[0].y
78       
79        bx = self.vertices[1].x
80        by = self.vertices[1].y
81       
82        cx = self.vertices[2].x
83        cy = self.vertices[2].y
84       
85        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
86   
87    def calcP(self):
88        #calculate the perimeter
89        ax = self.vertices[0].x
90        ay = self.vertices[0].y
91       
92        bx = self.vertices[1].x
93        by = self.vertices[1].y
94       
95        cx = self.vertices[2].x
96        cy = self.vertices[2].y
97
98        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
99        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
100        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
101
102        return a+b+c
103           
104    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
105        """
106        neighbor1 is the triangle opposite vertex1 and so on.
107        Null represents no neighbor
108        """
109        self.neighbors = [neighbor1, neighbor2, neighbor3]
110
111    def setAttribute(self,attribute):
112        """
113        neighbor1 is the triangle opposite vertex1 and so on.
114        Null represents no neighbor
115        """
116        self.attribute = attribute #this is a string
117       
118    def __repr__(self):
119        return "[%s,%s]" % (self.vertices,self.attribute)
120       
121
122    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
123        """
124        Draw a triangle, returning the objectID
125        """
126        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
127                                     scale*-1*(self.vertices[1].y + yoffset),
128                                     scale*(self.vertices[0].x + xoffset),
129                                     scale*-1*(self.vertices[0].y + yoffset),
130                                     scale*(self.vertices[2].x + xoffset),
131                                     scale*-1*(self.vertices[2].y + yoffset),
132                                     tags = tags,
133                                     outline = colour,fill = '')
134
135
136        self.setID={}
137        #a dictionary of names.
138        #multiple sets are allowed, but the gui does not yet
139        #support this
140       
141        self.setID['None']=0
142        #contains the names of the sets pointing to the indexes
143        #in the list.
144       
145        self.sets=[[]]
146        #Contains the lists of triangles (triangle sets)
147
148      ##################
149
150       
151    def refineSet(self,setName):
152        Triangles = self.sets[self.setID[setName]]
153        Refine(self,Triangles)
154
155    def selectAllTriangles(self):
156        A=[]
157        A.extend(self.meshTriangles)
158        if not('All' in self.setID.keys()):
159            self.setID['All']=len(self.sets)
160            self.sets.append(A)
161        else:
162            self.sets[self.setID['All']]=A
163        return 'All'
164        # and objectIDs
165
166
167    def clearSelection(self):
168        A = []
169        if not('None' in self.setID.keys()):
170            self.setID['None']=len(self.sets)
171            self.sets.append(A)
172        return 'None'
173
174    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
175    #FIXME Draws over previous triangles - may bloat canvas
176        Triangles = self.sets[self.setID[setName]]
177        for triangle in Triangles:
178            triangle.draw(canvas,1,
179                          scale = SCALE,
180                          colour = colour)
181           
182    def undrawSet(self,canvas,setName,SCALE,colour='green'):
183    #FIXME Draws over previous lines - may bloat canvas
184        Triangles = self.sets[self.setID[setName]]
185        for triangle in Triangles:
186            triangle.draw(canvas,1,
187                          scale = SCALE,
188                          colour = colour)
189
190    def weed(self,Vertices,Segments):
191        #Depreciated
192        #weed out existing duplicates
193        print 'len(self.getUserSegments())'
194        print len(self.getUserSegments())
195        print 'len(self.getUserVertices())'
196        print len(self.getUserVertices())
197
198        point_keys = {}
199        for vertex in Vertices:
200            point = (vertex.x,vertex.y)
201            point_keys[point]=vertex
202        #inlined would looks very ugly
203
204        line_keys = {}
205        for segment in Segments:
206            vertex1 = segment.vertices[0]
207            vertex2 = segment.vertices[1]
208            point1 = (vertex1.x,vertex1.y)
209            point2 = (vertex2.x,vertex2.y)
210            segment.vertices[0]=point_keys[point1]
211            segment.vertices[1]=point_keys[point2]
212            vertex1 = segment.vertices[0]
213            vertex2 = segment.vertices[1]
214            point1 = (vertex1.x,vertex1.y)
215            point2 = (vertex2.x,vertex2.y)
216            line1 = (point1,point2)
217            line2 = (point2,point1)
218            if not (line_keys.has_key(line1) \
219                 or line_keys.has_key(line2)):
220                 line_keys[line1]=segment
221        Vertices=point_keys.values()
222        Segments=line_keys.values()
223        return Vertices,Segments
224
225    def segs_to_dict(self,segments):
226        dict={}
227        for segment in segments:
228            vertex1 = segment.vertices[0]
229            vertex2 = segment.vertices[1]
230            point1 = (vertex1.x,vertex1.y)
231            point2 = (vertex2.x,vertex2.y)
232            line = (point1,point2)
233            dict[line]=segment
234        return dict
235
236    def seg2line(self,s):
237        return ((s.vertices[0].x,s.vertices[0].y,)\
238                (s.vertices[1].x,s.vertices[1].y))
239
240    def line2seg(self,line,tag=None):
241        point0 = self.point2ver(line[0])
242        point1 = self.point2ver(line[1])
243        return Segment(point0,point1,tag=tag)
244
245    def ver2point(self,vertex):
246        return (vertex.x,vertex.y)
247
248    def point2ver(self,point):
249        return Vertex(point[0],point[1])
250
251    def smooth_polySet(self,min_radius=0.05):
252        #for all pairs of connecting segments:
253        #    propose a new segment that replaces the 2
254
255        #    If the difference between the new segment
256        #    and the old lines is small: replace the
257        #    old lines.
258
259        seg2line = self.seg2line
260        ver2point= self.ver2point
261        line2seg = self.line2seg
262        point2ver= self.point2ver
263
264        #create dictionaries of lines -> segments
265        userSegments = self.segs_to_dict(self.userSegments)
266        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
267
268        #lump user and alpha segments
269        for key in alphaSegments.keys():
270            userSegments[key]=alphaSegments[key]
271
272        #point_keys = tuple -> vertex
273        #userVertices = vertex -> [line,line] - lines from that node
274        point_keys = {}
275        userVertices={}
276        for vertex in self.getUserVertices():
277            point = ver2point(vertex)
278            if not point_keys.has_key(point):
279                point_keys[point]=vertex
280                userVertices[vertex]=[]
281        for key in userSegments.keys():
282            line = key
283            point_0 = key[0]
284            point_1 = key[1]
285            userVertices[point_keys[point_0]].append(line)
286            userVertices[point_keys[point_1]].append(line)
287
288        for point in point_keys.keys():
289            try:
290            #removed keys can cause keyerrors
291                vertex = point_keys[point]
292                lines = userVertices[vertex]
293   
294                #if there are 2 lines on the node
295                if len(lines)==2:
296                    line_0 = lines[0]
297                    line_1 = lines[1]
298   
299                    #if the tags are the the same on the 2 lines
300                    if userSegments[line_0].tag == userSegments[line_1].tag:
301                        tag = userSegments[line_0].tag
302   
303                        #point_a is one of the next nodes, point_b is the other
304                        if point==line_0[0]:
305                            point_a = line_0[1]
306                        if point==line_0[1]:
307                            point_a = line_0[0]
308                        if point==line_1[0]:
309                            point_b = line_1[1]
310                        if point==line_1[1]:
311                            point_b = line_1[0]
312   
313   
314                        #line_2 is proposed
315                        line_2 = (point_a,point_b)
316
317                        #calculate the area of the triangle between
318                        #the two existing segments and the proposed
319                        #new segment
320                        ax = point_a[0]
321                        ay = point_a[1]
322                        bx = point_b[0]
323                        by = point_b[1]
324                        cx = point[0]
325                        cy = point[1]
326                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
327
328                        #calculate the perimeter
329                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
330                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
331                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
332                        perimeter = len_a+len_b+len_c
333
334                        #calculate the radius
335                        r = area/(2*perimeter)
336
337                        #if the radius is small: then replace the existing
338                        #segments with the new one
339                        if r < min_radius:
340                            if len_c < min_radius: append = False
341                            else: append = True
342                            #if the new seg is also time, don't add it
343                            if append:
344                                segment = self.line2seg(line_2,tag=tag)
345
346                            list_a=userVertices[point_keys[point_a]]
347                            list_b=userVertices[point_keys[point_b]]
348
349                            if line_0 in list_a:
350                                list_a.remove(line_0)
351                            else:
352                                list_a.remove(line_1)
353
354                            if line_0 in list_b:
355                                list_b.remove(line_0)
356                            else:
357                                list_b.remove(line_1)
358
359                            if append:
360                                list_a.append(line_2)
361                                list_b.append(line_2)
362                            else:
363                                if len(list_a)==0:
364                                    userVertices.pop(point_keys[point_a])
365                                    point_keys.pop(point_a)
366                                if len(list_b)==0:
367                                    userVertices.pop(point_keys[point_b])
368                                    point_keys.pop(point_b)
369
370                            userVertices.pop(point_keys[point])
371                            point_keys.pop(point)
372                            userSegments.pop(line_0)
373                            userSegments.pop(line_1)
374
375                            if append:
376                                userSegments[line_2]=segment
377            except:
378                pass
379
380        #self.userVerticies = userVertices.keys()
381        #self.userSegments = []
382        #for key in userSegments.keys():
383        #    self.userSegments.append(userSegments[key])
384        #self.alphaUserSegments = []
385
386        self.userVerticies = []
387        self.userSegments = []
388        self.alphaUserSegments = []
389
390        return userVertices,userSegments,alphaSegments
391
392    def triangles_to_polySet(self,setName):
393        #self.smooth_polySet()
394
395        seg2line = self.seg2line
396        ver2point= self.ver2point
397        line2seg = self.line2seg
398        point2ver= self.point2ver
399
400        from Numeric import array,allclose
401        #turn the triangles into a set
402        Triangles = self.sets[self.setID[setName]]
403        Triangles_dict = {}
404        for triangle in Triangles:
405            Triangles_dict[triangle]=None
406 
407
408        #create a dict of points to vertexes (tuple -> object)
409        #also create a set of vertexes (object -> True)
410        point_keys = {}
411        userVertices={}
412        for vertex in self.getUserVertices():
413            point = ver2point(vertex)
414            if not point_keys.has_key(point):
415                point_keys[point]=vertex
416                userVertices[vertex]=True
417
418        #create a dict of lines to segments (tuple -> object)
419        userSegments = self.segs_to_dict(self.userSegments)
420        #append the userlines in an affine linespace
421        affine_lines = Affine_Linespace()
422        for line in userSegments.keys():
423            affine_lines.append(line)
424        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
425        for line in alphaSegments.keys():
426            affine_lines.append(line)
427       
428        for triangle in Triangles:
429            for i in (0,1,2):
430                #for every triangles neighbour:
431                if not Triangles_dict.has_key(triangle.neighbors[i]):
432                #if the neighbour is not in the set:
433                    a = triangle.vertices[i-1]
434                    b = triangle.vertices[i-2]
435                    #Get possible matches:
436                    point_a = ver2point(a)
437                    point_b = ver2point(b)
438                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
439                    line = (point_a,point_b)
440                    tag = None
441
442
443                    #this bit checks for matching lines
444                    possible_lines = affine_lines[line] 
445                    possible_lines = unique(possible_lines)
446                    found = 0                           
447                    for user_line in possible_lines:
448                        if self.point_on_line(midpoint,user_line):
449                            found+=1
450                            assert found<2
451                            if userSegments.has_key(user_line):
452                                parent_segment = userSegments.pop(user_line)
453                            if alphaSegments.has_key(user_line):
454                                parent_segment = alphaSegments.pop(user_line)
455                            tag = parent_segment.tag
456                            offspring = [line]
457                            offspring.extend(self.subtract_line(user_line,
458                                                                line))
459                            affine_lines.remove(user_line)
460                            for newline in offspring:
461                                line_vertices = []
462                                for point in newline:
463                                    if point_keys.has_key(point):
464                                        vert = point_keys[point]
465                                    else:
466                                        vert = Vertex(point[0],point[1])
467                                        userVertices[vert]=True
468                                        point_keys[point]=vert
469                                    line_vertices.append(vert)
470                                segment = Segment(line_vertices[0],
471                                                  line_vertices[1],tag)
472                                userSegments[newline]=segment
473                                affine_lines.append(newline)
474                            #break
475                    assert found<2
476
477
478
479                    #if no matching lines
480                    if not found:
481                        line_vertices = []
482                        for point in line:
483                            if point_keys.has_key(point):
484                                vert = point_keys[point]
485                            else:
486                                vert = Vertex(point[0],point[1])
487                                userVertices[vert]=True
488                                point_keys[point]=vert
489                            line_vertices.append(vert)
490                        segment = Segment(line_vertices[0],
491                                          line_vertices[1],tag)
492                        userSegments[line]=segment
493                        affine_lines.append(line)
494       
495        self.userVerticies = []
496        self.userSegments = []
497        self.alphaUserSegments = []
498
499        return userVertices,userSegments,alphaSegments
500
501    def subtract_line(self,parent,child):
502        #Subtracts child from parent
503        #Requires that the child is a
504        #subline of parent to work.
505
506        from Numeric import allclose,dot,array
507        A= parent[0]
508        B= parent[1]
509        a = child[0]
510        b = child[1]
511
512        A_array = array(parent[0])
513        B_array = array(parent[1])
514        a_array   = array(child[0])
515        b_array   = array(child[1])
516
517        assert not A == B
518        assert not a == b
519
520        answer = []
521
522        #if the new line does not share a
523        #vertex with the old one
524        if not (allclose(A_array,a_array)\
525             or allclose(B_array,b_array)\
526             or allclose(A_array,b_array)\
527             or allclose(a_array,B_array)):
528            if dot(A_array-a_array,A_array-a_array) \
529            < dot(A_array-b_array,A_array-b_array):
530                sibling1 = (A,a)
531                sibling2 = (B,b)
532                return [sibling1,sibling2]
533            else:
534                sibling1 = (A,b)
535                sibling2 = (B,a)
536                return [sibling1,sibling2]
537
538        elif allclose(A_array,a_array):
539            if allclose(B_array,b_array):
540                return []
541            else:
542                sibling = (b,B)
543                return [sibling]
544        elif allclose(B_array,b_array):
545            sibling = (a,A)
546            return [sibling]
547
548        elif allclose(A_array,b_array):
549            if allclose(B,a):
550                return []
551            else:
552                sibling = (a,B)
553                return [sibling]
554        elif allclose(a_array,B_array):
555            sibling = (b,A)
556            return [sibling]
557
558    def point_on_line(self,point,line):
559        #returns true within a tolerance of 3 degrees
560        x=point[0]
561        y=point[1]
562        x0=line[0][0]
563        x1=line[1][0]
564        y0=line[0][1]
565        y1=line[1][1]
566        from Numeric import array, dot, allclose
567        from math import sqrt
568        tol = 3. #DEGREES
569        tol = tol*3.1415/180
570
571        a = array([x - x0, y - y0]) 
572        a_normal = array([a[1], -a[0]])     
573        len_a_normal = sqrt(sum(a_normal**2)) 
574
575        b = array([x1 - x0, y1 - y0])         
576        len_b = sqrt(sum(b**2)) 
577   
578        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
579            #Point is somewhere on the infinite extension of the line
580
581            len_a = sqrt(sum(a**2))                                     
582            if dot(a, b) >= 0 and len_a <= len_b:
583               return True
584            else:   
585               return False
586        else:
587          return False
588
589    def line_length(self,line):
590        x0=line[0][0]
591        x1=line[1][0]
592        y0=line[0][1]
593        y1=line[1][1]
594        return ((x1-x0)**2-(y1-y0)**2)**0.5     
595
596    def threshold(self,setName,min=None,max=None,attribute_name='elevation'):
597        """
598        threshold using  d
599        """
600        triangles = self.sets[self.setID[setName]]
601        A = []
602
603        if attribute_name in self.attributeTitles:
604            i = self.attributeTitles.index(attribute_name)
605        else: i = -1#no attribute
606        if not max == None:
607            for t in triangles:
608                if (min<self.av_att(t,i)<max):
609                    A.append(t)
610        else:
611            for t in triangles:
612                if (min<self.av_att(t,i)):
613                    A.append(t)
614        self.sets[self.setID[setName]] = A
615
616    def general_threshold(self,setName,min=None,max=None\
617              ,attribute_name = 'elevation',function=None):
618        """
619        Thresholds the triangles
620        """
621        from visual.graph import arange,ghistogram,color as colour
622        triangles = self.sets[self.setID[setName]]
623        A = []
624        data=[]
625        #data is for the graph
626
627        if attribute_name in self.attributeTitles:
628            i = self.attributeTitles.index(attribute_name)
629        else: i = -1
630        if not max == None:
631            for t in triangles:
632                value=function(t,i)
633                if (min<value<max):
634                    A.append(t)
635                data.append(value)
636        else:
637            for t in triangles:
638                value=function(t,i)
639                if (min<value):
640                    A.append(t)
641                data.append(value)
642        self.sets[self.setID[setName]] = A
643
644        if self.visualise_graph:
645            if len(data)>0:
646                max=data[0]
647                min=data[0]
648                for value in data:
649                    if value > max:
650                        max = value
651                    if value < min:
652                        min = value
653
654                inc = (max-min)/100
655
656                histogram = ghistogram(bins=arange(min,max,inc),\
657                             color = colour.red)
658                histogram.plot(data=data)
659       
660    def av_att(self,triangle,i):
661        if i==-1: return 1
662        else:
663            #evaluates the average attribute of the vertices of a triangle.
664            V = triangle.getVertices()
665            a0 = (V[0].attributes[i])
666            a1 = (V[1].attributes[i])
667            a2 = (V[2].attributes[i])
668            return (a0+a1+a2)/3
669
670    def Courant_ratio(self,triangle,index):
671        """
672        Uses the courant threshold
673        """
674        e = self.av_att(triangle,index)
675        A = triangle.calcArea()
676        P = triangle.calcP()
677        r = A/(2*P)
678        e = max(0.1,abs(e))
679        return r/e**0.5
680
681    def Gradient(self,triangle,index):
682        V = triangle.vertices
683        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
684        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
685        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
686            print ((grad_x**2)+(grad_y**2))**(0.5)
687        return ((grad_x**2)+(grad_y**2))**(0.5)
688   
689
690    def append_triangle(self,triangle):
691        self.meshTriangles.append(triangle)
692
693    def replace_triangle(self,triangle,replacement):
694        i = self.meshTriangles.index(triangle)
695        self.meshTriangles[i]=replacement
696        assert replacement in self.meshTriangles
697
698 
699"""Refines triangles
700
701   Implements the #triangular bisection?# algorithm.
702 
703   
704"""
705
706def Refine(mesh, triangles):
707    """
708    Given a general_mesh, and a triangle number, split
709    that triangle in the mesh in half. Then to prevent
710    vertices and edges from meeting, keep refining
711    neighbouring triangles until the mesh is clean.
712    """
713    state = BisectionState(mesh)
714    for triangle in triangles:
715        if not state.refined_triangles.has_key(triangle):
716            triangle.rotate_longest_side()
717            state.start(triangle)
718            Refine_mesh(mesh, state)
719
720def Refine_mesh(mesh, state):
721    """
722    """
723    state.getState(mesh)
724    refine_triangle(mesh,state)
725    state.evolve()
726    if not state.end:
727        Refine_mesh(mesh,state)
728
729def refine_triangle(mesh,state):
730    split(mesh,state.current_triangle,state.new_point)
731    if state.case == 'one':
732        state.r[3]=state.current_triangle#triangle 2
733
734        new_triangle_id = len(mesh.meshTriangles)-1
735        new_triangle = mesh.meshTriangles[new_triangle_id]
736
737        split(mesh,new_triangle,state.old_point)
738        state.r[2]=new_triangle#triangle 1.2
739        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
740        r = state.r
741        state.repairCaseOne()
742
743    if state.case == 'two':
744        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
745
746        new_triangle = state.current_triangle
747
748        split(mesh,new_triangle,state.old_point)
749
750        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
751        state.r[4]=new_triangle#triangle 2.2
752        r = state.r
753
754        state.repairCaseTwo()
755
756    if state.case == 'vertex':
757        state.r[2]=state.current_triangle#triangle 2
758        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
759        r = state.r
760        state.repairCaseVertex()
761       
762    if state.case == 'start':
763        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
764        state.r[3]=state.current_triangle#triangle 2
765
766    if state.next_case == 'boundary':
767        state.repairCaseBoundary()
768
769
770def split(mesh, triangle, new_point):
771        """
772        Given a mesh, triangle_id and a new point,
773        split the corrosponding triangle into two
774        new triangles and update the mesh.
775        """
776
777        new_triangle1 = Triangle(new_point,triangle.vertices[0],
778                                 triangle.vertices[1],
779                                 attribute = triangle.attribute,
780                                 neighbors = None)
781        new_triangle2 = Triangle(new_point,triangle.vertices[2],
782                                 triangle.vertices[0],
783                                 attribute = triangle.attribute,
784                                 neighbors = None)
785
786        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
787        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
788
789        mesh.meshTriangles.append(new_triangle1)
790
791        triangle.vertices = new_triangle2.vertices
792        triangle.neighbors = new_triangle2.neighbors
793
794
795class State:
796
797    def __init__(self):
798        pass
799
800class BisectionState(State):
801
802
803    def __init__(self,mesh):
804        self.len = len(mesh.meshTriangles)
805        self.refined_triangles = {}
806        self.mesh = mesh
807        self.current_triangle = None
808        self.case = 'start'
809        self.end = False
810        self.r = [None,None,None,None,None]
811
812    def start(self, triangle):
813        self.current_triangle = triangle
814        self.case = 'start'
815        self.end = False
816        self.r = [None,None,None,None,None]
817
818    def getState(self,mesh):
819        if not self.case == 'vertex':
820            self.new_point=self.getNewVertex(mesh, self.current_triangle)
821            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
822            self.neighbour = self.current_triangle.neighbors[0]
823            if not self.neighbour is None:
824                self.neighbour.rotate_longest_side()
825            self.next_case = self.get_next_case(mesh,self.neighbour,
826                                                self.current_triangle)
827        if self.case == 'vertex':
828            self.new_point=self.old_point
829
830
831    def evolve(self):
832        if self.case == 'vertex':
833            self.end = True
834
835        self.last_case = self.case
836        self.case = self.next_case
837
838        self.old_point = self.new_point
839        self.current_triangle = self.neighbour
840
841        if self.case == 'boundary':
842            self.end = True
843        self.refined_triangles[self.r[2]]=1
844        self.refined_triangles[self.r[3]]=1
845        if not self.r[4] is None:
846            self.refined_triangles[self.r[4]]=1
847        self.r[0]=self.r[2]
848        self.r[1]=self.r[3]
849
850
851    def getNewVertex(self,mesh,triangle):
852        coordinate1 = triangle.vertices[1]
853        coordinate2 = triangle.vertices[2]
854        a = ([coordinate1.x*1.,coordinate1.y*1.])
855        b = ([coordinate2.x*1.,coordinate2.y*1.])
856        attributes = []
857        for i in range(len(coordinate1.attributes)):
858            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
859            attributes.append(att)
860        new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
861        newVertex = Vertex(new_coordinate[0],new_coordinate[1],
862                           attributes = attributes)
863        mesh.maxVertexIndex+=1
864        newVertex.index = mesh.maxVertexIndex
865        mesh.meshVertices.append(newVertex)
866        return newVertex
867
868    def get_next_case(self, mesh,neighbour,triangle):
869        """
870        Given the locations of two neighbouring triangles,
871        examine the interior indices of their vertices (i.e.
872        0,1 or 2) to determine what how the neighbour needs
873        to be refined.
874        """
875        if (neighbour is None):
876                next_case = 'boundary'
877        else:
878                if triangle.vertices[1].x==neighbour.vertices[2].x:
879                    if triangle.vertices[1].y==neighbour.vertices[2].y:
880                        next_case = 'vertex'
881                if triangle.vertices[1].x==neighbour.vertices[0].x:
882                    if triangle.vertices[1].y==neighbour.vertices[0].y:
883                        next_case = 'two'
884                if triangle.vertices[1].x==neighbour.vertices[1].x:
885                    if triangle.vertices[1].y==neighbour.vertices[1].y:
886                        next_case = 'one'
887        return next_case
888
889
890
891    def repairCaseVertex(self):
892
893        r = self.r
894
895
896        self.link(r[0],r[2])
897        self.repair(r[0])
898
899        self.link(r[1],r[3])
900        self.repair(r[1])
901
902        self.repair(r[2])
903
904        self.repair(r[3])
905
906
907    def repairCaseOne(self):
908        r = self.rkey
909
910
911        self.link(r[0],r[2])
912        self.repair(r[0])
913
914        self.link(r[1],r[4])
915        self.repair(r[1])
916
917        self.repair(r[4])
918
919    def repairCaseTwo(self):
920        r = self.r
921
922        self.link(r[0],r[4])
923        self.repair(r[0])
924
925        self.link(r[1],r[3])
926        self.repair(r[1])
927
928        self.repair(r[4])
929
930    def repairCaseBoundary(self):
931        r = self.r
932        self.repair(r[2])
933        self.repair(r[3])
934
935
936
937    def repair(self,triangle):
938        """
939        Given a triangle that knows its neighbours, this will
940        force the neighbours to comply.
941
942        However, it needs to compare the vertices of triangles
943        for this implementation
944
945        But it doesn't work for invalid neighbour structures
946        """
947        n=triangle.neighbors
948        for i in (0,1,2):
949            if not n[i] is None:
950                for j in (0,1,2):#to find which side of the list is broken
951                    if not (n[i].vertices[j] in triangle.vertices):
952                    #ie if j is the side of n that needs fixing
953                        k = j
954                n[i].neighbors[k]=triangle
955
956
957
958    def link(self,triangle1,triangle2):
959        """
960        make triangle1 neighbors point to t
961                #count = 0riangle2
962        """
963        count = 0
964        for i in (0,1,2):#to find which side of the list is broken
965            if not (triangle1.vertices[i] in triangle2.vertices):
966                j = i
967                count+=1
968        assert count == 1
969        triangle1.neighbors[j]=triangle2
970
971class Discretised_Tuple_Set:
972    """
973    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
974    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
975    a[(10000)]=[] #NOT KEYERROR
976
977    a.append[(0.01)]
978    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
979
980    #NOT IMPLIMENTED
981    a.remove[(0.01)]
982    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
983
984    a.remove[(0.17)]
985    => {(0.0):[(0.02),(0.01)],0.2:[]}
986    #NOT IMPLIMENTED
987    at a.precision = 2:
988        a.round_up_rel[0.0]=
989        a.round_flat[0.0]=
990        a.round_down_rel[0.0]=
991
992        a.up((0.1,2.04))=
993
994    If t_rel = 0, nothing gets rounded into
995    two bins. If t_rel = 0.5, everything does.
996
997    Ideally, precision can be set high, so that multiple
998    entries are rarely in the same bin. And t_rel should
999    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
1000    so that it is rare to put items in mutiple bins.
1001
1002    Ex bins per entry = product(a,b,c...,n)
1003    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
1004    b = 1 or 2 ...
1005
1006    BUT!!! to avoid missing the right bin:
1007    (-10)**(precision+1)*t_rel must be greater than the
1008    greatest possible variation that an identical element
1009    can display.
1010
1011
1012    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
1013    but not .6 - this looks wrong, but note that *everything* will round,
1014    so .6 wont be missed as everything close to it will check in .7 and .5.
1015    """
1016    def __init__(self,p_rel = 6,t_rel = 0.01):
1017        self.__p_rel__ = p_rel
1018        self.__t_rel__ = t_rel
1019
1020        self.__p_abs__ = p_rel+1
1021        self.__t_abs__ = t_rel
1022
1023        assert t_rel <= 0.5
1024        self.__items__ = {}
1025        from math import frexp
1026        self.frexp = frexp
1027        roundings = [self.round_up_rel,\
1028        self.round_down_rel,self.round_flat_rel,\
1029        self.round_down_abs,self.round_up_abs,\
1030        self.round_flat_abs]#
1031
1032        self.roundings = roundings
1033
1034    def __repr__(self):
1035        return '%s'%self.__items__
1036
1037    def rounded_keys(self,key):
1038        key = tuple(key)
1039        keys = [key]
1040        keys = self.__rounded_keys__(key)
1041        return (keys)
1042
1043    def __rounded_keys__(self,key):
1044        keys = []
1045        rounded_key=list(key)
1046        rounded_values=list(key)
1047
1048        roundings = list(self.roundings)
1049
1050        #initialise rounded_values
1051        round = roundings.pop(0)
1052        for i in range(len(rounded_values)):
1053            rounded_key[i]=round(key[i])
1054            rounded_values[i]={}
1055            rounded_values[i][rounded_key[i]]=None
1056        keys.append(tuple(rounded_key))
1057       
1058        for round in roundings:
1059            for i in range(len(rounded_key)):
1060                rounded_value=round(key[i])
1061                if not rounded_values[i].has_key(rounded_value):
1062                    #ie unless round_up_rel = round_down_rel
1063                    #so the keys stay unique
1064                    for j in range(len(keys)):
1065                        rounded_key = list(keys[j])
1066                        rounded_key[i]=rounded_value
1067                        keys.append(tuple(rounded_key))
1068        return keys
1069
1070    def append(self,item):
1071        keys = self.rounded_keys(item)
1072        for key in keys:
1073            if self.__items__.has_key(key):
1074                self.__items__[key].append(item)
1075            else:
1076                self.__items__[key]=[item]
1077
1078    def __getitem__(self,key):
1079        answer = []
1080        keys = self.rounded_keys(key)
1081        for key in keys:
1082            if self.__items__.has_key(key):
1083                answer.extend(self.__items__[key])
1084        #if len(answer)==0:
1085        #    raise KeyError#FIXME or return KeyError
1086        #                  #FIXME or just return []?
1087        else:
1088            return answer #FIXME or unique(answer)?
1089
1090    def __delete__(self,item):
1091        keys = self.rounded_keys(item)
1092        answer = False
1093        #if any of the possible keys contains
1094        #a list, return true
1095        for key in keys:       
1096            if self.__items__.has_key(key):
1097                if item in self.__items__[key]:
1098                    self.__items__[key].remove(item)
1099
1100    def remove(self,item):
1101        self.__delete__(item)
1102
1103    def __contains__(self,item):
1104
1105        keys = self.rounded_keys(item)
1106        answer = False
1107        #if any of the possible keys contains
1108        #a list, return true
1109        for key in keys:       
1110            if self.__items__.has_key(key):
1111                if item in self.__items__[key]:
1112                    answer = True
1113        return answer
1114
1115
1116    def has_item(self,item):
1117        return self.__contains__(item)
1118
1119    def round_up_rel2(self,value):
1120         t_rel=self.__t_rel__
1121         #Rounding up the value
1122         m,e = self.frexp(value)
1123         m = m/2
1124         e = e + 1
1125         #m is the mantissa, e the exponent
1126         # 0.5 < |m| < 1.0
1127         m = m+t_rel*(10**-(self.__p_rel__))
1128         #bump m up
1129         m = round(m,self.__p_rel__)
1130         return m*(2.**e)
1131
1132    def round_down_rel2(self,value):
1133         t_rel=self.__t_rel__
1134         #Rounding down the value
1135         m,e = self.frexp(value)
1136         m = m/2
1137         e = e + 1
1138         #m is the mantissa, e the exponent
1139         # 0.5 < m < 1.0
1140         m = m-t_rel*(10**-(self.__p_rel__))
1141         #bump the |m| down, by 5% or whatever
1142         #self.p_rel dictates
1143         m = round(m,self.__p_rel__)
1144         return m*(2.**e)
1145
1146    def round_flat_rel2(self,value):
1147    #redundant
1148         m,e = self.frexp(value)
1149         m = m/2
1150         e = e + 1
1151         m = round(m,self.__p_rel__)
1152         return m*(2.**e)
1153
1154    def round_up_rel(self,value):
1155         t_rel=self.__t_rel__
1156         #Rounding up the value
1157         m,e = self.frexp(value)
1158         #m is the mantissa, e the exponent
1159         # 0.5 < |m| < 1.0
1160         m = m+t_rel*(10**-(self.__p_rel__))
1161         #bump m up
1162         m = round(m,self.__p_rel__)
1163         return m*(2.**e)
1164
1165    def round_down_rel(self,value):
1166         t_rel=self.__t_rel__
1167         #Rounding down the value
1168         m,e = self.frexp(value)
1169         #m is the mantissa, e the exponent
1170         # 0.5 < m < 1.0
1171         m = m-t_rel*(10**-(self.__p_rel__))
1172         #bump the |m| down, by 5% or whatever
1173         #self.p_rel dictates
1174         m = round(m,self.__p_rel__)
1175         return m*(2.**e)
1176
1177    def round_flat_rel(self,value):
1178    #redundant
1179         m,e = self.frexp(value)
1180         m = round(m,self.__p_rel__)
1181         return m*(2.**e)
1182
1183    def round_up_abs(self,value):
1184         t_abs=self.__t_abs__
1185         #Rounding up the value
1186         m = value+t_abs*(10**-(self.__p_abs__))
1187         #bump m up
1188         m = round(m,self.__p_abs__)
1189         return m
1190
1191    def round_down_abs(self,value):
1192         t_abs=self.__t_abs__
1193         #Rounding down the value
1194         m = value-t_abs*(10**-(self.__p_abs__))
1195         #bump the |m| down, by 5% or whatever
1196         #self.p_rel dictates
1197         m = round(m,self.__p_abs__)
1198         return m
1199
1200    def round_flat_abs(self,value):
1201    #redundant?
1202         m = round(value,self.__p_abs__)
1203         return m
1204
1205    def keys(self):
1206        return self.__items__.keys()
1207
1208
1209class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
1210    """
1211    This is a discretised tuple set, but
1212    based on a mapping. The mapping MUST
1213    return a sequence.
1214
1215    example:
1216    def weight(animal):
1217        return [animal.weight]
1218   
1219    a = Mapped_Discretised_Tuple_Set(weight)
1220    a.append[cow]
1221    a.append[fox]
1222    a.append[horse]
1223
1224    a[horse]    -> [cow,horse]
1225    a[dog]      -> [fox]
1226    a[elephant] -> []
1227    """
1228    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
1229        Discretised_Tuple_Set.__init__\
1230         (self,p_rel,t_rel = t_rel)
1231        self.mapping = mapping
1232
1233    def rounded_keys(self,key):
1234        mapped_key = tuple(self.mapping(key))
1235        keys = self.__rounded_keys__(mapped_key)
1236        return keys
1237
1238class Affine_Linespace(Mapped_Discretised_Tuple_Set):
1239    """
1240    The affine linespace creates a way to record and compare lines.
1241    Precision is a bit of a hack, but it creates a way to avoid
1242    misses caused by round offs (between two lines of different
1243    lenghts, the short one gets rounded off more).
1244    I am starting to think that a quadratic search would be faster.
1245    Nearly.
1246    """
1247    def __init__(self,p_rel=4,t_rel=0.2):
1248        Mapped_Discretised_Tuple_Set.__init__\
1249            (self,self.affine_line,\
1250            p_rel=p_rel,t_rel=t_rel)
1251
1252        roundings = \
1253        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
1254        self.roundings = roundings
1255        #roundings = \
1256        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
1257        #self.roundings = roundings
1258
1259    def affine_line(self,line):
1260        point_1 = line[0]
1261        point_2 = line[1]
1262        #returns the equation of a line
1263        #between two points, in the from
1264        #(a,b,-c), as in ax+by-c=0
1265        #or line *dot* (x,y,1) = (0,0,0)
1266
1267        #Note that it normalises the line
1268        #(a,b,-c) so Line*Line = 1.
1269        #This does not change the mathematical
1270        #properties, but it makes comparism
1271        #easier.
1272
1273        #There are probably better algorithms.
1274        x1 = point_1[0]
1275        x2 = point_2[0]
1276        y1 = point_1[1]
1277        y2 = point_2[1]
1278        dif_x = x1-x2
1279        dif_y = y1-y2
1280
1281        if dif_x == dif_y == 0:
1282            msg = 'points are the same'
1283            raise msg
1284        elif abs(dif_x)>=abs(dif_y):
1285            alpha = (-dif_y)/dif_x
1286            #a = alpha * b
1287            b = -1.
1288            c = (x1*alpha+x2*alpha+y1+y2)/2.
1289            a = alpha*b
1290        else:
1291            beta = dif_x/(-dif_y)
1292            #b = beta * a
1293            a = 1.
1294            c = (x1+x2+y1*beta+y2*beta)/2.
1295            b = beta*a
1296        mag = abs(a)+abs(b)
1297        #This does not change the mathematical
1298        #properties, but it makes comparism possible.
1299
1300        #note that the gradient is b/a, or (a/b)**-1.
1301        #so
1302
1303        #if a == 0:
1304        #    sign_a = 1.
1305        #else:
1306        #    sign_a = a/((a**2)**0.5)
1307        #if b == 0:
1308        #    sign_b = 1.
1309        #else:
1310        #    sign_b = b/((b**2)**0.5)
1311        #if c == 0:
1312        #    sign_c = 1.
1313        #else:
1314        #    sign_c = c/((c**2)**0.5)
1315        #a = a/mag*sign_a
1316        #b = b/mag*sign_b
1317        #c = c/mag*sign_c
1318        a = a/mag
1319        b = b/mag
1320        c = c/mag
1321        return a,b,c
Note: See TracBrowser for help on using the repository browser.