v0.13.1
print_rollers.py
Go to the documentation of this file.
1import sys
2import os
3from numpy.core.numeric import full
4import pyvista as pv
5import numpy as np
6import argparse
7
8from numpy.testing._private.utils import integer_repr
9
10MAX_NB_OF_ROLLERS = 10
11
13 x = txt.split(",")
14 nb_list = []
15 for i in x:
16 nb_list.append(float(i))
17 if len(nb_list) != 3:
18 exit("coords and orientation require 3 numbers: x,y,z")
19 return nb_list
20
21
22def is_not_vtk(file):
23 return not file.endswith('vtk')
24
25
27
28 a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 /
29 np.linalg.norm(vec2)).reshape(3)
30 v = np.cross(a, b)
31 c = np.dot(a, b)
32 s = np.linalg.norm(v)
33 if s <= 1e-12:
34 return np.array([0, 0, 0])
35 kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
36 R = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2))
37
38 sy = np.sqrt(R[0, 0] * R[0, 0] + R[1, 0] * R[1, 0])
39 angle0 = np.arctan2(R[2, 1], R[2, 2]) * 180 / np.pi
40 angle1 = np.arctan2(-R[2, 0], sy) * 180 / np.pi
41 angle2 = np.arctan2(R[1, 0], R[0, 0]) * 180 / np.pi
42
43 return np.array([angle0, angle1, angle2])
44
45def get_capped_cone(origin, orientation, radius1, angle, height, resolution):
46
47 orientation_np = np.array(orientation)
48 scale = np.linalg.norm(orientation_np)
49 norm_orientation = orientation_np / scale
50
51 full_h = radius1 / np.tan(angle * np.pi / 180)
52 new_origin = origin + norm_orientation * (full_h - height) / 2.
53 body_test = pv.Cone(center=new_origin, direction=orientation, height=full_h, radius=radius1, capping=False, angle=None, resolution=resolution)
54 cut_origin = origin + norm_orientation * (height / 2.)
55 slice = body_test.slice(normal=orientation,origin=cut_origin, generate_triangles=True)
56 slice = pv.PolyData(slice.points).delaunay_2d()
57 body_test = body_test.clip(normal=orientation, origin=cut_origin, invert=True)
58 body = pv.MultiBlock([body_test, slice]).combine()
59
60 return body
61
63 def __init__(self):
64 self.body_type = 'plane'
65 self.orientation = [0, 0, 1]
66 self.radius1 = 5
67 self.inner_radius = 1
68 self.height = 1
69 self.angle = 45
70 self.angle_a = 45
71 self.angle_b = 45
72 self.resolution = 64
73
74 def set_body_data_from_vertex(self, body_type_id, orientation_vec, roller_data):
75 types = {0: 'plane',
76 1: 'sphere',
77 2: 'cylinder',
78 3: 'cone',
79 4: 'torus',
80 5: 'roller',
81 6: 'stl',
82 7: 'nurbs',
83 8: 'last'
84 }
85
86 self.body_type = types[body_type_id]
87 self.orientation = orientation_vec
88 self.radius1 = roller_data[0]
89 self.inner_radius = roller_data[1]
90 # sometimes we use infinite cones or cylinders
91 self.height = roller_data[2] if not np.isinf(roller_data[2]) else self.radius1 * 3
92 self.angle = roller_data[3]
93 self.angle_a = roller_data[4]
94 self.angle_b = roller_data[5]
95
96 def get_body(self, origin):
97 # origin=[0, 0, 0]
98 orientation_np = np.array(self.orientation)
99 scale = np.linalg.norm(orientation_np)
100 normz_orien = orientation_np / scale
101
102 body = pv.Box()
103
104 if(self.body_type == 'cylinder'):
105 body = pv.Cylinder(center=origin, direction=self.orientation, radius=self.radius1,
106 height=self.height, resolution=self.resolution)
107 elif(self.body_type == 'sphere'):
108 body = pv.Sphere(center=origin, radius=self.radius1, phi_resolution=self.resolution)
109 elif (self.body_type == 'torus'):
110 pyvista_radius = self.radius1 - self.inner_radius
111 body1 = pv.ParametricTorus(ringradius=pyvista_radius, crosssectionradius=self.inner_radius)
112 body2 = pv.Cylinder(center=[0,0,0], direction=[0, 0, 1], radius=pyvista_radius,
113 height=self.inner_radius*2, resolution=self.resolution)
114 roller = [body1, body2]
115 body = pv.MultiBlock(roller).combine()
116
117
118 angles = get_euler_angles_from_vectors(np.array([0., 0., 1.]), normz_orien)
119 # print(angles[0], angles[1], angles[2])
120 body.rotate_x(angles[0])
121 body.rotate_y(angles[1])
122 body.rotate_z(angles[2])
123
124 transform_matrix = np.array(
125 [[1, 0, 0, origin[0]], [0, 1, 0, origin[1]], [0, 0, 1, origin[2]], [0, 0, 0, 1]])
126 body.transform(transform_matrix, inplace=True)
127
128 elif(self.body_type == 'plane'):
129 body = pv.Plane(center=origin, direction=self.orientation, i_resolution=self.resolution,
130 j_resolution=self.resolution, i_size=scale, j_size=scale)
131
132 elif(self.body_type == 'cone'):
133 body = get_capped_cone(origin, self.orientation, self.radius1, self.angle, self.height, self.resolution)
134 # body.plot(color="tan", smooth_shading=True)
135 # body.save("test_mesh_0.vtu")
136 elif(self.body_type == 'roller'):
137 # torus
138 pyvista_radius = self.radius1 - self.inner_radius
139 body1 = pv.ParametricTorus(ringradius=pyvista_radius, crosssectionradius=self.inner_radius)
140 angles = get_euler_angles_from_vectors(np.array([0., 0., 1.]), normz_orien)
141 body1.rotate_x(angles[0])
142 body1.rotate_y(angles[1])
143 body1.rotate_z(angles[2])
144 transform_matrix = np.array(
145 [[1, 0, 0, origin[0]], [0, 1, 0, origin[1]], [0, 0, 1, origin[2]], [0, 0, 0, 1]])
146 body1.transform(transform_matrix, inplace=True)
147
148 # cone1
149 slope_a = np.tan(self.angle_a * np.pi / 180)
150 offset_a = (self.inner_radius * slope_a / np.sqrt(1 + slope_a * slope_a)) + self.height / 2
151 new_origin = origin + normz_orien * offset_a
152 cone_rad = pyvista_radius + self.inner_radius / np.sqrt(1 + slope_a * slope_a)
153 body2 = get_capped_cone(new_origin, self.orientation, cone_rad, self.angle_a, self.height, self.resolution)
154
155 # cone2
156 slope_b = np.tan(self.angle_b * np.pi / 180)
157 offset_b = (self.inner_radius * slope_b / np.sqrt(1 + slope_b * slope_b)) + self.height / 2
158 new_origin = origin - normz_orien * offset_b
159 cone_rad = pyvista_radius + self.inner_radius / np.sqrt(1 + slope_b * slope_b)
160 orientation_flip = normz_orien * -1
161 body3 = get_capped_cone(new_origin, orientation_flip, cone_rad, self.angle_b, self.height, self.resolution)
162 body = pv.MultiBlock([body1, body2, body3]).combine()
163 else:
164 exit("this body type is not supported")
165
166 # body.plot()
167 return body
168
169
170
171if __name__ == '__main__':
172
173 # print("Number of arguments:", len(sys.argv), "arguments.")
174 # print("Argument List:", str(sys.argv))
175 resolution = 64
176 rigid_bodies_dict = {}
177
178 if "-resolution" in sys.argv:
179 resolution = int(sys.argv[sys.argv.index("-resolution") + 1])
180
181 # for x in range(MAX_NB_OF_ROLLERS):
182
183 # body_arg = "-body" + str(x)
184 # if body_arg in sys.argv:
185 # rigid_body = RigidBody()
186 # rigid_body.resolution = resolution
187 # rigid_body.body_type = sys.argv[sys.argv.index(body_arg) + 1]
188
189 # orientation_str = '-direction' + str(x)
190 # if orientation_str in sys.argv:
191 # orient_text = sys.argv[sys.argv.index(orientation_str) + 1]
192 # rigid_body.orientation = split_string_to_vector3(orient_text)
193
194 # if rigid_body.body_type == 'torus' or rigid_body.body_type == 'cylinder' or rigid_body.body_type == 'sphere' or rigid_body.body_type == 'cone' or rigid_body.body_type == 'roller':
195 # if "-radius" + str(x) in sys.argv:
196 # rigid_body.radius1 = float(
197 # sys.argv[sys.argv.index("-radius" + str(x)) + 1])
198 # else:
199 # exit("-radius x must be specified for cylinder, torus or sphere")
200 # if rigid_body.body_type == 'cone' or rigid_body.body_type == 'cylinder':
201 # if "-height" + str(x) in sys.argv:
202 # rigid_body.height = float(
203 # sys.argv[sys.argv.index("-height" + str(x)) + 1])
204 # if rigid_body.body_type == 'torus':
205 # if "-inner_radius" + str(x) in sys.argv:
206 # rigid_body.inner_radius = float(
207 # sys.argv[sys.argv.index("-inner_radius" + str(x)) + 1])
208 # else:
209 # exit(
210 # "-inner_radius x must be specified for cylinder, torus or sphere")
211 # if rigid_body.body_type == 'cone':
212 # if "-angle" + str(x) in sys.argv:
213 # rigid_body.angle = float(
214 # sys.argv[sys.argv.index("-angle" + str(x)) + 1])
215 # else:
216 # exit(
217 # "-angle x must be specified for cone")
218 # if rigid_body.body_type == 'roller':
219 # if "-angle_a" + str(x) in sys.argv:
220 # rigid_body.angle_a = float(
221 # sys.argv[sys.argv.index("-angle_a" + str(x)) + 1])
222 # else:
223 # exit(
224 # "-angle_a must be specified for roller")
225 # if "-angle_b" + str(x) in sys.argv:
226 # rigid_body.angle_b = float(
227 # sys.argv[sys.argv.index("-angle_b" + str(x)) + 1])
228 # else:
229 # exit(
230 # "-angle_b must be specified for roller")
231 # if "-fillet" + str(x) in sys.argv:
232 # rigid_body.inner_radius = float(
233 # sys.argv[sys.argv.index("-fillet" + str(x)) + 1])
234 # rigid_body.height = 3 * rigid_body.inner_radius
235 # else:
236 # exit(
237 # "-fillet must be specified for roller")
238 # if "-height" + str(x) in sys.argv:
239 # rigid_body.height = float(
240 # sys.argv[sys.argv.index("-height" + str(x)) + 1])
241 # rigid_bodies_dict[x] = rigid_body
242
243 # print("Number of bodies defined: " + str(len(rigid_bodies_dict)) + ".")
244 parser = argparse.ArgumentParser(
245 description="Convert multiple vtk roller files to vtk with rollers geometry", add_help=False)
246 parser.add_argument(
247 "-file_name", help="list of vtk files or a regexp mask", nargs='+')
248 args, unknown = parser.parse_known_args()
249
250 args.file_name.sort(key=os.path.getmtime)
251
252 # my_files = ["/Users/karollewandowski/mofem_install/users_modules_debug/multifield_plasticity/out_roller_2.vtk"]
253 # args.file_name = my_files
254
255 for file in args.file_name:
256 file_vtu = os.path.splitext(file)[0] + ".vtu"
257 if os.path.exists(file_vtu):
258 if os.path.getmtime(file) < os.path.getmtime(file_vtu):
259 print("skipping {} file".format(file))
260 continue
261 print("Processing " + file + " file.")
262 if not file.endswith('vtk'):
263 exit("Provide vtk file or files (-file_name)")
264 mesh = pv.read(file)
265
266 # if mesh.n_points != len(rigid_bodies_dict):
267 # print("vtk file does not have enough data points")
268 # continue
269
270 index = 0
271 all_bodies = []
272 for point in mesh.points:
273
274 id_array = mesh['_ROLLER_ID'][index]
275 body_type_id = mesh['_BODY_TYPE'][index]
276 orientation_vec = mesh['_ORIENTATION'][index]
277 roller_data = mesh['_ROLLER_DATA'][index]
278
279 rigid_body = RigidBody()
280 rigid_body.resolution = resolution
281 rigid_body.set_body_data_from_vertex(body_type_id, orientation_vec, roller_data);
282
283 # if id_array not in rigid_bodies_dict:
284 # exit("vtk files and input have inconsistent ids")
285 all_bodies.append(rigid_body.get_body(point))
286 index = index + 1
287
288 mesh = pv.MultiBlock(all_bodies).combine()
289 mesh.save(file[:-4] + ".vtu")
290 exit()
if(!static_cast< bool >(ifstream(param_file)))
def get_body(self, origin)
def set_body_data_from_vertex(self, body_type_id, orientation_vec, roller_data)
def is_not_vtk(file)
def split_string_to_vector3(txt)
def get_euler_angles_from_vectors(vec1, vec2)
def get_capped_cone(origin, orientation, radius1, angle, height, resolution)