Changeset 4906


Ignore:
Timestamp:
Jan 4, 2008, 2:41:09 PM (16 years ago)
Author:
duncan
Message:

removing dead code

Location:
anuga_core/source
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/pmesh/mesh.py

    r4905 r4906  
    5353    kinds = _kinds()
    5454   
    55 SET_COLOUR='red'
    56 
    5755#FIXME: this is not tested.
    5856from anuga.utilities.numerical_tools import gradient
     
    299297                                   self.getTag())
    300298       
    301 class Triangle(MeshObject):
    302     """
    303     A triangle element, defined by 3 vertices.
    304     Attributes based on the Triangle program.
    305     """
    306 
    307     def __init__(self, vertex1, vertex2, vertex3, attribute = None,
    308                  neighbors = None ):
    309         """
    310         Vertices, the initial arguments, are listed in counterclockwise order.
    311         """
    312         self.vertices= [vertex1,vertex2, vertex3 ]
    313        
    314         if attribute is None:
    315             self.attribute =""
    316         else:
    317             self.attribute = attribute #this is a string
    318            
    319         if neighbors is None:
    320             self.neighbors=[]
    321         else:
    322             self.neighbors=neighbors
    323 
    324     def replace(self,new_triangle):
    325         self = new_triangle
    326 
    327     def longestSideID(self):
    328         ax = self.vertices[0].x
    329         ay = self.vertices[0].y
    330        
    331         bx = self.vertices[1].x
    332         by = self.vertices[1].y
    333        
    334         cx = self.vertices[2].x
    335         cy = self.vertices[2].y
    336 
    337         lenA = ((cx-bx)**2+(cy-by)**2)**0.5
    338         lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
    339         lenC = ((bx-ax)**2+(by-ay)**2)**0.5
    340  
    341         len = [lenA,lenB,lenC]
    342         return len.index(max(len))
    343 
    344     def rotate(self,offset):
    345         """
    346         permute the order of the sides of the triangle
    347         offset must be 0,1 or 2
    348         """
    349 
    350         if offset == 0:
    351             pass
    352         else:
    353             if offset == 1:
    354                 self.vertices = [self.vertices[1],self.vertices[2],
    355                                  self.vertices[0]]
    356                 self.neighbors = [self.neighbors[1],self.neighbors[2],
    357                                   self.neighbors[0]]
    358             if offset == 2:
    359                 self.vertices = [self.vertices[2],self.vertices[0],
    360                                  self.vertices[1]]
    361                 self.neighbors = [self.neighbors[2],self.neighbors[0],
    362                                   self.neighbors[1]]
    363 
    364     def rotate_longest_side(self):
    365         self.rotate(self.longestSideID())
    366 
    367     def getVertices(self):
    368         return self.vertices
    369    
    370     def get_vertices(self):
    371         """
    372         Return a list of the vertices.  The x and y values will be relative
    373         Easting and Northings for the zone of the current geo_ref.
    374         """
    375         return self.vertices
    376    
    377     def calcArea(self):
    378         ax = self.vertices[0].x
    379         ay = self.vertices[0].y
    380        
    381         bx = self.vertices[1].x
    382         by = self.vertices[1].y
    383        
    384         cx = self.vertices[2].x
    385         cy = self.vertices[2].y
    386        
    387         return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
    388    
    389     def calcP(self):
    390         #calculate the perimeter
    391         ax = self.vertices[0].x
    392         ay = self.vertices[0].y
    393        
    394         bx = self.vertices[1].x
    395         by = self.vertices[1].y
    396        
    397         cx = self.vertices[2].x
    398         cy = self.vertices[2].y
    399 
    400         a =  ((cx-bx)**2+(cy-by)**2)**0.5
    401         b =  ((ax-cx)**2+(ay-cy)**2)**0.5
    402         c =  ((bx-ax)**2+(by-ay)**2)**0.5
    403 
    404         return a+b+c
    405            
    406     def setNeighbors(self,neighbor1=None, neighbor2=None, neighbor3=None):
    407         """
    408         neighbor1 is the triangle opposite vertex1 and so on.
    409         Null represents no neighbor
    410         """
    411         self.neighbors = [neighbor1, neighbor2, neighbor3]
    412 
    413     def setAttribute(self,attribute):
    414         """
    415         neighbor1 is the triangle opposite vertex1 and so on.
    416         Null represents no neighbor
    417         """
    418         self.attribute = attribute #this is a string
    419        
    420     def __repr__(self):
    421         return "[%s,%s]" % (self.vertices,self.attribute)
    422        
    423 
    424     def draw(self, canvas, tags, scale=1, xoffset=0, yoffset=0,
    425              colour="green"):
    426         """
    427         Draw a triangle, returning the objectID
    428         """
    429         return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
    430                                      scale*-1*(self.vertices[1].y + yoffset),
    431                                      scale*(self.vertices[0].x + xoffset),
    432                                      scale*-1*(self.vertices[0].y + yoffset),
    433                                      scale*(self.vertices[2].x + xoffset),
    434                                      scale*-1*(self.vertices[2].y + yoffset),
    435                                      tags = tags,
    436                                      outline = colour,fill = '')
    437 
    438299class Segment(MeshObject):
    439300    """
     
    631492        self.tri_mesh=None
    632493       
    633         self.setID={}
    634         #a dictionary of names.
    635         #multiple sets are allowed, but the gui does not yet
    636         #support this
    637        
    638         self.setID['None']=0
    639         #contains the names of the sets pointing to the indexes
    640         #in the list.
    641        
    642         self.sets=[[]]
    643         #Contains the lists of triangles (triangle sets)
    644 
    645        
     494       
    646495        self.visualise_graph = True
    647496
     
    23382187            self.regions.append(Object)
    23392188 
    2340 ############################################
    2341 
    2342        
    2343     def refineSet(self,setName):
    2344         Triangles = self.sets[self.setID[setName]]
    2345         Refine(self,Triangles)
    2346 
    2347     def selectAllTriangles(self):
    2348         A=[]
    2349         A.extend(self.meshTriangles)
    2350         if not('All' in self.setID.keys()):
    2351             self.setID['All']=len(self.sets)
    2352             self.sets.append(A)
    2353         else:
    2354             self.sets[self.setID['All']]=A
    2355         return 'All'
    2356         # and objectIDs
    2357 
    2358 
    2359     def clearSelection(self):
    2360         A = []
    2361         if not('None' in self.setID.keys()):
    2362             self.setID['None']=len(self.sets)
    2363             self.sets.append(A)
    2364         return 'None'
    2365 
    2366     def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
    2367     #FIXME Draws over previous triangles - may bloat canvas
    2368         Triangles = self.sets[self.setID[setName]]
    2369         for triangle in Triangles:
    2370             triangle.draw(canvas,1,
    2371                           scale = SCALE,
    2372                           colour = colour)
    2373            
    2374     def undrawSet(self,canvas,setName,SCALE,colour='green'):
    2375     #FIXME Draws over previous lines - may bloat canvas
    2376         Triangles = self.sets[self.setID[setName]]
    2377         for triangle in Triangles:
    2378             triangle.draw(canvas,1,
    2379                           scale = SCALE,
    2380                           colour = colour)
    2381 
    2382     def weed(self,Vertices,Segments):
    2383         #Depreciated
    2384         #weed out existing duplicates
    2385         print 'len(self.getUserSegments())'
    2386         print len(self.getUserSegments())
    2387         print 'len(self.getUserVertices())'
    2388         print len(self.getUserVertices())
    2389 
    2390         point_keys = {}
    2391         for vertex in Vertices:
    2392             point = (vertex.x,vertex.y)
    2393             point_keys[point]=vertex
    2394         #inlined would looks very ugly
    2395 
    2396         line_keys = {}
    2397         for segment in Segments:
    2398             vertex1 = segment.vertices[0]
    2399             vertex2 = segment.vertices[1]
    2400             point1 = (vertex1.x,vertex1.y)
    2401             point2 = (vertex2.x,vertex2.y)
    2402             segment.vertices[0]=point_keys[point1]
    2403             segment.vertices[1]=point_keys[point2]
    2404             vertex1 = segment.vertices[0]
    2405             vertex2 = segment.vertices[1]
    2406             point1 = (vertex1.x,vertex1.y)
    2407             point2 = (vertex2.x,vertex2.y)
    2408             line1 = (point1,point2)
    2409             line2 = (point2,point1)
    2410             if not (line_keys.has_key(line1) \
    2411                  or line_keys.has_key(line2)):
    2412                  line_keys[line1]=segment
    2413         Vertices=point_keys.values()
    2414         Segments=line_keys.values()
    2415         return Vertices,Segments
    2416 
    2417     def segs_to_dict(self,segments):
    2418         dict={}
    2419         for segment in segments:
    2420             vertex1 = segment.vertices[0]
    2421             vertex2 = segment.vertices[1]
    2422             point1 = (vertex1.x,vertex1.y)
    2423             point2 = (vertex2.x,vertex2.y)
    2424             line = (point1,point2)
    2425             dict[line]=segment
    2426         return dict
    2427 
    2428     def seg2line(self,s):
    2429         return ((s.vertices[0].x,s.vertices[0].y,)\
    2430                 (s.vertices[1].x,s.vertices[1].y))
    2431 
    2432     def line2seg(self,line,tag=None):
    2433         point0 = self.point2ver(line[0])
    2434         point1 = self.point2ver(line[1])
    2435         return Segment(point0,point1,tag=tag)
    2436 
    2437     def ver2point(self,vertex):
    2438         return (vertex.x,vertex.y)
    2439 
    2440     def point2ver(self,point):
    2441         return Vertex(point[0],point[1])
    2442 
    2443     def smooth_polySet(self,min_radius=0.05):
    2444         #for all pairs of connecting segments:
    2445         #    propose a new segment that replaces the 2
    2446 
    2447         #    If the difference between the new segment
    2448         #    and the old lines is small: replace the
    2449         #    old lines.
    2450 
    2451         seg2line = self.seg2line
    2452         ver2point= self.ver2point
    2453         line2seg = self.line2seg
    2454         point2ver= self.point2ver
    2455 
    2456         #create dictionaries of lines -> segments
    2457         userSegments = self.segs_to_dict(self.userSegments)
    2458         alphaSegments = self.segs_to_dict(self.alphaUserSegments)
    2459 
    2460         #lump user and alpha segments
    2461         for key in alphaSegments.keys():
    2462             userSegments[key]=alphaSegments[key]
    2463 
    2464         #point_keys = tuple -> vertex
    2465         #userVertices = vertex -> [line,line] - lines from that node
    2466         point_keys = {}
    2467         userVertices={}
    2468         for vertex in self.getUserVertices():
    2469             point = ver2point(vertex)
    2470             if not point_keys.has_key(point):
    2471                 point_keys[point]=vertex
    2472                 userVertices[vertex]=[]
    2473         for key in userSegments.keys():
    2474             line = key
    2475             point_0 = key[0]
    2476             point_1 = key[1]
    2477             userVertices[point_keys[point_0]].append(line)
    2478             userVertices[point_keys[point_1]].append(line)
    2479 
    2480         for point in point_keys.keys():
    2481             try:
    2482             #removed keys can cause keyerrors
    2483                 vertex = point_keys[point]
    2484                 lines = userVertices[vertex]
    2485    
    2486                 #if there are 2 lines on the node
    2487                 if len(lines)==2:
    2488                     line_0 = lines[0]
    2489                     line_1 = lines[1]
    2490    
    2491                     #if the tags are the the same on the 2 lines
    2492                     if userSegments[line_0].tag == userSegments[line_1].tag:
    2493                         tag = userSegments[line_0].tag
    2494    
    2495                         #point_a is one of the next nodes, point_b is the other
    2496                         if point==line_0[0]:
    2497                             point_a = line_0[1]
    2498                         if point==line_0[1]:
    2499                             point_a = line_0[0]
    2500                         if point==line_1[0]:
    2501                             point_b = line_1[1]
    2502                         if point==line_1[1]:
    2503                             point_b = line_1[0]
    2504    
    2505    
    2506                         #line_2 is proposed
    2507                         line_2 = (point_a,point_b)
    2508 
    2509                         #calculate the area of the triangle between
    2510                         #the two existing segments and the proposed
    2511                         #new segment
    2512                         ax = point_a[0]
    2513                         ay = point_a[1]
    2514                         bx = point_b[0]
    2515                         by = point_b[1]
    2516                         cx = point[0]
    2517                         cy = point[1]
    2518                         area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
    2519 
    2520                         #calculate the perimeter
    2521                         len_a =  ((cx-bx)**2+(cy-by)**2)**0.5
    2522                         len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5
    2523                         len_c =  ((bx-ax)**2+(by-ay)**2)**0.5
    2524                         perimeter = len_a+len_b+len_c
    2525 
    2526                         #calculate the radius
    2527                         r = area/(2*perimeter)
    2528 
    2529                         #if the radius is small: then replace the existing
    2530                         #segments with the new one
    2531                         if r < min_radius:
    2532                             if len_c < min_radius: append = False
    2533                             else: append = True
    2534                             #if the new seg is also time, don't add it
    2535                             if append:
    2536                                 segment = self.line2seg(line_2,tag=tag)
    2537 
    2538                             list_a=userVertices[point_keys[point_a]]
    2539                             list_b=userVertices[point_keys[point_b]]
    2540 
    2541                             if line_0 in list_a:
    2542                                 list_a.remove(line_0)
    2543                             else:
    2544                                 list_a.remove(line_1)
    2545 
    2546                             if line_0 in list_b:
    2547                                 list_b.remove(line_0)
    2548                             else:
    2549                                 list_b.remove(line_1)
    2550 
    2551                             if append:
    2552                                 list_a.append(line_2)
    2553                                 list_b.append(line_2)
    2554                             else:
    2555                                 if len(list_a)==0:
    2556                                     userVertices.pop(point_keys[point_a])
    2557                                     point_keys.pop(point_a)
    2558                                 if len(list_b)==0:
    2559                                     userVertices.pop(point_keys[point_b])
    2560                                     point_keys.pop(point_b)
    2561 
    2562                             userVertices.pop(point_keys[point])
    2563                             point_keys.pop(point)
    2564                             userSegments.pop(line_0)
    2565                             userSegments.pop(line_1)
    2566 
    2567                             if append:
    2568                                 userSegments[line_2]=segment
    2569             except:
    2570                 pass
    2571 
    2572         #self.userVerticies = userVertices.keys()
    2573         #self.userSegments = []
    2574         #for key in userSegments.keys():
    2575         #    self.userSegments.append(userSegments[key])
    2576         #self.alphaUserSegments = []
    2577 
    2578         self.userVerticies = []
    2579         self.userSegments = []
    2580         self.alphaUserSegments = []
    2581 
    2582         return userVertices,userSegments,alphaSegments
    2583 
    2584     def triangles_to_polySet(self,setName):
    2585         #self.smooth_polySet()
    2586 
    2587         seg2line = self.seg2line
    2588         ver2point= self.ver2point
    2589         line2seg = self.line2seg
    2590         point2ver= self.point2ver
    2591 
    2592         from Numeric import array,allclose
    2593         #turn the triangles into a set
    2594         Triangles = self.sets[self.setID[setName]]
    2595         Triangles_dict = {}
    2596         for triangle in Triangles:
    2597             Triangles_dict[triangle]=None
    2598  
    2599 
    2600         #create a dict of points to vertexes (tuple -> object)
    2601         #also create a set of vertexes (object -> True)
    2602         point_keys = {}
    2603         userVertices={}
    2604         for vertex in self.getUserVertices():
    2605             point = ver2point(vertex)
    2606             if not point_keys.has_key(point):
    2607                 point_keys[point]=vertex
    2608                 userVertices[vertex]=True
    2609 
    2610         #create a dict of lines to segments (tuple -> object)
    2611         userSegments = self.segs_to_dict(self.userSegments)
    2612         #append the userlines in an affine linespace
    2613         affine_lines = Affine_Linespace()
    2614         for line in userSegments.keys():
    2615             affine_lines.append(line)
    2616         alphaSegments = self.segs_to_dict(self.alphaUserSegments)
    2617         for line in alphaSegments.keys():
    2618             affine_lines.append(line)
    2619        
    2620         for triangle in Triangles:
    2621             for i in (0,1,2):
    2622                 #for every triangles neighbour:
    2623                 if not Triangles_dict.has_key(triangle.neighbors[i]):
    2624                 #if the neighbour is not in the set:
    2625                     a = triangle.vertices[i-1]
    2626                     b = triangle.vertices[i-2]
    2627                     #Get possible matches:
    2628                     point_a = ver2point(a)
    2629                     point_b = ver2point(b)
    2630                     midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
    2631                     line = (point_a,point_b)
    2632                     tag = None
    2633 
    2634 
    2635                     #this bit checks for matching lines
    2636                     possible_lines = affine_lines[line]
    2637                     possible_lines = unique(possible_lines)
    2638                     found = 0                           
    2639                     for user_line in possible_lines:
    2640                         if self.point_on_line(midpoint,user_line):
    2641                             found+=1
    2642                             assert found<2
    2643                             if userSegments.has_key(user_line):
    2644                                 parent_segment = userSegments.pop(user_line)
    2645                             if alphaSegments.has_key(user_line):
    2646                                 parent_segment = alphaSegments.pop(user_line)
    2647                             tag = parent_segment.tag
    2648                             offspring = [line]
    2649                             offspring.extend(self.subtract_line(user_line,
    2650                                                                 line))
    2651                             affine_lines.remove(user_line)
    2652                             for newline in offspring:
    2653                                 line_vertices = []
    2654                                 for point in newline:
    2655                                     if point_keys.has_key(point):
    2656                                         vert = point_keys[point]
    2657                                     else:
    2658                                         vert = Vertex(point[0],point[1])
    2659                                         userVertices[vert]=True
    2660                                         point_keys[point]=vert
    2661                                     line_vertices.append(vert)
    2662                                 segment = Segment(line_vertices[0],
    2663                                                   line_vertices[1],tag)
    2664                                 userSegments[newline]=segment
    2665                                 affine_lines.append(newline)
    2666                             #break
    2667                     assert found<2
    2668 
    2669 
    2670 
    2671                     #if no matching lines
    2672                     if not found:
    2673                         line_vertices = []
    2674                         for point in line:
    2675                             if point_keys.has_key(point):
    2676                                 vert = point_keys[point]
    2677                             else:
    2678                                 vert = Vertex(point[0],point[1])
    2679                                 userVertices[vert]=True
    2680                                 point_keys[point]=vert
    2681                             line_vertices.append(vert)
    2682                         segment = Segment(line_vertices[0],
    2683                                           line_vertices[1],tag)
    2684                         userSegments[line]=segment
    2685                         affine_lines.append(line)
    2686        
    2687         self.userVerticies = []
    2688         self.userSegments = []
    2689         self.alphaUserSegments = []
    2690 
    2691         return userVertices,userSegments,alphaSegments
    2692 
    2693     def subtract_line(self,parent,child):
    2694         #Subtracts child from parent
    2695         #Requires that the child is a
    2696         #subline of parent to work.
    2697 
    2698         from Numeric import allclose,dot,array
    2699         A= parent[0]
    2700         B= parent[1]
    2701         a = child[0]
    2702         b = child[1]
    2703 
    2704         A_array = array(parent[0])
    2705         B_array = array(parent[1])
    2706         a_array   = array(child[0])
    2707         b_array   = array(child[1])
    2708 
    2709         assert not A == B
    2710         assert not a == b
    2711 
    2712         answer = []
    2713 
    2714         #if the new line does not share a
    2715         #vertex with the old one
    2716         if not (allclose(A_array,a_array)\
    2717              or allclose(B_array,b_array)\
    2718              or allclose(A_array,b_array)\
    2719              or allclose(a_array,B_array)):
    2720             if dot(A_array-a_array,A_array-a_array) \
    2721             < dot(A_array-b_array,A_array-b_array):
    2722                 sibling1 = (A,a)
    2723                 sibling2 = (B,b)
    2724                 return [sibling1,sibling2]
    2725             else:
    2726                 sibling1 = (A,b)
    2727                 sibling2 = (B,a)
    2728                 return [sibling1,sibling2]
    2729 
    2730         elif allclose(A_array,a_array):
    2731             if allclose(B_array,b_array):
    2732                 return []
    2733             else:
    2734                 sibling = (b,B)
    2735                 return [sibling]
    2736         elif allclose(B_array,b_array):
    2737             sibling = (a,A)
    2738             return [sibling]
    2739 
    2740         elif allclose(A_array,b_array):
    2741             if allclose(B,a):
    2742                 return []
    2743             else:
    2744                 sibling = (a,B)
    2745                 return [sibling]
    2746         elif allclose(a_array,B_array):
    2747             sibling = (b,A)
    2748             return [sibling]
    2749 
    2750     def point_on_line(self,point,line):
    2751         #returns true within a tolerance of 3 degrees
    2752         x=point[0]
    2753         y=point[1]
    2754         x0=line[0][0]
    2755         x1=line[1][0]
    2756         y0=line[0][1]
    2757         y1=line[1][1]
    2758         from Numeric import array, dot, allclose
    2759         from math import sqrt
    2760         tol = 3. #DEGREES
    2761         tol = tol*3.1415/180
    2762 
    2763         a = array([x - x0, y - y0])
    2764         a_normal = array([a[1], -a[0]])     
    2765         len_a_normal = sqrt(sum(a_normal**2))
    2766 
    2767         b = array([x1 - x0, y1 - y0])         
    2768         len_b = sqrt(sum(b**2))
    2769    
    2770         if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
    2771             #Point is somewhere on the infinite extension of the line
    2772 
    2773             len_a = sqrt(sum(a**2))                                     
    2774             if dot(a, b) >= 0 and len_a <= len_b:
    2775                return True
    2776             else:   
    2777                return False
    2778         else:
    2779           return False
    2780 
    2781     def line_length(self,line):
    2782         x0=line[0][0]
    2783         x1=line[1][0]
    2784         y0=line[0][1]
    2785         y1=line[1][1]
    2786         return ((x1-x0)**2-(y1-y0)**2)**0.5     
    2787 
    2788     def threshold(self,setName,min=None,max=None,attribute_name='elevation'):
    2789         """
    2790         threshold using  d
    2791         """
    2792         triangles = self.sets[self.setID[setName]]
    2793         A = []
    2794 
    2795         if attribute_name in self.attributeTitles:
    2796             i = self.attributeTitles.index(attribute_name)
    2797         else: i = -1#no attribute
    2798         if not max == None:
    2799             for t in triangles:
    2800                 if (min<self.av_att(t,i)<max):
    2801                     A.append(t)
    2802         else:
    2803             for t in triangles:
    2804                 if (min<self.av_att(t,i)):
    2805                     A.append(t)
    2806         self.sets[self.setID[setName]] = A
    2807 
    2808     def general_threshold(self,setName,min=None,max=None\
    2809               ,attribute_name = 'elevation',function=None):
    2810         """
    2811         Thresholds the triangles
    2812         """
    2813         from visual.graph import arange,ghistogram,color as colour
    2814         triangles = self.sets[self.setID[setName]]
    2815         A = []
    2816         data=[]
    2817         #data is for the graph
    2818 
    2819         if attribute_name in self.attributeTitles:
    2820             i = self.attributeTitles.index(attribute_name)
    2821         else: i = -1
    2822         if not max == None:
    2823             for t in triangles:
    2824                 value=function(t,i)
    2825                 if (min<value<max):
    2826                     A.append(t)
    2827                 data.append(value)
    2828         else:
    2829             for t in triangles:
    2830                 value=function(t,i)
    2831                 if (min<value):
    2832                     A.append(t)
    2833                 data.append(value)
    2834         self.sets[self.setID[setName]] = A
    2835 
    2836         if self.visualise_graph:
    2837             if len(data)>0:
    2838                 max=data[0]
    2839                 min=data[0]
    2840                 for value in data:
    2841                     if value > max:
    2842                         max = value
    2843                     if value < min:
    2844                         min = value
    2845 
    2846                 inc = (max-min)/100
    2847 
    2848                 histogram = ghistogram(bins=arange(min,max,inc),\
    2849                              color = colour.red)
    2850                 histogram.plot(data=data)
    2851        
    2852     def av_att(self,triangle,i):
    2853         if i==-1: return 1
    2854         else:
    2855             #evaluates the average attribute of the vertices of a triangle.
    2856             V = triangle.getVertices()
    2857             a0 = (V[0].attributes[i])
    2858             a1 = (V[1].attributes[i])
    2859             a2 = (V[2].attributes[i])
    2860             return (a0+a1+a2)/3
    2861 
    2862     def Courant_ratio(self,triangle,index):
    2863         """
    2864         Uses the courant threshold
    2865         """
    2866         e = self.av_att(triangle,index)
    2867         A = triangle.calcArea()
    2868         P = triangle.calcP()
    2869         r = A/(2*P)
    2870         e = max(0.1,abs(e))
    2871         return r/e**0.5
    2872 
    2873     def Gradient(self,triangle,index):
    2874         V = triangle.vertices
    2875         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]
    2876         grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    2877         if ((grad_x**2)+(grad_y**2))**(0.5)<0:
    2878             print ((grad_x**2)+(grad_y**2))**(0.5)
    2879         return ((grad_x**2)+(grad_y**2))**(0.5)
    2880    
    2881 
    2882     def append_triangle(self,triangle):
    2883         self.meshTriangles.append(triangle)
    2884 
    2885     def replace_triangle(self,triangle,replacement):
    2886         i = self.meshTriangles.index(triangle)
    2887         self.meshTriangles[i]=replacement
    2888         assert replacement in self.meshTriangles
    2889 
    28902189def importUngenerateFile(ofile):
    28912190    """
     
    32102509    return u
    32112510
    3212 """Refines triangles
    3213 
    3214    Implements the #triangular bisection?# algorithm.
    3215  
    3216    
    3217 """
    3218 
    3219 def Refine(mesh, triangles):
    3220     """
    3221     Given a general_mesh, and a triangle number, split
    3222     that triangle in the mesh in half. Then to prevent
    3223     vertices and edges from meeting, keep refining
    3224     neighbouring triangles until the mesh is clean.
    3225     """
    3226     state = BisectionState(mesh)
    3227     for triangle in triangles:
    3228         if not state.refined_triangles.has_key(triangle):
    3229             triangle.rotate_longest_side()
    3230             state.start(triangle)
    3231             Refine_mesh(mesh, state)
    3232 
    3233 def Refine_mesh(mesh, state):
    3234     """
    3235     """
    3236     state.getState(mesh)
    3237     refine_triangle(mesh,state)
    3238     state.evolve()
    3239     if not state.end:
    3240         Refine_mesh(mesh,state)
    3241 
    3242 def refine_triangle(mesh,state):
    3243     split(mesh,state.current_triangle,state.new_point)
    3244     if state.case == 'one':
    3245         state.r[3]=state.current_triangle#triangle 2
    3246 
    3247         new_triangle_id = len(mesh.meshTriangles)-1
    3248         new_triangle = mesh.meshTriangles[new_triangle_id]
    3249 
    3250         split(mesh,new_triangle,state.old_point)
    3251         state.r[2]=new_triangle#triangle 1.2
    3252         state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
    3253         r = state.r
    3254         state.repairCaseOne()
    3255 
    3256     if state.case == 'two':
    3257         state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
    3258 
    3259         new_triangle = state.current_triangle
    3260 
    3261         split(mesh,new_triangle,state.old_point)
    3262 
    3263         state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
    3264         state.r[4]=new_triangle#triangle 2.2
    3265         r = state.r
    3266 
    3267         state.repairCaseTwo()
    3268 
    3269     if state.case == 'vertex':
    3270         state.r[2]=state.current_triangle#triangle 2
    3271         state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
    3272         r = state.r
    3273         state.repairCaseVertex()
    3274        
    3275     if state.case == 'start':
    3276         state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
    3277         state.r[3]=state.current_triangle#triangle 2
    3278 
    3279     if state.next_case == 'boundary':
    3280         state.repairCaseBoundary()
    3281 
    3282 
    3283 def split(mesh, triangle, new_point):
    3284         """
    3285         Given a mesh, triangle_id and a new point,
    3286         split the corrosponding triangle into two
    3287         new triangles and update the mesh.
    3288         """
    3289 
    3290         new_triangle1 = Triangle(new_point,triangle.vertices[0],
    3291                                  triangle.vertices[1],
    3292                                  attribute = triangle.attribute,
    3293                                  neighbors = None)
    3294         new_triangle2 = Triangle(new_point,triangle.vertices[2],
    3295                                  triangle.vertices[0],
    3296                                  attribute = triangle.attribute,
    3297                                  neighbors = None)
    3298 
    3299         new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
    3300         new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
    3301 
    3302         mesh.meshTriangles.append(new_triangle1)
    3303 
    3304         triangle.vertices = new_triangle2.vertices
    3305         triangle.neighbors = new_triangle2.neighbors
    3306 
    3307 
    3308 class State:
    3309 
    3310     def __init__(self):
    3311         pass
    3312 
    3313 class BisectionState(State):
    3314 
    3315 
    3316     def __init__(self,mesh):
    3317         self.len = len(mesh.meshTriangles)
    3318         self.refined_triangles = {}
    3319         self.mesh = mesh
    3320         self.current_triangle = None
    3321         self.case = 'start'
    3322         self.end = False
    3323         self.r = [None,None,None,None,None]
    3324 
    3325     def start(self, triangle):
    3326         self.current_triangle = triangle
    3327         self.case = 'start'
    3328         self.end = False
    3329         self.r = [None,None,None,None,None]
    3330 
    3331     def getState(self,mesh):
    3332         if not self.case == 'vertex':
    3333             self.new_point=self.getNewVertex(mesh, self.current_triangle)
    3334             #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
    3335             self.neighbour = self.current_triangle.neighbors[0]
    3336             if not self.neighbour is None:
    3337                 self.neighbour.rotate_longest_side()
    3338             self.next_case = self.get_next_case(mesh,self.neighbour,
    3339                                                 self.current_triangle)
    3340         if self.case == 'vertex':
    3341             self.new_point=self.old_point
    3342 
    3343 
    3344     def evolve(self):
    3345         if self.case == 'vertex':
    3346             self.end = True
    3347 
    3348         self.last_case = self.case
    3349         self.case = self.next_case
    3350 
    3351         self.old_point = self.new_point
    3352         self.current_triangle = self.neighbour
    3353 
    3354         if self.case == 'boundary':
    3355             self.end = True
    3356         self.refined_triangles[self.r[2]]=1
    3357         self.refined_triangles[self.r[3]]=1
    3358         if not self.r[4] is None:
    3359             self.refined_triangles[self.r[4]]=1
    3360         self.r[0]=self.r[2]
    3361         self.r[1]=self.r[3]
    3362 
    3363 
    3364     def getNewVertex(self,mesh,triangle):
    3365         coordinate1 = triangle.vertices[1]
    3366         coordinate2 = triangle.vertices[2]
    3367         a = ([coordinate1.x*1.,coordinate1.y*1.])
    3368         b = ([coordinate2.x*1.,coordinate2.y*1.])
    3369         attributes = []
    3370         for i in range(len(coordinate1.attributes)):
    3371             att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
    3372             attributes.append(att)
    3373         new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
    3374         newVertex = Vertex(new_coordinate[0],new_coordinate[1],
    3375                            attributes = attributes)
    3376         mesh.maxVertexIndex+=1
    3377         newVertex.index = mesh.maxVertexIndex
    3378         mesh.meshVertices.append(newVertex)
    3379         return newVertex
    3380 
    3381     def get_next_case(self, mesh,neighbour,triangle):
    3382         """
    3383         Given the locations of two neighbouring triangles,
    3384         examine the interior indices of their vertices (i.e.
    3385         0,1 or 2) to determine what how the neighbour needs
    3386         to be refined.
    3387         """
    3388         if (neighbour is None):
    3389                 next_case = 'boundary'
    3390         else:
    3391                 if triangle.vertices[1].x==neighbour.vertices[2].x:
    3392                     if triangle.vertices[1].y==neighbour.vertices[2].y:
    3393                         next_case = 'vertex'
    3394                 if triangle.vertices[1].x==neighbour.vertices[0].x:
    3395                     if triangle.vertices[1].y==neighbour.vertices[0].y:
    3396                         next_case = 'two'
    3397                 if triangle.vertices[1].x==neighbour.vertices[1].x:
    3398                     if triangle.vertices[1].y==neighbour.vertices[1].y:
    3399                         next_case = 'one'
    3400         return next_case
    3401 
    3402 
    3403 
    3404     def repairCaseVertex(self):
    3405 
    3406         r = self.r
    3407 
    3408 
    3409         self.link(r[0],r[2])
    3410         self.repair(r[0])
    3411 
    3412         self.link(r[1],r[3])
    3413         self.repair(r[1])
    3414 
    3415         self.repair(r[2])
    3416 
    3417         self.repair(r[3])
    3418 
    3419 
    3420     def repairCaseOne(self):
    3421         r = self.rkey
    3422 
    3423 
    3424         self.link(r[0],r[2])
    3425         self.repair(r[0])
    3426 
    3427         self.link(r[1],r[4])
    3428         self.repair(r[1])
    3429 
    3430         self.repair(r[4])
    3431 
    3432     def repairCaseTwo(self):
    3433         r = self.r
    3434 
    3435         self.link(r[0],r[4])
    3436         self.repair(r[0])
    3437 
    3438         self.link(r[1],r[3])
    3439         self.repair(r[1])
    3440 
    3441         self.repair(r[4])
    3442 
    3443     def repairCaseBoundary(self):
    3444         r = self.r
    3445         self.repair(r[2])
    3446         self.repair(r[3])
    3447 
    3448 
    3449 
    3450     def repair(self,triangle):
    3451         """
    3452         Given a triangle that knows its neighbours, this will
    3453         force the neighbours to comply.
    3454 
    3455         However, it needs to compare the vertices of triangles
    3456         for this implementation
    3457 
    3458         But it doesn't work for invalid neighbour structures
    3459         """
    3460         n=triangle.neighbors
    3461         for i in (0,1,2):
    3462             if not n[i] is None:
    3463                 for j in (0,1,2):#to find which side of the list is broken
    3464                     if not (n[i].vertices[j] in triangle.vertices):
    3465                     #ie if j is the side of n that needs fixing
    3466                         k = j
    3467                 n[i].neighbors[k]=triangle
    3468 
    3469 
    3470 
    3471     def link(self,triangle1,triangle2):
    3472         """
    3473         make triangle1 neighbors point to t
    3474                 #count = 0riangle2
    3475         """
    3476         count = 0
    3477         for i in (0,1,2):#to find which side of the list is broken
    3478             if not (triangle1.vertices[i] in triangle2.vertices):
    3479                 j = i
    3480                 count+=1
    3481         assert count == 1
    3482         triangle1.neighbors[j]=triangle2
    3483 
    3484 class Discretised_Tuple_Set:
    3485     """
    3486     if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
    3487     a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
    3488     a[(10000)]=[] #NOT KEYERROR
    3489 
    3490     a.append[(0.01)]
    3491     => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
    3492 
    3493     #NOT IMPLIMENTED
    3494     a.remove[(0.01)]
    3495     => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
    3496 
    3497     a.remove[(0.17)]
    3498     => {(0.0):[(0.02),(0.01)],0.2:[]}
    3499     #NOT IMPLIMENTED
    3500     at a.precision = 2:
    3501         a.round_up_rel[0.0]=
    3502         a.round_flat[0.0]=
    3503         a.round_down_rel[0.0]=
    3504 
    3505         a.up((0.1,2.04))=
    3506 
    3507     If t_rel = 0, nothing gets rounded into
    3508     two bins. If t_rel = 0.5, everything does.
    3509 
    3510     Ideally, precision can be set high, so that multiple
    3511     entries are rarely in the same bin. And t_rel should
    3512     be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
    3513     so that it is rare to put items in mutiple bins.
    3514 
    3515     Ex bins per entry = product(a,b,c...,n)
    3516     a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
    3517     b = 1 or 2 ...
    3518 
    3519     BUT!!! to avoid missing the right bin:
    3520     (-10)**(precision+1)*t_rel must be greater than the
    3521     greatest possible variation that an identical element
    3522     can display.
    3523 
    3524 
    3525     Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
    3526     but not .6 - this looks wrong, but note that *everything* will round,
    3527     so .6 wont be missed as everything close to it will check in .7 and .5.
    3528     """
    3529     def __init__(self,p_rel = 6,t_rel = 0.01):
    3530         self.__p_rel__ = p_rel
    3531         self.__t_rel__ = t_rel
    3532 
    3533         self.__p_abs__ = p_rel+1
    3534         self.__t_abs__ = t_rel
    3535 
    3536         assert t_rel <= 0.5
    3537         self.__items__ = {}
    3538         from math import frexp
    3539         self.frexp = frexp
    3540         roundings = [self.round_up_rel,\
    3541         self.round_down_rel,self.round_flat_rel,\
    3542         self.round_down_abs,self.round_up_abs,\
    3543         self.round_flat_abs]#
    3544 
    3545         self.roundings = roundings
    3546 
    3547     def __repr__(self):
    3548         return '%s'%self.__items__
    3549 
    3550     def rounded_keys(self,key):
    3551         key = tuple(key)
    3552         keys = [key]
    3553         keys = self.__rounded_keys__(key)
    3554         return (keys)
    3555 
    3556     def __rounded_keys__(self,key):
    3557         keys = []
    3558         rounded_key=list(key)
    3559         rounded_values=list(key)
    3560 
    3561         roundings = list(self.roundings)
    3562 
    3563         #initialise rounded_values
    3564         round = roundings.pop(0)
    3565         for i in range(len(rounded_values)):
    3566             rounded_key[i]=round(key[i])
    3567             rounded_values[i]={}
    3568             rounded_values[i][rounded_key[i]]=None
    3569         keys.append(tuple(rounded_key))
    3570        
    3571         for round in roundings:
    3572             for i in range(len(rounded_key)):
    3573                 rounded_value=round(key[i])
    3574                 if not rounded_values[i].has_key(rounded_value):
    3575                     #ie unless round_up_rel = round_down_rel
    3576                     #so the keys stay unique
    3577                     for j in range(len(keys)):
    3578                         rounded_key = list(keys[j])
    3579                         rounded_key[i]=rounded_value
    3580                         keys.append(tuple(rounded_key))
    3581         return keys
    3582 
    3583     def append(self,item):
    3584         keys = self.rounded_keys(item)
    3585         for key in keys:
    3586             if self.__items__.has_key(key):
    3587                 self.__items__[key].append(item)
    3588             else:
    3589                 self.__items__[key]=[item]
    3590 
    3591     def __getitem__(self,key):
    3592         answer = []
    3593         keys = self.rounded_keys(key)
    3594         for key in keys:
    3595             if self.__items__.has_key(key):
    3596                 answer.extend(self.__items__[key])
    3597         #if len(answer)==0:
    3598         #    raise KeyError#FIXME or return KeyError
    3599         #                  #FIXME or just return []?
    3600         else:
    3601             return answer #FIXME or unique(answer)?
    3602 
    3603     def __delete__(self,item):
    3604         keys = self.rounded_keys(item)
    3605         answer = False
    3606         #if any of the possible keys contains
    3607         #a list, return true
    3608         for key in keys:       
    3609             if self.__items__.has_key(key):
    3610                 if item in self.__items__[key]:
    3611                     self.__items__[key].remove(item)
    3612 
    3613     def remove(self,item):
    3614         self.__delete__(item)
    3615 
    3616     def __contains__(self,item):
    3617 
    3618         keys = self.rounded_keys(item)
    3619         answer = False
    3620         #if any of the possible keys contains
    3621         #a list, return true
    3622         for key in keys:       
    3623             if self.__items__.has_key(key):
    3624                 if item in self.__items__[key]:
    3625                     answer = True
    3626         return answer
    3627 
    3628 
    3629     def has_item(self,item):
    3630         return self.__contains__(item)
    3631 
    3632     def round_up_rel2(self,value):
    3633          t_rel=self.__t_rel__
    3634          #Rounding up the value
    3635          m,e = self.frexp(value)
    3636          m = m/2
    3637          e = e + 1
    3638          #m is the mantissa, e the exponent
    3639          # 0.5 < |m| < 1.0
    3640          m = m+t_rel*(10**-(self.__p_rel__))
    3641          #bump m up
    3642          m = round(m,self.__p_rel__)
    3643          return m*(2.**e)
    3644 
    3645     def round_down_rel2(self,value):
    3646          t_rel=self.__t_rel__
    3647          #Rounding down the value
    3648          m,e = self.frexp(value)
    3649          m = m/2
    3650          e = e + 1
    3651          #m is the mantissa, e the exponent
    3652          # 0.5 < m < 1.0
    3653          m = m-t_rel*(10**-(self.__p_rel__))
    3654          #bump the |m| down, by 5% or whatever
    3655          #self.p_rel dictates
    3656          m = round(m,self.__p_rel__)
    3657          return m*(2.**e)
    3658 
    3659     def round_flat_rel2(self,value):
    3660     #redundant
    3661          m,e = self.frexp(value)
    3662          m = m/2
    3663          e = e + 1
    3664          m = round(m,self.__p_rel__)
    3665          return m*(2.**e)
    3666 
    3667     def round_up_rel(self,value):
    3668          t_rel=self.__t_rel__
    3669          #Rounding up the value
    3670          m,e = self.frexp(value)
    3671          #m is the mantissa, e the exponent
    3672          # 0.5 < |m| < 1.0
    3673          m = m+t_rel*(10**-(self.__p_rel__))
    3674          #bump m up
    3675          m = round(m,self.__p_rel__)
    3676          return m*(2.**e)
    3677 
    3678     def round_down_rel(self,value):
    3679          t_rel=self.__t_rel__
    3680          #Rounding down the value
    3681          m,e = self.frexp(value)
    3682          #m is the mantissa, e the exponent
    3683          # 0.5 < m < 1.0
    3684          m = m-t_rel*(10**-(self.__p_rel__))
    3685          #bump the |m| down, by 5% or whatever
    3686          #self.p_rel dictates
    3687          m = round(m,self.__p_rel__)
    3688          return m*(2.**e)
    3689 
    3690     def round_flat_rel(self,value):
    3691     #redundant
    3692          m,e = self.frexp(value)
    3693          m = round(m,self.__p_rel__)
    3694          return m*(2.**e)
    3695 
    3696     def round_up_abs(self,value):
    3697          t_abs=self.__t_abs__
    3698          #Rounding up the value
    3699          m = value+t_abs*(10**-(self.__p_abs__))
    3700          #bump m up
    3701          m = round(m,self.__p_abs__)
    3702          return m
    3703 
    3704     def round_down_abs(self,value):
    3705          t_abs=self.__t_abs__
    3706          #Rounding down the value
    3707          m = value-t_abs*(10**-(self.__p_abs__))
    3708          #bump the |m| down, by 5% or whatever
    3709          #self.p_rel dictates
    3710          m = round(m,self.__p_abs__)
    3711          return m
    3712 
    3713     def round_flat_abs(self,value):
    3714     #redundant?
    3715          m = round(value,self.__p_abs__)
    3716          return m
    3717 
    3718     def keys(self):
    3719         return self.__items__.keys()
    3720 
    3721 
    3722 class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
    3723     """
    3724     This is a discretised tuple set, but
    3725     based on a mapping. The mapping MUST
    3726     return a sequence.
    3727 
    3728     example:
    3729     def weight(animal):
    3730         return [animal.weight]
    3731    
    3732     a = Mapped_Discretised_Tuple_Set(weight)
    3733     a.append[cow]
    3734     a.append[fox]
    3735     a.append[horse]
    3736 
    3737     a[horse]    -> [cow,horse]
    3738     a[dog]      -> [fox]
    3739     a[elephant] -> []
    3740     """
    3741     def __init__(self,mapping,p_rel = 6, t_rel=0.01):
    3742         Discretised_Tuple_Set.__init__\
    3743          (self,p_rel,t_rel = t_rel)
    3744         self.mapping = mapping
    3745 
    3746     def rounded_keys(self,key):
    3747         mapped_key = tuple(self.mapping(key))
    3748         keys = self.__rounded_keys__(mapped_key)
    3749         return keys
    3750 
    3751 class Affine_Linespace(Mapped_Discretised_Tuple_Set):
    3752     """
    3753     The affine linespace creates a way to record and compare lines.
    3754     Precision is a bit of a hack, but it creates a way to avoid
    3755     misses caused by round offs (between two lines of different
    3756     lenghts, the short one gets rounded off more).
    3757     I am starting to think that a quadratic search would be faster.
    3758     Nearly.
    3759     """
    3760     def __init__(self,p_rel=4,t_rel=0.2):
    3761         Mapped_Discretised_Tuple_Set.__init__\
    3762             (self,self.affine_line,\
    3763             p_rel=p_rel,t_rel=t_rel)
    3764 
    3765         roundings = \
    3766         [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
    3767         self.roundings = roundings
    3768         #roundings = \
    3769         #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
    3770         #self.roundings = roundings
    3771 
    3772     def affine_line(self,line):
    3773         point_1 = line[0]
    3774         point_2 = line[1]
    3775         #returns the equation of a line
    3776         #between two points, in the from
    3777         #(a,b,-c), as in ax+by-c=0
    3778         #or line *dot* (x,y,1) = (0,0,0)
    3779 
    3780         #Note that it normalises the line
    3781         #(a,b,-c) so Line*Line = 1.
    3782         #This does not change the mathematical
    3783         #properties, but it makes comparism
    3784         #easier.
    3785 
    3786         #There are probably better algorithms.
    3787         x1 = point_1[0]
    3788         x2 = point_2[0]
    3789         y1 = point_1[1]
    3790         y2 = point_2[1]
    3791         dif_x = x1-x2
    3792         dif_y = y1-y2
    3793 
    3794         if dif_x == dif_y == 0:
    3795             msg = 'points are the same'
    3796             raise msg
    3797         elif abs(dif_x)>=abs(dif_y):
    3798             alpha = (-dif_y)/dif_x
    3799             #a = alpha * b
    3800             b = -1.
    3801             c = (x1*alpha+x2*alpha+y1+y2)/2.
    3802             a = alpha*b
    3803         else:
    3804             beta = dif_x/(-dif_y)
    3805             #b = beta * a
    3806             a = 1.
    3807             c = (x1+x2+y1*beta+y2*beta)/2.
    3808             b = beta*a
    3809         mag = abs(a)+abs(b)
    3810         #This does not change the mathematical
    3811         #properties, but it makes comparism possible.
    3812 
    3813         #note that the gradient is b/a, or (a/b)**-1.
    3814         #so
    3815 
    3816         #if a == 0:
    3817         #    sign_a = 1.
    3818         #else:
    3819         #    sign_a = a/((a**2)**0.5)
    3820         #if b == 0:
    3821         #    sign_b = 1.
    3822         #else:
    3823         #    sign_b = b/((b**2)**0.5)
    3824         #if c == 0:
    3825         #    sign_c = 1.
    3826         #else:
    3827         #    sign_c = c/((c**2)**0.5)
    3828         #a = a/mag*sign_a
    3829         #b = b/mag*sign_b
    3830         #c = c/mag*sign_c
    3831         a = a/mag
    3832         b = b/mag
    3833         c = c/mag
    3834         return a,b,c
    3835 
    38362511
    38372512# FIXME (DSG-DSG)
  • anuga_core/source/anuga/pmesh/test_mesh.py

    r4905 r4906  
    3636                        'Point DistanceToPoint is wrong!')
    3737       
    38     def testTriangle(self):
    39         a = Vertex (0.0, 0.0)
    40         b = Vertex (0.0, 2.0)
    41         c = Vertex (2.0,0.0)
    42         d = Vertex (0.0, 4.0)
    43         e = Vertex (2.0, 2.0)
    44         f = Vertex (4.0,0.0)
    45        
    46         t1 = Triangle(b,a,c)       
    47         t2 = Triangle(b,c,e)     
    48         t3 = Triangle(e,c,f)     
    49         t4 = Triangle(d,b,e)
    50         t2.setNeighbors(t3,t4,t1)
    51        
    52         self.failUnless( t2.neighbors[2].vertices[0] == b, 'Triangle initialisation is wrong!')
    53    
    5438       
    5539    def testSegment(self):
     
    274258                        'Vertex deleted, instead of segment.')
    275259
    276     def testTriangleArea(self):
    277         a = Vertex (10.0, 10.0)
    278         b = Vertex (10.0, 20.0)
    279         c = Vertex (20.0,10.0)
    280        
    281         d = Vertex (-20.0, 0.0)
    282         e = Vertex (-20.0, -20.0)
    283         f = Vertex (20.0,-20.0)
    284        
    285         t1 = Triangle(b,a,c)     
    286         t2 = Triangle(e,d,f)
    287        
    288 #         print "t1", t1
    289 #         print "t1 area ", t1.calcArea()
    290 #         print "t2", t2
    291 #         print "t2 area ", t2.calcArea()
    292         self.failUnless( t1.calcArea() == 50 and t2.calcArea() == 400, 'Triangle area is wrong!')
    293260    def testisUserSegmentNew (self):
    294261        mesh = Mesh()
     
    17531720
    17541721
    1755 #___________beginning of Peters tests
    1756 
    1757     def test_set_stuff(self):
    1758         """
    1759         Documentation
    1760         """
    1761         #making a test mesh
    1762         p0=[0.,2.]
    1763         p1=[1.,2.]
    1764         p2=[0.,1.]
    1765         p3=[1.,1.]
    1766         p4=[0.,0.]
    1767         p5=[2.,0.]
    1768         p6=[-1.,1.]
    1769         point_list = [p0,p1,p2,p3,p4,p5,p6]
    1770 
    1771         a0=[0]
    1772         a1=[0]
    1773         a2=[100]
    1774         a3=[0]
    1775         a4=[0]
    1776         a5=[0]
    1777         a6=[0]
    1778         attribute=[a0,a1,a2,a3,a4,a5,a6]
    1779        
    1780         t0=[0,3,1]
    1781         t1=[0,2,3]
    1782         t2=[2,4,3]
    1783         t3=[4,5,3]
    1784         t4=[1,3,5]
    1785         t5=[2,6,4]
    1786 
    1787         n0=[4,-1,2]
    1788         n1=[2,0,-1]
    1789         n2=[3,1,5]
    1790         n3=[4,2,-1]
    1791         n4=[3,-1,0]
    1792         n5=[-1,2,-1]
    1793 
    1794         tri_list = [t0,t1,t2,t3,t4,t5]
    1795         n_list = [n0,n1,n2,n3,n4,n5]
    1796         for i in range(6):
    1797             for j in (0,1,2):
    1798                 a=attribute[tri_list[i][j]]
    1799                 tri_list[i][j]=point_list[tri_list[i][j]]
    1800                 tri_list[i][j]=Vertex(tri_list[i][j][0]\
    1801                                       ,tri_list[i][j][1],a)
    1802             neighbours=n_list[i]
    1803             tri_list[i]=Triangle(tri_list[i][0],\
    1804                                  tri_list[i][1],tri_list[i][2]\
    1805                                 ,neighbors=neighbours)
    1806 
    1807         #testing selectAll
    1808         mesh = Mesh()
    1809         mesh.attributeTitles=['attrib']
    1810 
    1811         mesh.meshTriangles=tri_list
    1812 
    1813         mesh.selectAllTriangles()
    1814         A=mesh.sets[mesh.setID['All']]
    1815         assert list_comp(tri_list,A)
    1816 
    1817        #testing threshold
    1818         mesh = Mesh()
    1819         mesh.attributeTitles=['attrib']
    1820 
    1821         mesh.meshTriangles=tri_list
    1822         mesh.selectAllTriangles()
    1823         mesh.threshold('All',min=30,max=35,attribute_name = 'attrib')
    1824         A = [tri_list[1],tri_list[2],tri_list[5]]
    1825         B = mesh.sets[mesh.setID['All']]
    1826         assert list_comp(A,B)
    1827        
    1828 
    1829         A = [tri_list[3],tri_list[2],tri_list[5]]
    1830         assert not list_comp(A,B)
    1831 
    1832         #testing
    1833 
    1834     def test_Discretised_Tuple_Set_rounding(self):
    1835         #This is the hardest bit of DST
    1836 
    1837         tol = 0.1
    1838         a=Discretised_Tuple_Set(p_rel=1,t_rel= tol)
    1839         m = 0.541
    1840         m_up = 0.6
    1841         m_down = 0.5
    1842         assert m_up == a.round_up_rel(m)
    1843         assert m_down == a.round_down_rel(m)
    1844 
    1845         tol = 0.1
    1846         a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
    1847         m = 0.539
    1848         m_up = 0.5
    1849         m_down = 0.5
    1850         assert m_up == a.round_up_rel(m)
    1851         assert m_down == a.round_down_rel(m)
    1852 
    1853         tol = 0.5
    1854         a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
    1855 
    1856 
    1857         m = 0.6
    1858         m_up = 0.7
    1859         m_down = 0.5
    1860         assert m_up == a.round_up_rel(m)
    1861         assert m_down == a.round_down_rel(m)
    1862 
    1863         m = 0.599
    1864         m_up = 0.6
    1865         m_down = 0.5
    1866         assert m_up == a.round_up_rel(m)
    1867         assert m_down == a.round_down_rel(m)
    1868 
    1869     def test_Discretised_Tuple_Set_get(self):
    1870        
    1871         tol = 0.25
    1872         a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
    1873         b = (1.1,1.1)
    1874         a.append(b)
    1875         list = [(1.2,1.),(1.,1.),(1.,1.2),(1.2,1.2)]
    1876         for key in list:
    1877             assert a[key][0]==b
    1878             assert len(a[key])==1
    1879        
    1880         c = (2.1,1.)
    1881         a.append(c)
    1882         assert a[(2.,1.)][0]==c
    1883         assert a[(2.2,1.)][0]==c
    1884 
    1885     def test_mapped_Discretised_Tuple_Set(self):
    1886 
    1887         def map(sequence):
    1888             return [len(sequence)]
    1889 
    1890         tol = 0.5
    1891         a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
    1892         b = range(20)
    1893         a.append(b)
    1894         assert b in a[range(17)]
    1895         assert b in a[range(22)]
    1896 
    1897         tol = 0.01
    1898         a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
    1899         b = range(20)
    1900         a.append(b)
    1901         assert b in a[range(20)]
    1902         assert b in a[range(19)]
    1903         assert not range(17) in a
    1904 
    1905 #___________end of Peters tests
    1906 
    19071722    def test_add_region_from_polygon(self):
    19081723        m=Mesh()
Note: See TracChangeset for help on using the changeset viewer.