Coverage for src / sdynpy / fileio / sdynpy_pdf3D.py: 8%

1# -*- coding: utf-8 -*- 

2""" 

3Created on Wed Jan 8 11:44:09 2025 

4 

5@author: dprohe 

6""" 

7 

8from ..core.sdynpy_geometry import (Geometry,CoordinateSystemArray,IGNORE_PLOTS, 

9 global_coord,_beam_elem_types,_face_element_types, 

10 _solid_element_types, _vtk_element_map, 

11 _element_types,split_list,GeometryPlotter) 

12from ..core.sdynpy_shape import ShapeArray 

13from ..core.sdynpy_colors import colormap 

14 

15import datetime 

16import warnings 

17import numpy as np 

18import pyvista as pv 

19 

20try: 

21 from vtk import vtkU3DExporter 

22except ImportError: 

23 vtkU3DExporter = None 

24 

25def _create_u3d_for_animation(output_name, geometry,node_size : int = 5, line_width: int = 1, opacity: float = 1.0, 

26 show_edges: bool =False): 

27 if IGNORE_PLOTS: 

28 raise ValueError('Cannot Create VTK Objects, No GPU Found') 

29 if vtkU3DExporter is None: 

30 raise ValueError('Cannot Import vtkU3DExporter. It must first be installed with `pip install vtk-u3dexporter`') 

31 

32 # Get part information 

33 nodes = geometry.node.flatten() 

34 css = geometry.coordinate_system.flatten() 

35 elems = geometry.element.flatten() 

36 tls = geometry.traceline.flatten() 

37 

38 coordinate_system = CoordinateSystemArray(nodes.shape) 

39 for key, node in nodes.ndenumerate(): 

40 coordinate_system[key] = css[np.where(geometry.coordinate_system.id == node.def_cs)] 

41 global_node_positions = global_coord(coordinate_system, nodes.coordinate) 

42 node_index_dict = {node.id[()]: index[0] for index, node in nodes.ndenumerate()} 

43 node_index_map = np.vectorize(node_index_dict.__getitem__) 

44 

45 # Now go through and get the element and line connectivity 

46 face_element_connectivity = [] 

47 face_element_colors = [] 

48 node_colors = [] 

49 line_connectivity = [] 

50 line_colors = [] 

51 for index, node in nodes.ndenumerate(): 

52 # element_connectivity.append(1) 

53 # element_connectivity.append(index[0]) 

54 # element_colors.append(node.color) 

55 node_colors.append(node.color) 

56 for index, element in elems.ndenumerate(): 

57 # Check which type of element it is 

58 if element.type in _beam_elem_types: # Beamlike element, use a line 

59 try: 

60 line_connectivity.append(node_index_map(element.connectivity)) 

61 except KeyError: 

62 raise KeyError( 

63 'Element {:} contains a node id not found in the node array'.format(element.id)) 

64 line_colors.append(element.color) 

65 elif element.type in _face_element_types: 

66 try: 

67 if len(element.connectivity) == 3: 

68 face_element_connectivity.append(node_index_map(element.connectivity)) 

69 face_element_colors.append(element.color) 

70 else: 

71 face_element_connectivity.append(node_index_map(element.connectivity[[0,1,3]])) 

72 face_element_colors.append(element.color) 

73 face_element_connectivity.append(node_index_map(element.connectivity[[1,2,3]])) 

74 face_element_colors.append(element.color) 

75 except KeyError: 

76 raise KeyError( 

77 'Element {:} contains a node id not found in the node array'.format(element.id)) 

78 elif element.type in _solid_element_types: 

79 warnings.warn('Solid Elements are currently not supported and will be skipped!') 

80 else: 

81 raise ValueError('Unknown element type {:}'.format(element.type)) 

82 

83 for index, tl in tls.ndenumerate(): 

84 for conn_group in split_list(tl.connectivity, 0): 

85 if len(conn_group) == 0: 

86 continue 

87 try: 

88 mapped_conn_group = node_index_map(conn_group) 

89 except KeyError: 

90 raise KeyError( 

91 'Traceline {:} contains a node id not found in the node array'.format(tl.id)) 

92 for indices in zip(mapped_conn_group[:-1],mapped_conn_group[1:]): 

93 line_connectivity.append(indices) 

94 line_colors.append(tl.color) 

95 

96 # Now we start to plot 

97 plotter = GeometryPlotter(editor=True) 

98 face_map = [] 

99 for conn,color in zip(face_element_connectivity,face_element_colors): 

100 node_indices,inverse = np.unique(conn,return_inverse=True) 

101 node_positions = global_node_positions[node_indices][inverse] 

102 node_connectivity = np.arange(node_positions.shape[0]) 

103 mesh = pv.PolyData(node_positions,faces=[len(node_connectivity)]+[c for c in node_connectivity]) 

104 mesh.cell_data['color'] = color 

105 mesh.SetObjectName('Elem: {:} {:} {:}'.format(*geometry.node.id[node_indices[inverse]])) 

106 plotter.add_mesh(mesh,scalars='color',cmap=colormap, clim = [0,15], 

107 show_edges=show_edges, show_scalar_bar = False, 

108 line_width=line_width,opacity=opacity,label='Elem: {:} {:} {:}'.format(*geometry.node.id[node_indices[inverse]])) 

109 face_map.append(geometry.node.id[node_indices[inverse]]) 

110 

111 line_map = [] 

112 for conn,color in zip(line_connectivity,line_colors): 

113 node_indices,inverse = np.unique(conn,return_inverse=True) 

114 node_positions = global_node_positions[node_indices][inverse] 

115 node_connectivity = np.arange(node_positions.shape[0]) 

116 mesh = pv.PolyData(node_positions,lines=[len(node_connectivity)]+[c for c in node_connectivity]) 

117 mesh.cell_data['color'] = color 

118 mesh.SetObjectName('Line: {:} {:}'.format(*geometry.node.id[node_indices[inverse]])) 

119 plotter.add_mesh(mesh,scalars='color',cmap=colormap, clim = [0,15], 

120 show_edges=show_edges, show_scalar_bar = False, 

121 line_width=line_width,opacity=opacity,label='Line: {:} {:}'.format(*geometry.node.id[node_indices[inverse]])) 

122 line_map.append(geometry.node.id[node_indices[inverse]]) 

123 

124 point_map = [] 

125 for node,position,color in zip(geometry.node.id,global_node_positions,geometry.node.color): 

126 if node_size > 0: 

127 mesh = pv.PolyData(position) 

128 mesh.cell_data['color'] = color 

129 plotter.add_mesh(mesh, scalars = color, cmap=colormap, clim=[0,15], 

130 show_edges=show_edges, show_scalar_bar=False, point_size=node_size, 

131 opacity=opacity) 

132 point_map.append(node) 

133 

134 exporter = vtkU3DExporter.vtkU3DExporter() 

135 exporter.SetFileName(output_name) 

136 exporter.SetInput(plotter.render_window) 

137 exporter.Write() 

138 

139 plotter.close() 

140 

141 return face_map, line_map, point_map 

142 

143def create_animated_modeshape_content(geometry,shapes = None,u3d_name='geometry', 

144 js_name = 'mode_shapes.js',one_js=True, 

145 node_size = 5, line_width = 2, 

146 opacity = 1, show_edges = False, 

147 debug_js = False, displacement_scale = 1.0): 

148 """ 

149 Creates files that can be embedded into a PDF to create interactive content 

150  

151 This function produces a Universal3D model containing Geometry information. 

152 This function also produces one or more JavaScript programs that are used 

153 by the PDF to animate the model data. The model and programs can be 

154 embedded into a PDF either by using Adobe Acrobat Professional, or by using 

155 LaTeX with the media9 package to build a PDF document. 

156 

157 Parameters 

158 ---------- 

159 geometry : Geometry 

160 A SDynPy Geometry object that will be exported to Universal3D format. 

161 shapes : ShapeArray, optional 

162 A SDynPy ShapeArray object that contains the mode shapes that will be 

163 used to animate the model in the final PDF. If not specified, only the 

164 geometry will be exported. 

165 u3d_name : str, optional 

166 The file path and name of the Universal3D file that will be created 

167 from the SDynPy Geometry. Note that the .u3d extension will be 

168 automatically appended to the file path. The default is 'geometry'. 

169 js_name : str, optional 

170 The file path and name of the JavaScript file(s) that will be created 

171 from the geometry and shape information. The final file name will be 

172 constructed from the code `js_name.format(i+1)` where `i` is the shape 

173 index. This allows for a different file name to be produced for each 

174 mode if `one_js` is `False`. The default is 'mode_shapes.js'. 

175 one_js : bool, optional 

176 If True, only one JavaScript file will be produced containing all mode 

177 shape information. If False, one JavaScript file will be produced per 

178 mode in the `shapes` argument, with each script having that shape as 

179 shape that is plotted first. Even if `one_js` is False, all JavaScript 

180 programs will have all modes in them. The default is True. 

181 node_size : int, optional 

182 The size of the rendered node in pixels. The default is 5. 

183 line_width : int, optional 

184 The width of lines rendered in the model. The default is 2. Note that 

185 lines with a width of 1 may not show their color well in the final PDF. 

186 opacity : float, optional 

187 The opacity of the model, from 0 (transparent) to 1 (opaque). The 

188 default is 1. 

189 show_edges : bool, optional 

190 If True, edges will be drawn on element faces. The default is False. 

191 debug_js : bool or str, optional 

192 If True, the JavaScript programs will have various print statements 

193 enabled to help users debug. If False, these will be commented out. 

194 String values of 'map', 'point', 'line', or 'face' can also be specified 

195 to only print debug text for those operations. The default is False. 

196 displacement_scale : float, optional 

197 A scale factor applied to the shape displacements to achieve a 

198 reasonable deflection size. The default is 1.0. 

199 

200 Notes 

201 ----- 

202 There is no support currently for Solid elements (e.g. tets or hexes). 

203 Only surface elements (tris and quads), lines (beam elements or tracelines), 

204 and points (nodes) are currently supported. 

205 

206 Returns 

207 ------- 

208 None. 

209 

210 """ 

211 face_map, line_map, point_map = _create_u3d_for_animation(u3d_name, 

212 geometry, 

213 node_size, 

214 line_width, 

215 opacity, 

216 show_edges) 

217 

218 if shapes is None: 

219 return 

220 

221 # Need to get the shape information in the global coordinate system 

222 node_displacement_dict = {} 

223 for node_id in geometry.node.id: 

224 node_displacement_dict[node_id] = np.zeros((shapes.size,3)) 

225 for shape_index,shape in enumerate(shapes.flatten()): 

226 coords = abs( 

227 shape.coordinate[ 

228 (shape.coordinate.node == node_id) 

229 & np.isin(abs(shape.coordinate.direction), [1, 2, 3]) 

230 ] 

231 ) 

232 shape_data = shape[coords] 

233 coord_directions = geometry.global_deflection(coords) 

234 global_deflection = np.sum(shape_data[:,np.newaxis]*coord_directions,axis=0) 

235 node_displacement_dict[node_id][shape_index] = global_deflection 

236 

237 node_position_dict = {node_id:position for node_id,position in zip(geometry.node.id,geometry.global_node_coordinate())} 

238 

239 if one_js: 

240 js_indices = [0] 

241 else: 

242 js_indices = np.arange(shapes.size) 

243 

244 for i in js_indices: 

245 js_output_name = js_name.format(i+1) 

246 with open(js_output_name,'w') as f: 

247 f.write( 

248 """//Written by sdynpy_pdf3D.py on {date:} 

249  

250 function pause(){{if(speed)lastspeed=speed; speed=0;}} 

251 function play(){{speed=lastspeed;}} 

252 function scaleSpeed(s){{lastspeed*=s; if(speed) speed=lastspeed;}} 

253 function scale(s){{scale*=s;}} 

254  

255 var omega0=Math.PI; // init. angular frequency (half turn per second) 

256 var speed=1; // speed multiplier 

257 var lastspeed=1; 

258 var theta=0; 

259 var scale={displacement_scale:}; 

260 var shape={shape:}; 

261 var num_shapes = {num_shapes:}; 

262 var panel = 1 

263 var node_positions = {{}}; 

264 var node_displacements = {{}}; 

265  

266 // Node and Shape Information 

267 """.format(date=datetime.datetime.now(), 

268 shape=i+1, 

269 num_shapes=shapes.size, 

270 displacement_scale=displacement_scale)) 

271 for node_id in node_position_dict: 

272 f.write('\nnode_positions[{:}] = Vector3({:},{:},{:});\n'.format(node_id,*node_position_dict[node_id])) 

273 f.write('node_displacements[{:}] = [\n'.format(node_id)) 

274 for shape_vector in node_displacement_dict[node_id]: 

275 f.write(' Vector3({:},{:},{:}),\n'.format(*shape_vector)) 

276 f.write(' ];\n') 

277 f.write('\n// Frequency Information\n') 

278 f.write('var mode_frequencies = [];\n') 

279 for i,shape in enumerate(shapes.flatten()): 

280 f.write('mode_frequencies[{:}] = {:0.2f};\n'.format(i,shape.frequency)) 

281 # Now we need to go through and collect the node information and 

282 # Map it back to node id numbers 

283 f.write('\n// Connectivity Maps\n') 

284 f.write('face_map = [\n') 

285 for face in face_map: 

286 f.write(' {:},\n'.format(str(list(face)))) 

287 f.write(' ];\n') 

288 f.write('line_map = [\n') 

289 for line in line_map: 

290 f.write(' {:},\n'.format(str(list(line)))) 

291 f.write(' ];\n') 

292 f.write('point_map = [\n') 

293 for point in point_map: 

294 f.write(' {:},\n'.format(point)) 

295 f.write(' ];\n') 

296 f.write(""" 

297 var scene = this.scene; 

298 var face_nodes = []; 

299 var line_nodes = []; 

300 var point_nodes = []; 

301 for (i=0; i < scene.nodes.count; i++){{ 

302 name = scene.nodes.getByIndex(i).name 

303 if (name.slice(0,4) === "Mesh"){{ 

304 {debug:}host.console.println("Mesh found at node "+i+" "+ name) 

305 face_nodes.push(i); 

306 }} else if (name.slice(0,4) === "Line"){{ 

307 {debug:}host.console.println("Line found at node "+i+" "+ name) 

308 line_nodes.push(i); 

309 }} else if (name.slice(0,4) === "Poin"){{ 

310 {debug:}host.console.println("Point found at node "+i+" "+ name) 

311 point_nodes.push(i); 

312 }} else {{ 

313 {debug:}host.console.println("Unknown node "+i+" "+name) 

314 }} 

315 }} 

316 """.format(debug='' if debug_js is True or debug_js=='map' else '// ')) 

317 

318 # Now we need to go through and apply displacements 

319 f.write(''' 

320 // Translate Nodes 

321 for (i=0; i < point_nodes.length; i++){{ 

322 geometry_node_id = point_map[i]; 

323 model_node = scene.nodes.getByIndex(point_nodes[i]); 

324 {debug:}host.console.println("Point Index "+i+" corresponds to model node "+model_node.name+" and geometry node id "+geometry_node_id); 

325 model_node.transform.setIdentity(); 

326 model_node.transform.translateInPlace(node_displacements[geometry_node_id][shape-1].scale(scale)); 

327 {debug:}host.console.println(model_node.transform.toString()+"\\n") 

328 }}'''.format(debug='' if debug_js is True or debug_js=='point' else '// ')) 

329 f.write(''' 

330 // Translate Lines 

331  

332 // Computation variables 

333 var p1; // position of node 1 

334 var p2; // position of node 2 

335 var d1; // displacement of node 1 

336 var d2; // displacement of node 2 

337 var v1; // traceline vector from p1 to p2 (p2-p1) 

338 var v2; // traceline vector after modification ((d2+p2)-(d1+p1)) 

339 var vc1; // traceline center ((p2+p1)/2) 

340 var vc2; // traceline center after modification ((p2+d2+p1+d1)/2) 

341 var axis; // axis of rotation v1 x v2 unit vector 

342 var angle; // angle of rotation cos(angle) = v1.v2/(|v1||v2|) 

343 var c; // scale factor |v2|/|v1| 

344 var t; // translate vc2 - c*vc1 

345  

346 // Go through traceline nodes 

347 for (i=0; i < line_nodes.length; i++){{ 

348 model_node = scene.nodes.getByIndex(line_nodes[i]); 

349 geometry_node_ids = line_map[i]; 

350 {debug:}host.console.println("Line Index "+i+" corresponds to model node "+model_node.name+" and geometry node ids "+geometry_node_ids); 

351 model_node.transform.setIdentity(); 

352 p1 = node_positions[geometry_node_ids[0]]; 

353 p2 = node_positions[geometry_node_ids[1]]; 

354 d1 = node_displacements[geometry_node_ids[0]][shape-1].scale(scale); 

355 d2 = node_displacements[geometry_node_ids[1]][shape-1].scale(scale); 

356 v1 = p2.add(p1.scale(-1)); 

357 v2 = p2.add(d2).add(p1.add(d1).scale(-1)); 

358 vc1 = p2.add(p1).scale(0.5); 

359 vc2 = p2.add(d2).add(p1).add(d1).scale(0.5); 

360 axis = v1.cross(v2); 

361 angle = Math.acos(v1.dot(v2)/(Math.sqrt(v1.dot(v1))*Math.sqrt(v2.dot(v2)))); 

362 c = Math.sqrt(v2.dot(v2)/v1.dot(v1)); 

363 t = vc2.add(vc1.scale(-c)); 

364 {debug:}host.console.println("Transforming TL: "+geometry_node_ids); 

365 {debug:}host.console.println(" p1: "+p1) 

366 {debug:}host.console.println(" p2: "+p2) 

367 {debug:}host.console.println(" d1: "+d1) 

368 {debug:}host.console.println(" d2: "+d2) 

369 {debug:}host.console.println(" v1: "+v1) 

370 {debug:}host.console.println(" v2: "+v2) 

371 {debug:}host.console.println(" vc1: "+vc1) 

372 {debug:}host.console.println(" vc2: "+vc2) 

373 {debug:}host.console.println(" axis: "+axis) 

374 {debug:}host.console.println(" angle: "+angle) 

375 {debug:}host.console.println(" c: "+c) 

376 {debug:}host.console.println(" t: "+t) 

377 {debug:}host.console.println("\\n Transform prior to anything:") 

378 {debug:}host.console.println(model_node.transform.toString()) 

379 //model_node.transform.rotateAboutLineInPlace(angle,vc1,vc1.add(axis)); //There seems to be a bug with this command... 

380 model_node.transform.translateInPlace(vc1.scale(-1)); 

381 model_node.transform.rotateAboutVectorInPlace(angle,axis); 

382 model_node.transform.translateInPlace(vc1); 

383 {debug:}host.console.println("\\n Transform after rotation:") 

384 {debug:}host.console.println(model_node.transform.toString()) 

385 model_node.transform.scaleInPlace(c,c,c); 

386 {debug:}host.console.println("\\n Transform after scaling:") 

387 {debug:}host.console.println(model_node.transform.toString()) 

388 model_node.transform.translateInPlace(t); 

389 {debug:}host.console.println("\\n Transform after translation:") 

390 {debug:}host.console.println(model_node.transform.toString()) 

391 }}'''.format(debug='' if debug_js is True or debug_js=='line' else '// ')) 

392 

393 f.write(''' 

394 // We want to compute the SVD of the transformation, so we will define a number of functions here 

395 var _gamma = 5.828427124; 

396 var _cstar = 0.923879532; 

397 var _sstar = 0.3826834323; 

398 var EPSILON = 1e-6; 

399  

400 function condSwap(c, X, Y) 

401 { 

402 // used in step 2 

403 // var Z = X; 

404 // X = c ? Y : X; 

405 // Y = c ? Z : Y; 

406 return (c ? [Y,X] : [X,Y]) 

407 } 

408  

409 function condNegSwap(c, X, Y) 

410 { 

411 // used in step 2 and 3 

412 // var Z = -X; 

413 // X = c ? Y : X; 

414 // Y = c ? Z : Y; 

415 return (c ? [Y,-X] : [X,Y]) 

416 } 

417  

418 // matrix multiplication M = A * B 

419 function multAB(a11, a12, a13,a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33) 

420 { 

421 var m11=a11*b11 + a12*b21 + a13*b31; var m12=a11*b12 + a12*b22 + a13*b32; var m13=a11*b13 + a12*b23 + a13*b33; 

422 var m21=a21*b11 + a22*b21 + a23*b31; var m22=a21*b12 + a22*b22 + a23*b32; var m23=a21*b13 + a22*b23 + a23*b33; 

423 var m31=a31*b11 + a32*b21 + a33*b31; var m32=a31*b12 + a32*b22 + a33*b32; var m33=a31*b13 + a32*b23 + a33*b33; 

424 return [m11,m12,m13,m21,m22,m23,m31,m32,m33]; 

425 } 

426  

427 // matrix multiplication M = Transpose[A] * B 

428 function multAtB(a11, a12, a13, 

429 a21, a22, a23, 

430 a31, a32, a33, 

431 b11, b12, b13, 

432 b21, b22, b23, 

433 b31, b32, b33) 

434 { 

435 var m11=a11*b11 + a21*b21 + a31*b31; var m12=a11*b12 + a21*b22 + a31*b32; var m13=a11*b13 + a21*b23 + a31*b33; 

436 var m21=a12*b11 + a22*b21 + a32*b31; var m22=a12*b12 + a22*b22 + a32*b32; var m23=a12*b13 + a22*b23 + a32*b33; 

437 var m31=a13*b11 + a23*b21 + a33*b31; var m32=a13*b12 + a23*b22 + a33*b32; var m33=a13*b13 + a23*b23 + a33*b33; 

438 return [m11,m12,m13,m21,m22,m23,m31,m32,m33]; 

439 } 

440  

441 function quatToMat3(x,y,z,w) 

442 { 

443 var qxx = x*x; 

444 var qyy = y*y; 

445 var qzz = z*z; 

446 var qxz = x*z; 

447 var qxy = x*y; 

448 var qyz = y*z; 

449 var qwx = w*x; 

450 var qwy = w*y; 

451 var qwz = w*z; 

452  

453 var m11=1 - 2*(qyy + qzz); var m12=2*(qxy - qwz); var m13=2*(qxz + qwy); 

454 var m21=2*(qxy + qwz); var m22=1 - 2*(qxx + qzz); var m23=2*(qyz - qwx); 

455 var m31=2*(qxz - qwy); var m32=2*(qyz + qwx); var m33=1 - 2*(qxx + qyy); 

456 return [m11,m12,m13,m21,m22,m23,m31,m32,m33]; 

457 } 

458  

459 function approximateGivensQuaternion(a11, a12, a22) 

460 { 

461 /* 

462 * Given givens angle computed by approximateGivensAngles, 

463 * compute the corresponding rotation quaternion. 

464 */ 

465 var ch = 2*(a11-a22); 

466 var sh = a12; 

467 var b = _gamma*sh*sh < ch*ch; 

468 var w = 1/Math.sqrt(ch*ch+sh*sh); 

469 ch=b?w*ch:_cstar; 

470 sh=b?w*sh:_sstar; 

471 return [ch,sh] 

472 } 

473  

474 function jacobiConjugation(x, y, z, s11, s21, s22, s31, s32, s33, qV) 

475 { 

476 var return_val = approximateGivensQuaternion(s11,s21,s22); 

477 var ch = return_val[0]; 

478 var sh = return_val[1]; 

479  

480 var scale = ch*ch+sh*sh; 

481 var a = (ch*ch-sh*sh)/scale; 

482 var b = (2*sh*ch)/scale; 

483  

484 // make temp copy of S 

485 var _s11 = s11; 

486 var _s21 = s21; var _s22 = s22; 

487 var _s31 = s31; var _s32 = s32; var _s33 = s33; 

488  

489 // perform conjugation S = Q''*S*Q 

490 // Q already implicitly solved from a, b 

491 s11 = a*(a*_s11 + b*_s21) + b*(a*_s21 + b*_s22); 

492 s21 = a*(-b*_s11 + a*_s21) + b*(-b*_s21 + a*_s22); s22=-b*(-b*_s11 + a*_s21) + a*(-b*_s21 + a*_s22); 

493 s31 = a*_s31 + b*_s32; s32=-b*_s31 + a*_s32; s33=_s33; 

494  

495 // update cumulative rotation qV 

496 var tmp = []; 

497 tmp[0]=qV[0]*sh; 

498 tmp[1]=qV[1]*sh; 

499 tmp[2]=qV[2]*sh; 

500 sh *= qV[3]; 

501  

502 qV[0] *= ch; 

503 qV[1] *= ch; 

504 qV[2] *= ch; 

505 qV[3] *= ch; 

506  

507 // (x,y,z) corresponds to ((0,1,2),(1,2,0),(2,0,1)) 

508 // for (p,q) = ((0,1),(1,2),(0,2)) 

509 qV[z] += sh; 

510 qV[3] -= tmp[z]; // w 

511 qV[x] += tmp[y]; 

512 qV[y] -= tmp[x]; 

513  

514 // re-arrange matrix for next iteration 

515 _s11 = s22; 

516 _s21 = s32; _s22 = s33; 

517 _s31 = s21; _s32 = s31; _s33 = s11; 

518 s11 = _s11; 

519 s21 = _s21; s22 = _s22; 

520 s31 = _s31; s32 = _s32; s33 = _s33; 

521  

522 return [s11,s21,s22,s31,s32,s33,qV] 

523  

524 } 

525  

526 function dist2(x, y, z) 

527 { 

528 return x*x+y*y+z*z; 

529 } 

530  

531 function jacobiEigenanlysis(s11,s21,s22,s31,s32,s33) 

532 { 

533 var return_val; 

534 var qV = [] 

535 qV[3]=1; qV[0]=0;qV[1]=0;qV[2]=0; // follow same indexing convention as GLM 

536 for (var i=0;i<4;i++) 

537 { 

538 // we wish to eliminate the maximum off-diagonal element 

539 // on every iteration, but cycling over all 3 possible rotations 

540 // in fixed order (p,q) = (1,2) , (2,3), (1,3) still retains 

541 // asymptotic convergence 

542 return_val = jacobiConjugation(0,1,2,s11,s21,s22,s31,s32,s33,qV); // p,q = 0,1 

543 s11 = return_val[0]; s21 = return_val[1]; s22 = return_val[2]; s31 = return_val[3]; 

544 s32 = return_val[4]; s33 = return_val[5]; qV = return_val[6]; 

545 return_val = jacobiConjugation(1,2,0,s11,s21,s22,s31,s32,s33,qV); // p,q = 1,2 

546 s11 = return_val[0]; s21 = return_val[1]; s22 = return_val[2]; s31 = return_val[3]; 

547 s32 = return_val[4]; s33 = return_val[5]; qV = return_val[6]; 

548 return_val = jacobiConjugation(2,0,1,s11,s21,s22,s31,s32,s33,qV); // p,q = 0,2 

549 s11 = return_val[0]; s21 = return_val[1]; s22 = return_val[2]; s31 = return_val[3]; 

550 s32 = return_val[4]; s33 = return_val[5]; qV = return_val[6]; 

551 } 

552 return [s11,s21,s22,s31,s32,s33,qV] 

553 } 

554  

555 function sortSingularValues(b11, b12, b13, 

556 b21, b22, b23, 

557 b31, b32, b33, 

558 v11, v12, v13, 

559 v21, v22, v23, 

560 v31, v32, v33) 

561 { 

562 var return_val; 

563 var rho1 = dist2(b11,b21,b31); 

564 var rho2 = dist2(b12,b22,b32); 

565 var rho3 = dist2(b13,b23,b33); 

566 var c; 

567 c = rho1 < rho2; 

568 return_val = condNegSwap(c,b11,b12); 

569 b11 = return_val[0]; b12 = return_val[1]; 

570 return_val = condNegSwap(c,v11,v12); 

571 v11 = return_val[0]; v12 = return_val[1]; 

572 return_val = condNegSwap(c,b21,b22);  

573 b21 = return_val[0]; b22 = return_val[1]; 

574 return_val = condNegSwap(c,v21,v22);  

575 v21 = return_val[0]; v22 = return_val[1]; 

576 return_val = condNegSwap(c,b31,b32);  

577 b31 = return_val[0]; b32 = return_val[1]; 

578 return_val = condNegSwap(c,v31,v32);  

579 v31 = return_val[0]; v32 = return_val[1]; 

580 return_val = condSwap(c,rho1,rho2);  

581 rho1 = return_val[0]; rho2 = return_val[1]; 

582 c = rho1 < rho3; 

583 return_val = condNegSwap(c,b11,b13);  

584 b11 = return_val[0]; b13 = return_val[1];  

585 return_val = condNegSwap(c,v11,v13);  

586 v11 = return_val[0]; v13 = return_val[1]; 

587 return_val = condNegSwap(c,b21,b23);  

588 b21 = return_val[0]; b23 = return_val[1];  

589 return_val = condNegSwap(c,v21,v23);  

590 v21 = return_val[0]; v23 = return_val[1]; 

591 return_val = condNegSwap(c,b31,b33);  

592 b31 = return_val[0]; b33 = return_val[1];  

593 return_val = condNegSwap(c,v31,v33);  

594 v31 = return_val[0]; v33 = return_val[1]; 

595 return_val = condSwap(c,rho1,rho3);  

596 rho1 = return_val[0]; rho3 = return_val[1]; 

597 c = rho2 < rho3; 

598 return_val = condNegSwap(c,b12,b13);  

599 b12 = return_val[0]; b13 = return_val[1];  

600 return_val = condNegSwap(c,v12,v13);  

601 v12 = return_val[0]; v13 = return_val[1]; 

602 return_val = condNegSwap(c,b22,b23);  

603 b22 = return_val[0]; b23 = return_val[1];  

604 return_val = condNegSwap(c,v22,v23);  

605 v22 = return_val[0]; v23 = return_val[1]; 

606 return_val = condNegSwap(c,b32,b33);  

607 b32 = return_val[0]; b33 = return_val[1];  

608 return_val = condNegSwap(c,v32,v33);  

609 v32 = return_val[0]; v33 = return_val[1]; 

610 return [b11, b12, b13, b21, b22, b23, b31, b32, b33, v11, v12, v13, v21, v22, v23, v31, v32, v33]; 

611 } 

612  

613 function QRGivensQuaternion(a1, a2) 

614 { 

615 // a1 = pivot point on diagonal 

616 // a2 = lower triangular entry we want to annihilate 

617 var rho = Math.sqrt(a1*a1 + a2*a2); 

618  

619 var sh = rho > EPSILON ? a2 : 0; 

620 var ch = Math.abs(a1) + Math.max(rho,EPSILON); 

621 var b = a1 < 0; 

622 var return_val = condSwap(b,sh,ch); 

623 sh = return_val[0]; ch = return_val[1]; 

624 var w = 1/Math.sqrt(ch*ch+sh*sh); 

625 ch *= w; 

626 sh *= w; 

627 return [ch,sh]; 

628 } 

629  

630 function QRDecomposition(b11, b12, b13, 

631 b21, b22, b23, 

632 b31, b32, b33) 

633 {  

634 var ch1,sh1,ch2,sh2,ch3,sh3; 

635 var a,b; 

636 var return_val; 

637  

638 // first givens rotation (ch,0,0,sh) 

639 return_val = QRGivensQuaternion(b11,b21); 

640 ch1 = return_val[0]; sh1 = return_val[1] 

641 a=1-2*sh1*sh1; 

642 b=2*ch1*sh1; 

643 // apply B = Q'' * B 

644 var r11=a*b11+b*b21; var r12=a*b12+b*b22; var r13=a*b13+b*b23; 

645 var r21=-b*b11+a*b21; var r22=-b*b12+a*b22; var r23=-b*b13+a*b23; 

646 var r31=b31; var r32=b32; var r33=b33; 

647  

648 // second givens rotation (ch,0,-sh,0) 

649 return_val = QRGivensQuaternion(r11,r31); 

650 ch2 = return_val[0]; sh2 = return_val[1]; 

651 a=1-2*sh2*sh2; 

652 b=2*ch2*sh2; 

653 // apply B = Q'' * B; 

654 b11=a*r11+b*r31; b12=a*r12+b*r32; b13=a*r13+b*r33; 

655 b21=r21; b22=r22; b23=r23; 

656 b31=-b*r11+a*r31; b32=-b*r12+a*r32; b33=-b*r13+a*r33; 

657  

658 // third givens rotation (ch,sh,0,0) 

659 return_val = QRGivensQuaternion(b22,b32); 

660 ch3 = return_val[0]; sh3 = return_val[1]; 

661 a=1-2*sh3*sh3; 

662 b=2*ch3*sh3; 

663 // R is now set to desired value 

664 r11=b11; r12=b12; r13=b13; 

665 r21=a*b21+b*b31; r22=a*b22+b*b32; r23=a*b23+b*b33; 

666 r31=-b*b21+a*b31; r32=-b*b22+a*b32; r33=-b*b23+a*b33; 

667  

668 // construct the cumulative rotation Q=Q1 * Q2 * Q3 

669 // the number of floating point operations for three quaternion multiplications 

670 // is more or less comparable to the explicit form of the joined matrix. 

671 // certainly more memory-efficient! 

672 var sh12=sh1*sh1; 

673 var sh22=sh2*sh2; 

674 var sh32=sh3*sh3; 

675  

676 var q11=(-1+2*sh12)*(-1+2*sh22);  

677 var q12=4*ch2*ch3*(-1+2*sh12)*sh2*sh3+2*ch1*sh1*(-1+2*sh32);  

678 var q13=4*ch1*ch3*sh1*sh3-2*ch2*(-1+2*sh12)*sh2*(-1+2*sh32); 

679  

680 var q21=2*ch1*sh1*(1-2*sh22);  

681 var q22=-8*ch1*ch2*ch3*sh1*sh2*sh3+(-1+2*sh12)*(-1+2*sh32);  

682 var q23=-2*ch3*sh3+4*sh1*(ch3*sh1*sh3+ch1*ch2*sh2*(-1+2*sh32)); 

683  

684 var q31=2*ch2*sh2;  

685 var q32=2*ch3*(1-2*sh22)*sh3;  

686 var q33=(-1+2*sh22)*(-1+2*sh32); 

687  

688 return [q11,q12,q13,q21,q22,q23,q31,q32,q33,r11,r12,r13,r21,r22,r23,r31,r32,r33]; 

689 } 

690  

691 function svd(a11, a12, a13, 

692 a21, a22, a23, 

693 a31, a32, a33) 

694 { 

695 var ata_return_val; 

696 var jacobi_return_val; 

697 var quat2mat_return_val; 

698 var b_return_vals; 

699 var ssv_return_vals; 

700  

701 // normal equations matrix 

702 var ATA11, ATA12, ATA13; 

703 var ATA21, ATA22, ATA23; 

704 var ATA31, ATA32, ATA33; 

705 // host.console.println("Multiplying Matrices ATA "); 

706 ata_return_val = multAtB(a11,a12,a13,a21,a22,a23,a31,a32,a33, 

707 a11,a12,a13,a21,a22,a23,a31,a32,a33); 

708 // host.console.println("result = " + ata_return_val); 

709 ATA11=ata_return_val[0]; 

710 ATA12=ata_return_val[1]; 

711 ATA13=ata_return_val[2]; 

712 ATA21=ata_return_val[3]; 

713 ATA22=ata_return_val[4]; 

714 ATA23=ata_return_val[5]; 

715 ATA31=ata_return_val[6]; 

716 ATA32=ata_return_val[7]; 

717 ATA33=ata_return_val[8]; 

718  

719 // symmetric eigenalysis 

720 // host.console.println("Jacobi Eigenanalysis"); 

721 jacobi_return_val = jacobiEigenanlysis( ATA11,ATA21,ATA22, ATA31,ATA32,ATA33); 

722 // host.console.println("result = " + jacobi_return_val); 

723 ATA11 = jacobi_return_val[0]; 

724 ATA21 = jacobi_return_val[1]; 

725 ATA22 = jacobi_return_val[2]; 

726 ATA31 = jacobi_return_val[3]; 

727 ATA32 = jacobi_return_val[4]; 

728 ATA33 = jacobi_return_val[5]; 

729 var qV = jacobi_return_val[6]; 

730  

731 // host.console.println("Quat2Mat3"); 

732 quat2mat_return_val = quatToMat3(qV[0],qV[1],qV[2],qV[3]); 

733 // host.console.println("result = " + jacobi_return_val); 

734 v11=quat2mat_return_val[0]; 

735 v12=quat2mat_return_val[1]; 

736 v13=quat2mat_return_val[2]; 

737 v21=quat2mat_return_val[3]; 

738 v22=quat2mat_return_val[4]; 

739 v23=quat2mat_return_val[5]; 

740 v31=quat2mat_return_val[6]; 

741 v32=quat2mat_return_val[7]; 

742 v33=quat2mat_return_val[8]; 

743  

744 var b11, b12, b13; 

745 var b21, b22, b23; 

746 var b31, b32, b33; 

747 b_return_vals = multAB(a11,a12,a13,a21,a22,a23,a31,a32,a33, 

748 v11,v12,v13,v21,v22,v23,v31,v32,v33); 

749 b11=b_return_vals[0];  

750 b12=b_return_vals[1];  

751 b13=b_return_vals[2];  

752 b21=b_return_vals[3];  

753 b22=b_return_vals[4];  

754 b23=b_return_vals[5];  

755 b31=b_return_vals[6];  

756 b32=b_return_vals[7];  

757 b33=b_return_vals[8]; 

758  

759 // sort singular values and find V 

760 ssv_return_vals = sortSingularValues(b11, b12, b13, b21, b22, b23, b31, b32, b33, 

761 v11,v12,v13,v21,v22,v23,v31,v32,v33); 

762 b11=ssv_return_vals[0]; 

763 b12=ssv_return_vals[1]; 

764 b13=ssv_return_vals[2]; 

765 b21=ssv_return_vals[3]; 

766 b22=ssv_return_vals[4]; 

767 b23=ssv_return_vals[5]; 

768 b31=ssv_return_vals[6]; 

769 b32=ssv_return_vals[7]; 

770 b33=ssv_return_vals[8]; 

771 v11=ssv_return_vals[9]; 

772 v12=ssv_return_vals[10]; 

773 v13=ssv_return_vals[11]; 

774 v21=ssv_return_vals[12]; 

775 v22=ssv_return_vals[13]; 

776 v23=ssv_return_vals[14]; 

777 v31=ssv_return_vals[15]; 

778 v32=ssv_return_vals[16]; 

779 v33=ssv_return_vals[17]; 

780  

781 // QR decomposition 

782 qr_return_vals = QRDecomposition(b11, b12, b13, b21, b22, b23, b31, b32, b33); 

783 u11=qr_return_vals[0]; 

784 u12=qr_return_vals[1]; 

785 u13=qr_return_vals[2]; 

786 u21=qr_return_vals[3]; 

787 u22=qr_return_vals[4]; 

788 u23=qr_return_vals[5]; 

789 u31=qr_return_vals[6]; 

790 u32=qr_return_vals[7]; 

791 u33=qr_return_vals[8]; 

792 s11=qr_return_vals[9]; 

793 s12=qr_return_vals[10]; 

794 s13=qr_return_vals[11]; 

795 s21=qr_return_vals[12]; 

796 s22=qr_return_vals[13]; 

797 s23=qr_return_vals[14]; 

798 s31=qr_return_vals[15]; 

799 s32=qr_return_vals[16]; 

800 s33=qr_return_vals[17]; 

801  

802 return [u11, u12, u13, 

803 u21, u22, u23, 

804 u31, u32, u33, 

805 s11, s12, s13, 

806 s21, s22, s23, 

807 s31, s32, s33, 

808 v11, v12, v13, 

809 v21, v22, v23, 

810 v31, v32, v33]; 

811 } 

812  

813 function mat_inv(n00,n01,n02,n10,n11,n12,n20,n21,n22) 

814 { 

815 var v1 = (n11 * n22 - n21 * n12); 

816 var v2 = (n10 * n22 - n12 * n20); 

817 var v3 = (n10 * n21 - n11 * n20); 

818  

819 var det = (n00 * v1) - (n01 * v2) + (n02 * v3); 

820  

821 var invdet = 1/det; 

822  

823 var ninv00 = v1 * invdet; 

824 var ninv01 = (n02 * n21 - n01 * n22) * invdet; 

825 var ninv02 = (n01 * n12 - n02 * n11) * invdet; 

826 var ninv10 = -v2 * invdet; 

827 var ninv11 = (n00 * n22 - n02 * n20) * invdet; 

828 var ninv12 = (n10 * n02 - n00 * n12) * invdet; 

829 var ninv20 = v3 * invdet; 

830 var ninv21 = (n20 * n01 - n00 * n21) * invdet; 

831 var ninv22 = (n00 * n11 - n10 * n01) * invdet; 

832  

833 return [ninv00,ninv01,ninv02,ninv10,ninv11,ninv12,ninv20,ninv21,ninv22] 

834 } 

835  

836 function AmultAT3x4(a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34) 

837 { 

838 var AAT11 = a11*a11 + a12*a12 + a13*a13 + a14*a14; 

839 var AAT12 = a11*a21 + a12*a22 + a13*a23 + a14*a24; 

840 var AAT13 = a11*a31 + a12*a32 + a13*a33 + a14*a34; 

841 var AAT21 = a11*a21 + a12*a22 + a13*a23 + a14*a24; 

842 var AAT22 = a21*a21 + a22*a22 + a23*a23 + a24*a24; 

843 var AAT23 = a21*a31 + a22*a32 + a23*a33 + a24*a34; 

844 var AAT31 = a11*a31 + a12*a32 + a13*a33 + a14*a34; 

845 var AAT32 = a21*a31 + a22*a32 + a23*a33 + a24*a34; 

846 var AAT33 = a31*a31 + a32*a32 + a33*a33 + a34*a34; 

847  

848 return [AAT11,AAT12,AAT13,AAT21,AAT22,AAT23,AAT31,AAT32,AAT33]; 

849 } 

850  

851 function matmul4x3x3(a11,a12,a13,a21,a22,a23,a31,a32,a33,a41,a42,a43, 

852 b11,b12,b13,b21,b22,b23,b31,b32,b33) 

853 { 

854 var C11 = a11*b11 + a12*b21 + a13*b31; 

855 var C12 = a11*b12 + a12*b22 + a13*b32; 

856 var C13 = a11*b13 + a12*b23 + a13*b33; 

857 var C21 = a21*b11 + a22*b21 + a23*b31; 

858 var C22 = a21*b12 + a22*b22 + a23*b32; 

859 var C23 = a21*b13 + a22*b23 + a23*b33; 

860 var C31 = a31*b11 + a32*b21 + a33*b31; 

861 var C32 = a31*b12 + a32*b22 + a33*b32; 

862 var C33 = a31*b13 + a32*b23 + a33*b33; 

863 var C41 = a41*b11 + a42*b21 + a43*b31; 

864 var C42 = a41*b12 + a42*b22 + a43*b32; 

865 var C43 = a41*b13 + a42*b23 + a43*b33; 

866  

867 return [C11,C12,C13,C21,C22,C23,C31,C32,C33,C41,C42,C43]; 

868 } 

869  

870 function matrix_axis_angle(a,b,c,d,e,f,g,h,i) 

871 { 

872 var axis = Vector3(h-f,c-g,d-b); 

873 var angle = Math.atan2(axis.length,(a+e+i-1)); 

874 axis = axis.scale(1/axis.length); 

875 return [axis,angle] 

876 } 

877 ''') 

878 

879 f.write(''' 

880 // Go through element nodes 

881 for (i=0; i < face_nodes.length; i++){{ 

882 model_node = scene.nodes.getByIndex(face_nodes[i]); 

883 geometry_node_ids = face_map[i]; 

884 {debug:}host.console.println("Element Index "+i+" corresponds to model node "+model_node.name+" and geometry node ids "+geometry_node_ids); 

885 model_node.transform.setIdentity(); 

886 p1 = node_positions[geometry_node_ids[0]]; 

887 p2 = node_positions[geometry_node_ids[1]]; 

888 p3 = node_positions[geometry_node_ids[2]]; 

889 centroid = p1.add(p2.add(p3)).scale(.3333333333333333); 

890 {debug:}host.console.println("\\n Initial Transformation:"); 

891 {debug:}host.console.println(model_node.transform.toString()); 

892 model_node.transform.translateInPlace(centroid.scale(-1)); 

893 {debug:}host.console.println("\\n Transform after removing centroid "+centroid+":"); 

894 {debug:}host.console.println(model_node.transform.toString()); 

895 p1 = p1.subtract(centroid) 

896 p2 = p2.subtract(centroid) 

897 p3 = p3.subtract(centroid) 

898 d1 = node_displacements[geometry_node_ids[0]][shape-1].scale(scale); 

899 d2 = node_displacements[geometry_node_ids[1]][shape-1].scale(scale); 

900 d3 = node_displacements[geometry_node_ids[2]][shape-1].scale(scale); 

901 pd1 = p1.add(d1); 

902 pd2 = p2.add(d2); 

903 pd3 = p3.add(d3); 

904 {debug:}host.console.println("Transforming Element: "+geometry_node_ids); 

905 {debug:}host.console.println(" p1: "+p1) 

906 {debug:}host.console.println(" p2: "+p2) 

907 {debug:}host.console.println(" p3: "+p3) 

908 {debug:}host.console.println(" pd1: "+pd1) 

909 {debug:}host.console.println(" pd2: "+pd2) 

910 {debug:}host.console.println(" pd3: "+pd3) 

911 var AAT = AmultAT3x4(p1.x,p1.y,p1.z,1,p2.x,p2.y,p2.z,1,p3.x,p3.y,p3.z,1); 

912 {debug:}host.console.println(" AAT: "+AAT) 

913 var pinvA = matmul4x3x3.apply(this,[p1.x,p2.x,p3.x,p1.y,p2.y,p3.y,p1.z,p2.z,p3.z,1,1,1].concat(mat_inv.apply(this,AAT))); 

914 {debug:}host.console.println(" pinvA: "+pinvA) 

915 var b = [pd1.x,pd1.y,pd1.z,pd2.x,pd2.y,pd2.z,pd3.x,pd3.y,pd3.z]; 

916 {debug:}host.console.println(" b: "+b) 

917 var x = matmul4x3x3.apply(this,pinvA.concat(b)); 

918 {debug:}host.console.println(" x: "+x) 

919 T = [x[0],x[3],x[6],x[1],x[4],x[7],x[2],x[5],x[8]]; 

920 t = Vector3(x[9],x[10],x[11]); 

921 {debug:}host.console.println(" T: "+T) 

922 {debug:}host.console.println(" t: "+t) 

923 svd_res = svd.apply(this,T); 

924 U = svd_res.slice(0,9); 

925 S = svd_res.slice(9,18); 

926 V = svd_res.slice(18); 

927 {debug:}host.console.println(" U: "+U); 

928 {debug:}host.console.println(" S: "+S); 

929 {debug:}host.console.println(" V: "+V); 

930 uvar = matrix_axis_angle.apply(this,U); 

931 uaxis = uvar[0]; uangle = uvar[1]; 

932 vvar = matrix_axis_angle.apply(this,V); 

933 vaxis = vvar[0]; vangle = vvar[1]; 

934 {debug:}host.console.println(" U axis: "+uaxis); 

935 {debug:}host.console.println(" U angle: "+uangle); 

936 {debug:}host.console.println(" V axis: "+vaxis); 

937 {debug:}host.console.println(" V angle: "+vangle); 

938 model_node.transform.rotateAboutVectorInPlace(-vangle,vaxis); 

939 {debug:}host.console.println("\\n Transform after rotation about "+vaxis+" by "+(-vangle)+":"); 

940 {debug:}host.console.println(model_node.transform.toString()); 

941 model_node.transform.scaleInPlace(S[0],S[4],1); 

942 {debug:}host.console.println("\\n Transform after scaling by "+[S[0],S[4],1]+":"); 

943 {debug:}host.console.println(model_node.transform.toString()); 

944 model_node.transform.rotateAboutVectorInPlace(uangle,uaxis); 

945 {debug:}host.console.println("\\n Transform after rotation about "+uaxis+" by "+(-uangle)+":"); 

946 {debug:}host.console.println(model_node.transform.toString()); 

947 model_node.transform.translateInPlace(t.add(centroid)); 

948 {debug:}host.console.println("\\n Transform after translation by t and centroid "+t.add(centroid)+":"); 

949 {debug:}host.console.println(model_node.transform.toString()); 

950  

951 }} 

952 '''.format(debug='' if debug_js is True or debug_js=='face' else '// ')) 

953 

954 # Note: there should be some way to build the transformation matrix without an 

955 # SVD. However, it results in one of the dimensions of the transformation 

956 # getting squashed to zero, which doesn't have any affect on the resulting 

957 # geometry, but it does mess with the shading due to the surface normal being 

958 # kind of screwed up. I'm not sure exactly the solution apart from computing an 

959 # SVD and making all singular values positive (which is exactly what I do now) 

960 # function matmul4x4x4( 

961 # a00, a01, a02, a03, 

962 # a10, a11, a12, a13, 

963 # a20, a21, a22, a23, 

964 # a30, a31, a32, a33, 

965 # b00, b01, b02, b03, 

966 # b10, b11, b12, b13, 

967 # b20, b21, b22, b23, 

968 # b30, b31, b32, b33 

969 # ) { 

970 # const result = [ 

971 # a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30, 

972 # a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31, 

973 # a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32, 

974 # a00 * b03 + a01 * b13 + a02 * b23 + a03 * b33, 

975 # a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30, 

976 # a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31, 

977 # a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32, 

978 # a10 * b03 + a11 * b13 + a12 * b23 + a13 * b33, 

979 # a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30, 

980 # a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31, 

981 # a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32, 

982 # a20 * b03 + a21 * b13 + a22 * b23 + a23 * b33, 

983 # a30 * b00 + a31 * b10 + a32 * b20 + a33 * b30, 

984 # a30 * b01 + a31 * b11 + a32 * b21 + a33 * b31, 

985 # a30 * b02 + a31 * b12 + a32 * b22 + a33 * b32, 

986 # a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33, 

987 # ]; 

988 

989 # return result; 

990 # } 

991 # // Try the same thing by just building the transformation without SVDs 

992 # T_raw = [1,0,0,0, //indices 0-3 

993 # 0,1,0,0, //indices 4-7 

994 # 0,0,1,0, //indices 8-11 

995 # 0,0,0,1]; //indices 12-15 

996 # {debug:}host.console.println("\\n Initial Transformation:"); 

997 # {debug:}host.console.println(T_raw); 

998 # T_raw[12] -= centroid.x; 

999 # T_raw[13] -= centroid.y; 

1000 # T_raw[14] -= centroid.z; 

1001 # {debug:}host.console.println("\\n With Removal of Centroid:"); 

1002 # {debug:}host.console.println(T_raw); 

1003 # X = [x[0],x[1],x[2],0, 

1004 # x[3],x[4],x[5],0, 

1005 # x[6],x[7],x[8],0, 

1006 # x[9],x[10],x[11],1]; 

1007 # T_raw = matmul4x4x4.apply(this,T_raw.concat(X)); 

1008 # {debug:}host.console.println("\\n After Rotation, Scaling, and Translation:"); 

1009 # {debug:}host.console.println(T_raw); 

1010 # T_raw[12] += centroid.x; 

1011 # T_raw[13] += centroid.y; 

1012 # T_raw[14] += centroid.z; 

1013 # {debug:}host.console.println("\\n After Addition of Centroid:"); 

1014 # {debug:}host.console.println(T_raw); 

1015 # model_node.transform.set(T_raw); 

1016 # break; 

1017 

1018 f.write(''' 

1019 scene.update(); 

1020  

1021  

1022 timeEvHnd=new TimeEventHandler(); 

1023 timeEvHnd.onEvent=function(event) { 

1024 var dtheta=speed*omega0*event.deltaTime; 

1025 if (dtheta!=0){ 

1026 ''') 

1027 

1028 f.write(''' 

1029 // Translate Nodes 

1030 for (i=0; i < point_nodes.length; i++){{ 

1031 geometry_node_id = point_map[i]; 

1032 model_node = scene.nodes.getByIndex(point_nodes[i]); 

1033 model_node.transform.setIdentity(); 

1034 model_node.transform.translateInPlace(node_displacements[geometry_node_id][shape-1].scale(scale*Math.cos(theta))); 

1035 }}'''.format()) 

1036 f.write(''' 

1037 // Translate and Rotate Lines 

1038 for (i=0; i < line_nodes.length; i++){{ 

1039 model_node = scene.nodes.getByIndex(line_nodes[i]); 

1040 geometry_node_ids = line_map[i]; 

1041 model_node.transform.setIdentity(); 

1042 p1 = node_positions[geometry_node_ids[0]]; 

1043 p2 = node_positions[geometry_node_ids[1]]; 

1044 d1 = node_displacements[geometry_node_ids[0]][shape-1].scale(scale*Math.cos(theta)); 

1045 d2 = node_displacements[geometry_node_ids[1]][shape-1].scale(scale*Math.cos(theta)); 

1046 v1 = p2.add(p1.scale(-1)); 

1047 v2 = p2.add(d2).add(p1.add(d1).scale(-1)); 

1048 vc1 = p2.add(p1).scale(0.5); 

1049 vc2 = p2.add(d2).add(p1).add(d1).scale(0.5); 

1050 axis = v1.cross(v2); 

1051 angle = Math.acos(v1.dot(v2)/(Math.sqrt(v1.dot(v1))*Math.sqrt(v2.dot(v2)))); 

1052 c = Math.sqrt(v2.dot(v2)/v1.dot(v1)); 

1053 t = vc2.add(vc1.scale(-c)); 

1054 //model_node.transform.rotateAboutLineInPlace(angle,vc1,vc1.add(axis)); //There seems to be a bug with this command... 

1055 if (angle > 0.00175) {{ 

1056 model_node.transform.translateInPlace(vc1.scale(-1)); 

1057 model_node.transform.rotateAboutVectorInPlace(angle,axis); 

1058 model_node.transform.translateInPlace(vc1); 

1059 }} 

1060 model_node.transform.scaleInPlace(c,c,c); 

1061 model_node.transform.translateInPlace(t); 

1062 }}'''.format()) 

1063 f.write(''' 

1064 // Go through element nodes 

1065 for (i=0; i < face_nodes.length; i++){{ 

1066 model_node = scene.nodes.getByIndex(face_nodes[i]); 

1067 geometry_node_ids = face_map[i]; 

1068 model_node.transform.setIdentity(); 

1069 p1 = node_positions[geometry_node_ids[0]]; 

1070 p2 = node_positions[geometry_node_ids[1]]; 

1071 p3 = node_positions[geometry_node_ids[2]]; 

1072 centroid = p1.add(p2.add(p3)).scale(.3333333333333333); 

1073 model_node.transform.translateInPlace(centroid.scale(-1)); 

1074 p1 = p1.subtract(centroid) 

1075 p2 = p2.subtract(centroid) 

1076 p3 = p3.subtract(centroid) 

1077 d1 = node_displacements[geometry_node_ids[0]][shape-1].scale(scale*Math.cos(theta)); 

1078 d2 = node_displacements[geometry_node_ids[1]][shape-1].scale(scale*Math.cos(theta)); 

1079 d3 = node_displacements[geometry_node_ids[2]][shape-1].scale(scale*Math.cos(theta)); 

1080 pd1 = p1.add(d1); 

1081 pd2 = p2.add(d2); 

1082 pd3 = p3.add(d3); 

1083 var AAT = AmultAT3x4(p1.x,p1.y,p1.z,1,p2.x,p2.y,p2.z,1,p3.x,p3.y,p3.z,1); 

1084 var pinvA = matmul4x3x3.apply(this,[p1.x,p2.x,p3.x,p1.y,p2.y,p3.y,p1.z,p2.z,p3.z,1,1,1].concat(mat_inv.apply(this,AAT))); 

1085 var b = [pd1.x,pd1.y,pd1.z,pd2.x,pd2.y,pd2.z,pd3.x,pd3.y,pd3.z]; 

1086 var x = matmul4x3x3.apply(this,pinvA.concat(b));