v0.13.0
print_rollers.py
Go to the documentation of this file.
1 import sys
2 import os
3 from numpy.core.numeric import full
4 import pyvista as pv
5 import numpy as np
6 import argparse
7 
8 from numpy.testing._private.utils import integer_repr
9 
10 MAX_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 
22 def 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 
45 def 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 
62 class RigidBody:
63  def __init__(self):
64  self.body_typebody_type = 'plane'
65  self.orientationorientation = [0, 0, 1]
66  self.radius1radius1 = 5
67  self.inner_radiusinner_radius = 1
68  self.heightheight = 1
69  self.angleangle = 45
70  self.angle_aangle_a = 45
71  self.angle_bangle_b = 45
72  self.resolutionresolution = 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_typebody_type = types[body_type_id]
87  self.orientationorientation = orientation_vec
88  self.radius1radius1 = roller_data[0]
89  self.inner_radiusinner_radius = roller_data[1]
90  # sometimes we use infinite cones or cylinders
91  self.heightheight = roller_data[2] if not np.isinf(roller_data[2]) else self.radius1radius1 * 3
92  self.angleangle = roller_data[3]
93  self.angle_aangle_a = roller_data[4]
94  self.angle_bangle_b = roller_data[5]
95 
96  def get_body(self, origin):
97  # origin=[0, 0, 0]
98  orientation_np = np.array(self.orientationorientation)
99  scale = np.linalg.norm(orientation_np)
100  normz_orien = orientation_np / scale
101 
102  body = pv.Box()
103 
104  if(self.body_typebody_type == 'cylinder'):
105  body = pv.Cylinder(center=origin, direction=self.orientationorientation, radius=self.radius1radius1,
106  height=self.heightheight, resolution=self.resolutionresolution)
107  elif(self.body_typebody_type == 'sphere'):
108  body = pv.Sphere(center=origin, radius=self.radius1radius1, phi_resolution=self.resolutionresolution)
109  elif (self.body_typebody_type == 'torus'):
110  pyvista_radius = self.radius1radius1 - self.inner_radiusinner_radius
111  body1 = pv.ParametricTorus(ringradius=pyvista_radius, crosssectionradius=self.inner_radiusinner_radius)
112  body2 = pv.Cylinder(center=[0,0,0], direction=[0, 0, 1], radius=pyvista_radius,
113  height=self.inner_radiusinner_radius*2, resolution=self.resolutionresolution)
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_typebody_type == 'plane'):
129  body = pv.Plane(center=origin, direction=self.orientationorientation, i_resolution=self.resolutionresolution,
130  j_resolution=self.resolutionresolution, i_size=scale, j_size=scale)
131 
132  elif(self.body_typebody_type == 'cone'):
133  body = get_capped_cone(origin, self.orientationorientation, self.radius1radius1, self.angleangle, self.heightheight, self.resolutionresolution)
134  # body.plot(color="tan", smooth_shading=True)
135  # body.save("test_mesh_0.vtu")
136  elif(self.body_typebody_type == 'roller'):
137  # torus
138  pyvista_radius = self.radius1radius1 - self.inner_radiusinner_radius
139  body1 = pv.ParametricTorus(ringradius=pyvista_radius, crosssectionradius=self.inner_radiusinner_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_aangle_a * np.pi / 180)
150  offset_a = (self.inner_radiusinner_radius * slope_a / np.sqrt(1 + slope_a * slope_a)) + self.heightheight / 2
151  new_origin = origin + normz_orien * offset_a
152  cone_rad = pyvista_radius + self.inner_radiusinner_radius / np.sqrt(1 + slope_a * slope_a)
153  body2 = get_capped_cone(new_origin, self.orientationorientation, cone_rad, self.angle_aangle_a, self.heightheight, self.resolutionresolution)
154 
155  # cone2
156  slope_b = np.tan(self.angle_bangle_b * np.pi / 180)
157  offset_b = (self.inner_radiusinner_radius * slope_b / np.sqrt(1 + slope_b * slope_b)) + self.heightheight / 2
158  new_origin = origin - normz_orien * offset_b
159  cone_rad = pyvista_radius + self.inner_radiusinner_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_bangle_b, self.heightheight, self.resolutionresolution)
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 
171 if __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 split_string_to_vector3(txt)
def get_capped_cone(origin, orientation, radius1, angle, height, resolution)
def get_euler_angles_from_vectors(vec1, vec2)
def is_not_vtk(file)