303_LaplaceEquation.py 1.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677
  1. # Add the igl library to the modules search path
  2. import sys, os
  3. sys.path.insert(0, os.getcwd() + "/../")
  4. import pyigl as igl
  5. V = igl.eigen.MatrixXd()
  6. F = igl.eigen.MatrixXi()
  7. igl.readOFF("../../tutorial/shared/camelhead.off",V,F)
  8. # Find boundary edges
  9. E = igl.eigen.MatrixXi()
  10. igl.boundary_facets(F,E);
  11. # Find boundary vertices
  12. b = igl.eigen.MatrixXi()
  13. IA = igl.eigen.MatrixXi()
  14. IC = igl.eigen.MatrixXi()
  15. igl.unique(E,b,IA,IC);
  16. # List of all vertex indices
  17. vall = igl.eigen.MatrixXi()
  18. vin = igl.eigen.MatrixXi()
  19. igl.coloni(0,V.rows()-1,vall)
  20. # List of interior indices
  21. igl.setdiff(vall,b,vin,IA)
  22. # Construct and slice up Laplacian
  23. L = igl.eigen.SparseMatrixd()
  24. L_in_in = igl.eigen.SparseMatrixd()
  25. L_in_b = igl.eigen.SparseMatrixd()
  26. igl.cotmatrix(V,F,L)
  27. igl.slice(L,vin,vin,L_in_in)
  28. igl.slice(L,vin,b,L_in_b)
  29. # Dirichlet boundary conditions from z-coordinate
  30. bc = igl.eigen.MatrixXd()
  31. Z = V.col(2)
  32. igl.slice(Z,b,bc)
  33. # Solve PDE
  34. solver = igl.eigen.SimplicialLLTsparse(-L_in_in)
  35. Z_in = solver.solve(L_in_b*bc)
  36. # slice into solution
  37. igl.slice_into(Z_in,vin,Z)
  38. # Alternative, short hand
  39. mqwf = igl.min_quad_with_fixed_data()
  40. # Linear term is 0
  41. B = igl.eigen.MatrixXd()
  42. B.setZero(V.rows(),1);
  43. # Empty constraints
  44. Beq = igl.eigen.MatrixXd()
  45. Aeq = igl.eigen.SparseMatrixd()
  46. # Our cotmatrix is _negative_ definite, so flip sign
  47. igl.min_quad_with_fixed_precompute(-L,b,Aeq,True,mqwf)
  48. igl.min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z)
  49. # Pseudo-color based on solution
  50. C = igl.eigen.MatrixXd()
  51. igl.jet(Z,True,C)
  52. # Plot the mesh with pseudocolors
  53. viewer = igl.viewer.Viewer()
  54. viewer.data.set_mesh(V, F)
  55. viewer.core.show_lines = False
  56. viewer.data.set_colors(C)
  57. viewer.launch()