Commit 70f486ac by Jeremy BLEYER

 # # .. _demo_poisson_equation: # # Poisson equation # ================ # # This demo is implemented in a single Python file, # :download:demo_poisson.py, which contains both the variational forms # and the solver. # # This demo illustrates how to: # # * Solve a linear partial differential equation # * Create and apply Dirichlet boundary conditions # * Define Expressions # * Define a FunctionSpace # * Create a SubDomain # # The solution for :math:u in this demo will look as follows: # # .. image:: poisson_u.png # :scale: 75 % # # # Equation and problem definition # ------------------------------- # # The Poisson equation is the canonical elliptic partial differential # equation. For a domain :math:\Omega \subset \mathbb{R}^n with # boundary :math:\partial \Omega = \Gamma_{D} \cup \Gamma_{N}, the # Poisson equation with particular boundary conditions reads: # # .. math:: # - \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\ # u &= 0 \quad {\rm on} \ \Gamma_{D}, \\ # \nabla u \cdot n &= g \quad {\rm on} \ \Gamma_{N}. \\ # # Here, :math:f and :math:g are input data and :math:n denotes the # outward directed boundary normal. The most standard variational form # of Poisson equation reads: find :math:u \in V such that # # .. math:: # a(u, v) = L(v) \quad \forall \ v \in V, # # where :math:V is a suitable function space and # # .. math:: # a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\ # L(v) &= \int_{\Omega} f v \, {\rm d} x # + \int_{\Gamma_{N}} g v \, {\rm d} s. # # The expression :math:a(u, v) is the bilinear form and :math:L(v) # is the linear form. It is assumed that all functions in :math:V # satisfy the Dirichlet boundary conditions (:math:u = 0 \ {\rm on} \ # \Gamma_{D}). # # In this demo, we shall consider the following definitions of the input # functions, the domain, and the boundaries: # # * :math:\Omega = [0,1] \times [0,1] (a unit square) # * :math:\Gamma_{D} = \{(0, y) \cup (1, y) \subset \partial \Omega\} # (Dirichlet boundary) # * :math:\Gamma_{N} = \{(x, 0) \cup (x, 1) \subset \partial \Omega\} # (Neumann boundary) # * :math:g = \sin(5x) (normal derivative) # * :math:f = 10\exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02) (source # term) # # # Implementation # -------------- # # This description goes through the implementation (in # :download:demo_poisson.py) of a solver for the above described # Poisson equation step-by-step. # # First, the :py:mod:dolfin module is imported: :: from dolfin import * # We begin by defining a mesh of the domain and a finite element # function space :math:V relative to this mesh. As the unit square is # a very standard domain, we can use a built-in mesh provided by the # class :py:class:UnitSquareMesh . In order # to create a mesh consisting of 32 x 32 squares with each square # divided into two triangles, we do as follows :: # Create mesh and define function space mesh = UnitSquareMesh(32, 32) V = FunctionSpace(mesh, "Lagrange", 1) # The second argument to :py:class:FunctionSpace #  is the finite element # family, while the third argument specifies the polynomial # degree. Thus, in this case, our space V consists of first-order, # continuous Lagrange finite element functions (or in order words, # continuous piecewise linear polynomials). # # Next, we want to consider the Dirichlet boundary condition. A simple # Python function, returning a boolean, can be used to define the # subdomain for the Dirichlet boundary condition (:math:\Gamma_D). The # function should return True for those points inside the subdomain # and False for the points outside. In our case, we want to say that # the points :math:(x, y) such that :math:x = 0 or :math:x = 1 are # inside on the inside of :math:\Gamma_D. (Note that because of # rounding-off errors, it is often wise to instead specify :math:x < # \epsilon or :math:x > 1 - \epsilon where :math:\epsilon is a # small number (such as machine precision).) :: # Define Dirichlet boundary (x = 0 or x = 1) def boundary(x): return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS # Now, the Dirichlet boundary condition can be created using the class # :py:class:DirichletBC . A # :py:class:DirichletBC  takes three # arguments: the function space the boundary condition applies to, the # value of the boundary condition, and the part of the boundary on which # the condition applies. In our example, the function space is V, # the value of the boundary condition (0.0) can represented using a # :py:class:Constant  and the # Dirichlet boundary is defined immediately above. The definition of the # Dirichlet boundary condition then looks as follows: :: # Define boundary condition u0 = Constant(0.0) bc = DirichletBC(V, u0, boundary) # Next, we want to express the variational problem. First, we need to # specify the trial function :math:u and the test function :math:v, # both living in the function space :math:V. We do this by defining a # :py:class:TrialFunction  # and a :py:class:TestFunction #  on the previously defined # :py:class:FunctionSpace #  V. # # Further, the source :math:f and the boundary normal derivative # :math:g are involved in the variational forms, and hence we must # specify these. Both :math:f and :math:g are given by simple # mathematical formulas, and can be easily declared using the # :py:class:Expression  class. # Note that the strings defining f and g use C++ syntax since, # for efficiency, DOLFIN will generate and compile C++ code for these # expressions at run-time. # # With these ingredients, we can write down the bilinear form a and # the linear form L (using UFL operators). In summary, this reads :: # Define variational problem u = TrialFunction(V) v = TestFunction(V) f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2) g = Expression("sin(5*x[0])", degree=2) a = inner(grad(u), grad(v))*dx L = f*v*dx + g*v*ds # Now, we have specified the variational forms and can consider the # solution of the variational problem. First, we need to define a # :py:class:Function  u to # represent the solution. (Upon initialization, it is simply set to the # zero function.) A :py:class:Function #  represents a function living in # a finite element function space. Next, we can call the :py:func:solve #  function with the arguments a == L, # u and bc as follows: :: # Compute solution u = Function(V) solve(a == L, u, bc) # The function u will be modified during the call to solve. The # default settings for solving a variational problem have been # used. However, the solution process can be controlled in much more # detail if desired. # # A :py:class:Function  can be # manipulated in various ways, in particular, it can be plotted and # saved to file. Here, we output the solution to a VTK file (using # the suffix .pvd) for later visualization and also plot it using # the :py:func:plot  command: :: # Save solution in VTK format file = File("poisson.pvd") file << u # Plot solution import matplotlib.pyplot as plt plot(u) plt.show()