Package meshpy :: Module tools
[hide private]
[frames] | no frames]

Source Code for Module meshpy.tools

1 -def uniform_refine_triangles(points, elements, factor=2):
2 new_points = points[:] 3 new_elements = [] 4 old_face_to_new_faces = {} 5 face_point_dict = {} 6 7 points_per_edge = factor+1 8 9 def get_refined_face(a, b): 10 if a > b: 11 a, b = b, a 12 flipped = True 13 else: 14 flipped = False 15 16 try: 17 face_points = face_point_dict[a,b] 18 except KeyError: 19 a_pt, b_pt = [points[idx] for idx in [a, b]] 20 dx = (b_pt - a_pt)/factor 21 22 # build subdivided facet 23 face_points = [a] 24 25 for i in range(1, points_per_edge-1): 26 face_points.append(len(new_points)) 27 new_points.append(a_pt + dx*i) 28 29 face_points.append(b) 30 31 face_point_dict[a, b] = face_points 32 33 # build old_face_to_new_faces 34 old_face_to_new_faces[frozenset([a,b])] = [ 35 (face_points[i], face_points[i+1]) 36 for i in range(factor)] 37 38 if flipped: 39 return face_points[::-1] 40 else: 41 return face_points
42 43 for a, b, c in elements: 44 a_pt, b_pt, c_pt = [points[idx] for idx in [a, b, c]] 45 dr = (b_pt - a_pt)/factor 46 ds = (c_pt - a_pt)/factor 47 48 ab_refined, bc_refined, ac_refined = [ 49 get_refined_face(*pt_indices) 50 for pt_indices in [(a,b), (b,c), (a,c)]] 51 52 el_point_dict = {} 53 54 # fill out edges of el_point_dict 55 for i in range(points_per_edge): 56 el_point_dict[i,0] = ab_refined[i] 57 el_point_dict[0,i] = ac_refined[i] 58 el_point_dict[points_per_edge-1-i,i] = bc_refined[i] 59 60 # fill out interior of el_point_dict 61 for i in range(1, points_per_edge-1): 62 for j in range(1, points_per_edge-1-i): 63 el_point_dict[i,j] = len(new_points) 64 new_points.append(a_pt + dr*i + ds*j) 65 66 # generate elements 67 for i in range(0, points_per_edge-1): 68 for j in range(0, points_per_edge-1-i): 69 new_elements.append(( 70 el_point_dict[i,j], 71 el_point_dict[i+1,j], 72 el_point_dict[i,j+1], 73 )) 74 if i+1+j+1 <= factor: 75 new_elements.append(( 76 el_point_dict[i+1,j+1], 77 el_point_dict[i+1,j], 78 el_point_dict[i,j+1], 79 )) 80 81 from meshpy.triangle import MeshInfo 82 mi = MeshInfo() 83 mi.set_points(new_points) 84 mi.elements.resize(len(new_elements)) 85 for i, el in enumerate(new_elements): 86 mi.elements[i] = el 87 from meshpy.triangle import write_gnuplot_mesh 88 write_gnuplot_mesh("mesh.dat", mi) 89 90 return new_points, new_elements, old_face_to_new_faces 91 92 93 94
95 -def make_swizzle_matrix(spec):
96 import numpy 97 axes = ["x", "y", "z"] 98 99 mapping = dict((axis, axis) for axis in axes) 100 for one_spec in spec.split(","): 101 import_axis, final_axis = one_spec.split(":") 102 mapping[import_axis] = final_axis 103 104 assert set(mapping.keys()) == set(axes), \ 105 "axis mapping not complete" 106 assert set(axis.lstrip("-") for axis in mapping.itervalues()) == set(axes), \ 107 "Axis mapping not onto" 108 109 n = len(axes) 110 result = numpy.zeros((n, n), dtype=int) 111 112 for imp_axis, final_axis in mapping.iteritems(): 113 imp_axis = axes.index(imp_axis) 114 115 sign = 1 116 while final_axis.startswith("-"): 117 sign *= -1 118 final_axis = final_axis[1:] 119 final_axis = axes.index(final_axis) 120 121 result[final_axis, imp_axis] = sign 122 123 return result
124