demo_eigenvalue.py 1.0 KB

1234567891011121314151617181920212223242526272829303132333435363738394041
  1. # A simple eigenvalue solver
  2. # ==========================
  3. from dolfin import *
  4. from vedo import download
  5. from vedo.dolfin import plot
  6. # Define mesh, function space
  7. fpath = download("https://vedo.embl.es/examples/data/box_with_dent.xml.gz")
  8. mesh = Mesh(fpath)
  9. V = FunctionSpace(mesh, "Lagrange", 1)
  10. # Define basis and bilinear form
  11. u = TrialFunction(V)
  12. v = TestFunction(V)
  13. a = dot(grad(u), grad(v))*dx
  14. # Assemble stiffness form
  15. A = PETScMatrix()
  16. assemble(a, tensor=A)
  17. # Create eigensolver
  18. eigensolver = SLEPcEigenSolver(A)
  19. # Compute all eigenvalues of A x = \lambda x
  20. print("Computing eigenvalues. This can take a minute.")
  21. eigensolver.solve()
  22. # Extract largest (first) eigenpair
  23. r, c, rx, cx = eigensolver.get_eigenpair(0)
  24. print("Largest eigenvalue: ", r)
  25. # Initialize function and assign eigenvector
  26. u = Function(V)
  27. u.vector()[:] = rx
  28. # plot eigenfunction on mesh as colored points (ps=point size)
  29. plot(u, mode='mesh', ps=12, cmap='gist_earth')
  30. #or as wireframe
  31. plot(u, mode='mesh', wireframe=True, cmap='magma')