.. Documentation for the singular Poisson demo from DOLFIN. .. _demo_singular_poisson_cpp_documentation: Singular Poisson ================ This demo is implemented in a single Python file, :download:`demo_singular.py`, which contains both the variational forms and the solver. .. include:: ../common.txt Implementation -------------- This description goes through the implementation (in :download:`demo_singular.py`) of a solver for the above described Poisson equation step-by-step. First, the :py:mod:`dolfin` module is imported: .. code-block:: python from dolfin import * Then, we check that dolfin is configured with the backend called PETSc, since it provides us with a wide range of methods used by :py:class:`KrylovSolver `. We set PETSc as our backend for linear algebra. .. code-block:: python # Test for PETSc if not has_linear_algebra_backend("PETSc"): info("DOLFIN has not been configured with PETSc. Exiting.") exit() parameters["linear_algebra_backend"] = "PETSc" We begin by defining a mesh of the domain and a finite element function space :math:`V` relative to this mesh. We use a built-in mesh provided by the class :py:class:`UnitSquareMesh `. In order to create a mesh consisting of :math:`64 \times 64` squares with each square divided into two triangles, we do as follows: .. code-block:: python # Create mesh and define function space mesh = UnitSquareMesh(64, 64) V = FunctionSpace(mesh, "CG", 1) Now, we need to specify the trial functions (the unknowns) and the test functions on the space :math:`V`. This can be done using a :py:class:`TrialFunction ` and a :py:class:`TestFunction ` as follows: .. code-block:: python u = TrialFunction(V) v = TestFunction(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. .. code-block:: python f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") g = Expression("-sin(5*x[0])") With :math:`u,v,f` and :math:`g`, we can write down the bilinear form :math:`a` and the linear form :math:`L` (using UFL operators). .. code-block:: python a = inner(grad(u), grad(v))*dx L = f*v*dx + g*v*ds In order to transform our variational problem into a linear system we need to assemble the coefficient matrix ``A`` and the right-side vector ``b``. We do this using the function :py:meth:`assemble `: .. code-block:: python # Assemble system A = assemble(a) b = assemble(L) We specify a Vector for storing the result by defining a :py:class:`Function `. .. code-block:: python # Solution Function u = Function(V) Next, we specify the iterative solver we want to use, in this case a :py:class:`KrylovSolver `. The first argument is the left-hand side matrix, and the second argument specifies the method used. In this case we use the Generalized Minimum Residual (GMRES) method. .. code-block:: python # Create Krylov solver solver = KrylovSolver(A, "gmres") We impose our additional constraint by removing the null space component from the solution vector. In order to do this we need a basis for the null space. This is done by creating a vector that spans the null space, and then defining a basis from it. The basis is then attached to the :py:class:`KrylovSolver ` as its null space. .. code-block:: python # Create vector that spans the null space null_vec = Vector(u.vector()) V.dofmap().set(null_vec, 1.0) null_vec *= 1.0/null_vec.norm("l2") # Create null space basis object and attach to Krylov solver null_space = VectorSpaceBasis([null_vec]) solver.set_nullspace(null_space) Orthogonalization of ``b`` with respect to the null space makes sure that it doesn't contain any component in the null space. .. code-block:: python null_space.orthogonalize(b); Finally we are able to solve our linear system .. code-block:: python solver.solve(u.vector(), b) and plot the solution .. code-block:: python plot(u, interactive=True) Complete code ------------- .. literalinclude:: demo_singular.py :start-after: # Begin demo