poisson.py 1.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354
  1. """Poisson equation with Dirichlet conditions.
  2. -Laplace(u) = f in the unit square
  3. u = uD on the boundary
  4. uD = 1 + x^2 + 2*y^2
  5. (f = -6)
  6. """
  7. ########################################################### fenics
  8. import numpy as np
  9. from fenics import *
  10. # Create mesh and define function space
  11. mesh = UnitSquareMesh(8, 8)
  12. V = FunctionSpace(mesh, "P", 1)
  13. # Define boundary condition
  14. uD = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
  15. bc = DirichletBC(V, uD, "on_boundary")
  16. # Define variational problem
  17. w = TrialFunction(V)
  18. v = TestFunction(V)
  19. u = Function(V)
  20. f = Constant(-6.0)
  21. # Compute solution
  22. solve( dot(grad(w), grad(v))*dx == f*v*dx, u, bc)
  23. f = r'-\nabla^{2} u=f'
  24. ########################################################### vedo
  25. from vedo.dolfin import plot
  26. from vedo.pyplot import histogram
  27. from vedo import Latex
  28. l = Latex(f, s=0.2, c='k3').shift([.3,.6,.1])
  29. plot(u, l, cmap='jet', scalarbar='h', text=__doc__).clear()
  30. # Now show uD values on the boundary of a much finer mesh
  31. bmesh = BoundaryMesh(UnitSquareMesh(80, 80), "exterior")
  32. plot(uD, bmesh, cmap='cool', ps=5, legend='boundary') # ps = point size
  33. # now make some nonsense plot with the same plot() function
  34. yvals = u.compute_vertex_values(mesh)
  35. xvals = np.arange(len(yvals))
  36. plt = plot(xvals, yvals, 'go-')
  37. plt.show(new=True)
  38. # and a histogram
  39. hst = histogram(yvals)
  40. hst.show(new=True)