303_LaplaceEquation.py 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889
  1. # This file is part of libigl, a simple c++ geometry processing library.
  2. #
  3. # Copyright (C) 2017 Sebastian Koch <s.koch@tu-berlin.de> and Daniele Panozzo <daniele.panozzo@gmail.com>
  4. #
  5. # This Source Code Form is subject to the terms of the Mozilla Public License
  6. # v. 2.0. If a copy of the MPL was not distributed with this file, You can
  7. # obtain one at http://mozilla.org/MPL/2.0/.
  8. import sys, os
  9. # Add the igl library to the modules search path
  10. sys.path.insert(0, os.getcwd() + "/../")
  11. import pyigl as igl
  12. from shared import TUTORIAL_SHARED_PATH, check_dependencies
  13. dependencies = ["viewer"]
  14. check_dependencies(dependencies)
  15. V = igl.eigen.MatrixXd()
  16. F = igl.eigen.MatrixXi()
  17. igl.readOFF(TUTORIAL_SHARED_PATH + "camelhead.off", V, F)
  18. # Find boundary edges
  19. E = igl.eigen.MatrixXi()
  20. igl.boundary_facets(F, E)
  21. # Find boundary vertices
  22. b = igl.eigen.MatrixXi()
  23. IA = igl.eigen.MatrixXi()
  24. IC = igl.eigen.MatrixXi()
  25. igl.unique(E, b, IA, IC)
  26. # List of all vertex indices
  27. vall = igl.eigen.MatrixXi()
  28. vin = igl.eigen.MatrixXi()
  29. igl.coloni(0, V.rows() - 1, vall)
  30. # List of interior indices
  31. igl.setdiff(vall, b, vin, IA)
  32. # Construct and slice up Laplacian
  33. L = igl.eigen.SparseMatrixd()
  34. L_in_in = igl.eigen.SparseMatrixd()
  35. L_in_b = igl.eigen.SparseMatrixd()
  36. igl.cotmatrix(V, F, L)
  37. igl.slice(L, vin, vin, L_in_in)
  38. igl.slice(L, vin, b, L_in_b)
  39. # Dirichlet boundary conditions from z-coordinate
  40. bc = igl.eigen.MatrixXd()
  41. Z = V.col(2)
  42. igl.slice(Z, b, bc)
  43. # Solve PDE
  44. solver = igl.eigen.SimplicialLLTsparse(-L_in_in)
  45. Z_in = solver.solve(L_in_b * bc)
  46. # slice into solution
  47. igl.slice_into(Z_in, vin, Z)
  48. # Alternative, short hand
  49. mqwf = igl.min_quad_with_fixed_data()
  50. # Linear term is 0
  51. B = igl.eigen.MatrixXd()
  52. B.setZero(V.rows(), 1)
  53. # Empty constraints
  54. Beq = igl.eigen.MatrixXd()
  55. Aeq = igl.eigen.SparseMatrixd()
  56. # Our cotmatrix is _negative_ definite, so flip sign
  57. igl.min_quad_with_fixed_precompute(-L, b, Aeq, True, mqwf)
  58. igl.min_quad_with_fixed_solve(mqwf, B, bc, Beq, Z)
  59. # Pseudo-color based on solution
  60. C = igl.eigen.MatrixXd()
  61. igl.jet(Z, True, C)
  62. # Plot the mesh with pseudocolors
  63. viewer = igl.viewer.Viewer()
  64. viewer.data.set_mesh(V, F)
  65. viewer.core.show_lines = False
  66. viewer.data.set_colors(C)
  67. viewer.launch()