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
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
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
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
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