305_QuadraticProgramming.py 2.0 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. import igl
  2. b = igl.eigen.MatrixXi()
  3. B = igl.eigen.MatrixXd()
  4. bc = igl.eigen.MatrixXd()
  5. lx = igl.eigen.MatrixXd()
  6. ux = igl.eigen.MatrixXd()
  7. Beq = igl.eigen.MatrixXd()
  8. Bieq = igl.eigen.MatrixXd()
  9. Z = igl.eigen.MatrixXd()
  10. Q = igl.eigen.SparseMatrixd()
  11. Aeq = igl.eigen.SparseMatrixd()
  12. Aieq = igl.eigen.SparseMatrixd()
  13. def solve(viewer):
  14. global Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,Z
  15. params = igl.active_set_params()
  16. params.max_iter = 8
  17. igl.active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,params,Z)
  18. C = igl.eigen.MatrixXd()
  19. igl.jet(Z,0,1,C)
  20. viewer.data.set_colors(C)
  21. def key_down(viewer, key, mod):
  22. global Beq,solve
  23. if key == ord('.'):
  24. Beq[0,0] = Beq[0,0] * 2.0
  25. solve(viewer)
  26. return True
  27. elif key == ord(','):
  28. Beq[0,0] = Beq[0,0] / 2.0
  29. solve(viewer)
  30. return True
  31. elif key == ord(' '):
  32. solve(viewer)
  33. return True
  34. return False;
  35. V = igl.eigen.MatrixXd()
  36. F = igl.eigen.MatrixXi()
  37. igl.readOFF("../tutorial/shared/cheburashka.off",V,F)
  38. # Plot the mesh
  39. viewer = igl.viewer.Viewer()
  40. viewer.data.set_mesh(V, F)
  41. viewer.core.show_lines = False
  42. viewer.callback_key_down = key_down
  43. # One fixed point on belly
  44. b = igl.eigen.MatrixXi([[2556]])
  45. bc = igl.eigen.MatrixXd([[1]])
  46. # Construct Laplacian and mass matrix
  47. L = igl.eigen.SparseMatrixd()
  48. M = igl.eigen.SparseMatrixd()
  49. Minv = igl.eigen.SparseMatrixd()
  50. igl.cotmatrix(V,F,L)
  51. igl.massmatrix(V,F,igl.MASSMATRIX_TYPE_VORONOI,M);
  52. igl.invert_diag(M,Minv)
  53. # Bi-Laplacian
  54. Q = L.transpose() * (Minv * L)
  55. # Zero linear term
  56. B = igl.eigen.MatrixXd.Zero(V.rows(),1)
  57. # Lower and upper bound
  58. lx = igl.eigen.MatrixXd.Zero(V.rows(),1)
  59. ux = igl.eigen.MatrixXd.Ones(V.rows(),1)
  60. # Equality constraint constrain solution to sum to 1
  61. Beq = igl.eigen.MatrixXd([[0.08]])
  62. Aeq = M.diagonal().transpose().sparseView()
  63. # (Empty inequality constraints)
  64. solve(viewer)
  65. print("Press '.' to increase scale and resolve.")
  66. print("Press ',' to decrease scale and resolve.")
  67. viewer.launch()