Changeset 4906
- Timestamp:
- Jan 4, 2008, 2:41:09 PM (16 years ago)
- Location:
- anuga_core/source
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/pmesh/mesh.py
r4905 r4906 53 53 kinds = _kinds() 54 54 55 SET_COLOUR='red'56 57 55 #FIXME: this is not tested. 58 56 from anuga.utilities.numerical_tools import gradient … … 299 297 self.getTag()) 300 298 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 string318 319 if neighbors is None:320 self.neighbors=[]321 else:322 self.neighbors=neighbors323 324 def replace(self,new_triangle):325 self = new_triangle326 327 def longestSideID(self):328 ax = self.vertices[0].x329 ay = self.vertices[0].y330 331 bx = self.vertices[1].x332 by = self.vertices[1].y333 334 cx = self.vertices[2].x335 cy = self.vertices[2].y336 337 lenA = ((cx-bx)**2+(cy-by)**2)**0.5338 lenB = ((ax-cx)**2+(ay-cy)**2)**0.5339 lenC = ((bx-ax)**2+(by-ay)**2)**0.5340 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 triangle347 offset must be 0,1 or 2348 """349 350 if offset == 0:351 pass352 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.vertices369 370 def get_vertices(self):371 """372 Return a list of the vertices. The x and y values will be relative373 Easting and Northings for the zone of the current geo_ref.374 """375 return self.vertices376 377 def calcArea(self):378 ax = self.vertices[0].x379 ay = self.vertices[0].y380 381 bx = self.vertices[1].x382 by = self.vertices[1].y383 384 cx = self.vertices[2].x385 cy = self.vertices[2].y386 387 return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2388 389 def calcP(self):390 #calculate the perimeter391 ax = self.vertices[0].x392 ay = self.vertices[0].y393 394 bx = self.vertices[1].x395 by = self.vertices[1].y396 397 cx = self.vertices[2].x398 cy = self.vertices[2].y399 400 a = ((cx-bx)**2+(cy-by)**2)**0.5401 b = ((ax-cx)**2+(ay-cy)**2)**0.5402 c = ((bx-ax)**2+(by-ay)**2)**0.5403 404 return a+b+c405 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 neighbor410 """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 neighbor417 """418 self.attribute = attribute #this is a string419 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 objectID428 """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 438 299 class Segment(MeshObject): 439 300 """ … … 631 492 self.tri_mesh=None 632 493 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 646 495 self.visualise_graph = True 647 496 … … 2338 2187 self.regions.append(Object) 2339 2188 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']]=A2355 return 'All'2356 # and objectIDs2357 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 canvas2368 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 canvas2376 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 #Depreciated2384 #weed out existing duplicates2385 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]=vertex2394 #inlined would looks very ugly2395 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]=segment2413 Vertices=point_keys.values()2414 Segments=line_keys.values()2415 return Vertices,Segments2416 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]=segment2426 return dict2427 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 22446 2447 # If the difference between the new segment2448 # and the old lines is small: replace the2449 # old lines.2450 2451 seg2line = self.seg2line2452 ver2point= self.ver2point2453 line2seg = self.line2seg2454 point2ver= self.point2ver2455 2456 #create dictionaries of lines -> segments2457 userSegments = self.segs_to_dict(self.userSegments)2458 alphaSegments = self.segs_to_dict(self.alphaUserSegments)2459 2460 #lump user and alpha segments2461 for key in alphaSegments.keys():2462 userSegments[key]=alphaSegments[key]2463 2464 #point_keys = tuple -> vertex2465 #userVertices = vertex -> [line,line] - lines from that node2466 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]=vertex2472 userVertices[vertex]=[]2473 for key in userSegments.keys():2474 line = key2475 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 keyerrors2483 vertex = point_keys[point]2484 lines = userVertices[vertex]2485 2486 #if there are 2 lines on the node2487 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 lines2492 if userSegments[line_0].tag == userSegments[line_1].tag:2493 tag = userSegments[line_0].tag2494 2495 #point_a is one of the next nodes, point_b is the other2496 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 proposed2507 line_2 = (point_a,point_b)2508 2509 #calculate the area of the triangle between2510 #the two existing segments and the proposed2511 #new segment2512 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))/22519 2520 #calculate the perimeter2521 len_a = ((cx-bx)**2+(cy-by)**2)**0.52522 len_b = ((ax-cx)**2+(ay-cy)**2)**0.52523 len_c = ((bx-ax)**2+(by-ay)**2)**0.52524 perimeter = len_a+len_b+len_c2525 2526 #calculate the radius2527 r = area/(2*perimeter)2528 2529 #if the radius is small: then replace the existing2530 #segments with the new one2531 if r < min_radius:2532 if len_c < min_radius: append = False2533 else: append = True2534 #if the new seg is also time, don't add it2535 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]=segment2569 except:2570 pass2571 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,alphaSegments2583 2584 def triangles_to_polySet(self,setName):2585 #self.smooth_polySet()2586 2587 seg2line = self.seg2line2588 ver2point= self.ver2point2589 line2seg = self.line2seg2590 point2ver= self.point2ver2591 2592 from Numeric import array,allclose2593 #turn the triangles into a set2594 Triangles = self.sets[self.setID[setName]]2595 Triangles_dict = {}2596 for triangle in Triangles:2597 Triangles_dict[triangle]=None2598 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]=vertex2608 userVertices[vertex]=True2609 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 linespace2613 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 = None2633 2634 2635 #this bit checks for matching lines2636 possible_lines = affine_lines[line]2637 possible_lines = unique(possible_lines)2638 found = 02639 for user_line in possible_lines:2640 if self.point_on_line(midpoint,user_line):2641 found+=12642 assert found<22643 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.tag2648 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]=True2660 point_keys[point]=vert2661 line_vertices.append(vert)2662 segment = Segment(line_vertices[0],2663 line_vertices[1],tag)2664 userSegments[newline]=segment2665 affine_lines.append(newline)2666 #break2667 assert found<22668 2669 2670 2671 #if no matching lines2672 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]=True2680 point_keys[point]=vert2681 line_vertices.append(vert)2682 segment = Segment(line_vertices[0],2683 line_vertices[1],tag)2684 userSegments[line]=segment2685 affine_lines.append(line)2686 2687 self.userVerticies = []2688 self.userSegments = []2689 self.alphaUserSegments = []2690 2691 return userVertices,userSegments,alphaSegments2692 2693 def subtract_line(self,parent,child):2694 #Subtracts child from parent2695 #Requires that the child is a2696 #subline of parent to work.2697 2698 from Numeric import allclose,dot,array2699 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 == B2710 assert not a == b2711 2712 answer = []2713 2714 #if the new line does not share a2715 #vertex with the old one2716 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 degrees2752 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, allclose2759 from math import sqrt2760 tol = 3. #DEGREES2761 tol = tol*3.1415/1802762 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 line2772 2773 len_a = sqrt(sum(a**2))2774 if dot(a, b) >= 0 and len_a <= len_b:2775 return True2776 else:2777 return False2778 else:2779 return False2780 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.52787 2788 def threshold(self,setName,min=None,max=None,attribute_name='elevation'):2789 """2790 threshold using d2791 """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 attribute2798 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]] = A2807 2808 def general_threshold(self,setName,min=None,max=None\2809 ,attribute_name = 'elevation',function=None):2810 """2811 Thresholds the triangles2812 """2813 from visual.graph import arange,ghistogram,color as colour2814 triangles = self.sets[self.setID[setName]]2815 A = []2816 data=[]2817 #data is for the graph2818 2819 if attribute_name in self.attributeTitles:2820 i = self.attributeTitles.index(attribute_name)2821 else: i = -12822 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]] = A2835 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 = value2843 if value < min:2844 min = value2845 2846 inc = (max-min)/1002847 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 12854 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)/32861 2862 def Courant_ratio(self,triangle,index):2863 """2864 Uses the courant threshold2865 """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.52872 2873 def Gradient(self,triangle,index):2874 V = triangle.vertices2875 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]=replacement2888 assert replacement in self.meshTriangles2889 2890 2189 def importUngenerateFile(ofile): 2891 2190 """ … … 3210 2509 return u 3211 2510 3212 """Refines triangles3213 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, split3222 that triangle in the mesh in half. Then to prevent3223 vertices and edges from meeting, keep refining3224 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 23246 3247 new_triangle_id = len(mesh.meshTriangles)-13248 new_triangle = mesh.meshTriangles[new_triangle_id]3249 3250 split(mesh,new_triangle,state.old_point)3251 state.r[2]=new_triangle#triangle 1.23252 state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.13253 r = state.r3254 state.repairCaseOne()3255 3256 if state.case == 'two':3257 state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 13258 3259 new_triangle = state.current_triangle3260 3261 split(mesh,new_triangle,state.old_point)3262 3263 state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.13264 state.r[4]=new_triangle#triangle 2.23265 r = state.r3266 3267 state.repairCaseTwo()3268 3269 if state.case == 'vertex':3270 state.r[2]=state.current_triangle#triangle 23271 state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 13272 r = state.r3273 state.repairCaseVertex()3274 3275 if state.case == 'start':3276 state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 13277 state.r[3]=state.current_triangle#triangle 23278 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 two3287 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.vertices3305 triangle.neighbors = new_triangle2.neighbors3306 3307 3308 class State:3309 3310 def __init__(self):3311 pass3312 3313 class BisectionState(State):3314 3315 3316 def __init__(self,mesh):3317 self.len = len(mesh.meshTriangles)3318 self.refined_triangles = {}3319 self.mesh = mesh3320 self.current_triangle = None3321 self.case = 'start'3322 self.end = False3323 self.r = [None,None,None,None,None]3324 3325 def start(self, triangle):3326 self.current_triangle = triangle3327 self.case = 'start'3328 self.end = False3329 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_point3342 3343 3344 def evolve(self):3345 if self.case == 'vertex':3346 self.end = True3347 3348 self.last_case = self.case3349 self.case = self.next_case3350 3351 self.old_point = self.new_point3352 self.current_triangle = self.neighbour3353 3354 if self.case == 'boundary':3355 self.end = True3356 self.refined_triangles[self.r[2]]=13357 self.refined_triangles[self.r[3]]=13358 if not self.r[4] is None:3359 self.refined_triangles[self.r[4]]=13360 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])/23372 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+=13377 newVertex.index = mesh.maxVertexIndex3378 mesh.meshVertices.append(newVertex)3379 return newVertex3380 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 needs3386 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_case3401 3402 3403 3404 def repairCaseVertex(self):3405 3406 r = self.r3407 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.rkey3422 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.r3434 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.r3445 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 will3453 force the neighbours to comply.3454 3455 However, it needs to compare the vertices of triangles3456 for this implementation3457 3458 But it doesn't work for invalid neighbour structures3459 """3460 n=triangle.neighbors3461 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 broken3464 if not (n[i].vertices[j] in triangle.vertices):3465 #ie if j is the side of n that needs fixing3466 k = j3467 n[i].neighbors[k]=triangle3468 3469 3470 3471 def link(self,triangle1,triangle2):3472 """3473 make triangle1 neighbors point to t3474 #count = 0riangle23475 """3476 count = 03477 for i in (0,1,2):#to find which side of the list is broken3478 if not (triangle1.vertices[i] in triangle2.vertices):3479 j = i3480 count+=13481 assert count == 13482 triangle1.neighbors[j]=triangle23483 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 KEYERROR3489 3490 a.append[(0.01)]3491 => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}3492 3493 #NOT IMPLIMENTED3494 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 IMPLIMENTED3500 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 into3508 two bins. If t_rel = 0.5, everything does.3509 3510 Ideally, precision can be set high, so that multiple3511 entries are rarely in the same bin. And t_rel should3512 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_rel3517 b = 1 or 2 ...3518 3519 BUT!!! to avoid missing the right bin:3520 (-10)**(precision+1)*t_rel must be greater than the3521 greatest possible variation that an identical element3522 can display.3523 3524 3525 Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .53526 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_rel3531 self.__t_rel__ = t_rel3532 3533 self.__p_abs__ = p_rel+13534 self.__t_abs__ = t_rel3535 3536 assert t_rel <= 0.53537 self.__items__ = {}3538 from math import frexp3539 self.frexp = frexp3540 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 = roundings3546 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_values3564 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]]=None3569 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_rel3576 #so the keys stay unique3577 for j in range(len(keys)):3578 rounded_key = list(keys[j])3579 rounded_key[i]=rounded_value3580 keys.append(tuple(rounded_key))3581 return keys3582 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 KeyError3599 # #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 = False3606 #if any of the possible keys contains3607 #a list, return true3608 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 = False3620 #if any of the possible keys contains3621 #a list, return true3622 for key in keys:3623 if self.__items__.has_key(key):3624 if item in self.__items__[key]:3625 answer = True3626 return answer3627 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 value3635 m,e = self.frexp(value)3636 m = m/23637 e = e + 13638 #m is the mantissa, e the exponent3639 # 0.5 < |m| < 1.03640 m = m+t_rel*(10**-(self.__p_rel__))3641 #bump m up3642 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 value3648 m,e = self.frexp(value)3649 m = m/23650 e = e + 13651 #m is the mantissa, e the exponent3652 # 0.5 < m < 1.03653 m = m-t_rel*(10**-(self.__p_rel__))3654 #bump the |m| down, by 5% or whatever3655 #self.p_rel dictates3656 m = round(m,self.__p_rel__)3657 return m*(2.**e)3658 3659 def round_flat_rel2(self,value):3660 #redundant3661 m,e = self.frexp(value)3662 m = m/23663 e = e + 13664 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 value3670 m,e = self.frexp(value)3671 #m is the mantissa, e the exponent3672 # 0.5 < |m| < 1.03673 m = m+t_rel*(10**-(self.__p_rel__))3674 #bump m up3675 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 value3681 m,e = self.frexp(value)3682 #m is the mantissa, e the exponent3683 # 0.5 < m < 1.03684 m = m-t_rel*(10**-(self.__p_rel__))3685 #bump the |m| down, by 5% or whatever3686 #self.p_rel dictates3687 m = round(m,self.__p_rel__)3688 return m*(2.**e)3689 3690 def round_flat_rel(self,value):3691 #redundant3692 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 value3699 m = value+t_abs*(10**-(self.__p_abs__))3700 #bump m up3701 m = round(m,self.__p_abs__)3702 return m3703 3704 def round_down_abs(self,value):3705 t_abs=self.__t_abs__3706 #Rounding down the value3707 m = value-t_abs*(10**-(self.__p_abs__))3708 #bump the |m| down, by 5% or whatever3709 #self.p_rel dictates3710 m = round(m,self.__p_abs__)3711 return m3712 3713 def round_flat_abs(self,value):3714 #redundant?3715 m = round(value,self.__p_abs__)3716 return m3717 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, but3725 based on a mapping. The mapping MUST3726 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 = mapping3745 3746 def rounded_keys(self,key):3747 mapped_key = tuple(self.mapping(key))3748 keys = self.__rounded_keys__(mapped_key)3749 return keys3750 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 avoid3755 misses caused by round offs (between two lines of different3756 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 = roundings3768 #roundings = \3769 #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]3770 #self.roundings = roundings3771 3772 def affine_line(self,line):3773 point_1 = line[0]3774 point_2 = line[1]3775 #returns the equation of a line3776 #between two points, in the from3777 #(a,b,-c), as in ax+by-c=03778 #or line *dot* (x,y,1) = (0,0,0)3779 3780 #Note that it normalises the line3781 #(a,b,-c) so Line*Line = 1.3782 #This does not change the mathematical3783 #properties, but it makes comparism3784 #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-x23792 dif_y = y1-y23793 3794 if dif_x == dif_y == 0:3795 msg = 'points are the same'3796 raise msg3797 elif abs(dif_x)>=abs(dif_y):3798 alpha = (-dif_y)/dif_x3799 #a = alpha * b3800 b = -1.3801 c = (x1*alpha+x2*alpha+y1+y2)/2.3802 a = alpha*b3803 else:3804 beta = dif_x/(-dif_y)3805 #b = beta * a3806 a = 1.3807 c = (x1+x2+y1*beta+y2*beta)/2.3808 b = beta*a3809 mag = abs(a)+abs(b)3810 #This does not change the mathematical3811 #properties, but it makes comparism possible.3812 3813 #note that the gradient is b/a, or (a/b)**-1.3814 #so3815 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_a3829 #b = b/mag*sign_b3830 #c = c/mag*sign_c3831 a = a/mag3832 b = b/mag3833 c = c/mag3834 return a,b,c3835 3836 2511 3837 2512 # FIXME (DSG-DSG) -
anuga_core/source/anuga/pmesh/test_mesh.py
r4905 r4906 36 36 'Point DistanceToPoint is wrong!') 37 37 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 54 38 55 39 def testSegment(self): … … 274 258 'Vertex deleted, instead of segment.') 275 259 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", t1289 # print "t1 area ", t1.calcArea()290 # print "t2", t2291 # print "t2 area ", t2.calcArea()292 self.failUnless( t1.calcArea() == 50 and t2.calcArea() == 400, 'Triangle area is wrong!')293 260 def testisUserSegmentNew (self): 294 261 mesh = Mesh() … … 1753 1720 1754 1721 1755 #___________beginning of Peters tests1756 1757 def test_set_stuff(self):1758 """1759 Documentation1760 """1761 #making a test mesh1762 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 selectAll1808 mesh = Mesh()1809 mesh.attributeTitles=['attrib']1810 1811 mesh.meshTriangles=tri_list1812 1813 mesh.selectAllTriangles()1814 A=mesh.sets[mesh.setID['All']]1815 assert list_comp(tri_list,A)1816 1817 #testing threshold1818 mesh = Mesh()1819 mesh.attributeTitles=['attrib']1820 1821 mesh.meshTriangles=tri_list1822 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 #testing1833 1834 def test_Discretised_Tuple_Set_rounding(self):1835 #This is the hardest bit of DST1836 1837 tol = 0.11838 a=Discretised_Tuple_Set(p_rel=1,t_rel= tol)1839 m = 0.5411840 m_up = 0.61841 m_down = 0.51842 assert m_up == a.round_up_rel(m)1843 assert m_down == a.round_down_rel(m)1844 1845 tol = 0.11846 a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)1847 m = 0.5391848 m_up = 0.51849 m_down = 0.51850 assert m_up == a.round_up_rel(m)1851 assert m_down == a.round_down_rel(m)1852 1853 tol = 0.51854 a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)1855 1856 1857 m = 0.61858 m_up = 0.71859 m_down = 0.51860 assert m_up == a.round_up_rel(m)1861 assert m_down == a.round_down_rel(m)1862 1863 m = 0.5991864 m_up = 0.61865 m_down = 0.51866 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.251872 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]==b1878 assert len(a[key])==11879 1880 c = (2.1,1.)1881 a.append(c)1882 assert a[(2.,1.)][0]==c1883 assert a[(2.2,1.)][0]==c1884 1885 def test_mapped_Discretised_Tuple_Set(self):1886 1887 def map(sequence):1888 return [len(sequence)]1889 1890 tol = 0.51891 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.011898 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 a1904 1905 #___________end of Peters tests1906 1907 1722 def test_add_region_from_polygon(self): 1908 1723 m=Mesh()
Note: See TracChangeset
for help on using the changeset viewer.