.. Documentation for the tensor weighted Poisson demo from DOLFIN.

.. _demo_pde_tensor-weighted-poisson_python_documentation:

Tensor-weighted Poisson
=======================

This demo is implemented in two files; one file,
:download:`generate_data.py` , for generating data, and one file,
:download:`demo_tensor-weighted-poisson.py` , which contains both the
vaiational form and the solver.

.. include:: ../common.txt

Implementation
--------------

Implementation of generate_data.py
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

First, the :py:mod:`dolfin` module is imported:

.. code-block:: python

    from dolfin import *

Then, we define a mesh of the domain. We use the built-in mesh,
provided by the class :py:class:`UnitSquareMesh
<dolfin.cpp.mesh.UnitSquareMesh>`. In order to create a mesh
consisting of :math:`32 \times 32` squares with each square divided
into two triangles, we do as follows:

.. code-block:: python

    # Create mesh
    mesh = UnitSquareMesh(32, 32)

Now, we create mesh functions to store the values of the conductivity
matrix as it varies over the domain.  Since the matrix is symmetric,
we only create mesh functions for the upper triangular part of the
matrix. In :py:class:`MeshFunction <dolfin.cpp.mesh.MeshFunction>`
the first argument specifies the type of the mesh function, here we
use "double".  Other types allowed are "int", "size_t" and "bool".
The two following arguments are optional; the first gives the mesh the
:py:class:`MeshFunction <dolfin.cpp.mesh.MeshFunction>` is defined on,
and the second the topological dimension of the mesh function.

.. code-block:: python

    # Create mesh functions for c00, c01, c11
    c00 = MeshFunction("double", mesh, 2)
    c01 = MeshFunction("double", mesh, 2)
    c11 = MeshFunction("double", mesh, 2)

To set the values of the mesh functions, we go through all the cells
in the mesh and check whether the midpoint value of the cell in the
:math:`x`-direction is less than 0.5 or not (in practice this means
that we are checking which half of the unit square the cell is
in). Then we set the correct values of the mesh functions, depending
on which half we are in.

.. code-block:: python

    # Iterate over mesh and set values
    for cell in cells(mesh):
        if cell.midpoint().x() < 0.5:
            c00[cell] = 1.0
            c01[cell] = 0.3
            c11[cell] = 2.0
        else:
            c00[cell] = 3.0
            c01[cell] = 0.5
            c11[cell] = 4.0

Create files to store data and save to file:

.. code-block:: python

    # Store to file
    mesh_file = File("../unitsquare_32_32.xml.gz")
    c00_file = File("../unitsquare_32_32_c00.xml.gz")
    c01_file = File("../unitsquare_32_32_c01.xml.gz")
    c11_file = File("../unitsquare_32_32_c11.xml.gz")

    mesh_file << mesh
    c00_file << c00
    c01_file << c01
    c11_file << c11

Plot the mesh functions using the ``plot`` function:

.. code-block:: python

    # Plot mesh functions
    plot(c00, title="C00")
    plot(c01, title="C01")
    plot(c11, title="C11")

    interactive()

Implementation of tensor-weighted-poisson.py
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

This description goes through the implementation (in
:download:`demo_tensor-weighted-poisson.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 *

We proceed by defining a mesh of the domain and a finite element
function space :math:`V` relative to this mesh. We read the mesh file
generated by :download:`generate_data.py` and create the function
space in the following way:

.. code-block:: python

    # Read mesh from file and create function space
    mesh = Mesh("../unitsquare_32_32.xml.gz")
    V = FunctionSpace(mesh, "Lagrange", 1)

The second argument to :py:class:`FunctionSpace
<dolfin.cpp.function.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).)

.. code-block:: python

    # 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 <dolfin.cpp.fem.DirichletBC>`.  A
:py:class:`DirichletBC <dolfin.cpp.fem.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
:math:`V`. The value of the boundary condition :math:`(0.0)` can be
represented using a :py:class:`Constant
<dolfin.functions.constant.Constant>` and the Dirichlet boundary is
defined immediately above. The definition of the Dirichlet boundary
condition then looks as follows:

.. code-block:: python

    # Define boundary condition
    u0 = Constant(0.0)
    bc = DirichletBC(V, u0, boundary)

Before we define the conductivity matrix, we create a string
containing C++ code for evaluation of the conductivity. Later we will
use this string when we create an :py:class:`Expression
<dolfin.cpp.function.Expression>` containing the entries of the
matrix.

.. code-block:: python

    # Code for C++ evaluation of conductivity
    conductivity_code = """

    class Conductivity : public Expression
    {
    public:

      // Create expression with 3 components
      Conductivity() : Expression(3) {}

      // Function for evaluating expression on each cell
      void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
      {
        const uint D = cell.topological_dimension;
        const uint cell_index = cell.index;
        values[0] = (*c00)[cell_index];
        values[1] = (*c01)[cell_index];
        values[2] = (*c11)[cell_index];
      }

      // The data stored in mesh functions
      std::shared_ptr<MeshFunction<double> > c00;
      std::shared_ptr<MeshFunction<double> > c01;
      std::shared_ptr<MeshFunction<double> > c11;

    };
    """

We define the conductivity matrix by first creating mesh functions
from the files we stored in :download:`generate_data.py`.  Here, the
third argument in :py:class:`MeshFunction
<dolfin.cpp.mesh.MeshFunction>` is the path to the data files.  Then,
we define an expression for the entries in the matrix where we give
the C++ code as an argument for optimalization. Finally, we use the
UFL function ``as_matrix`` to create the matrix consisting of the
expressions.

.. code-block:: python

    # Define conductivity expression and matrix
    c00 = MeshFunction("double", mesh, "../unitsquare_32_32_c00.xml.gz")
    c01 = MeshFunction("double", mesh, "../unitsquare_32_32_c01.xml.gz")
    c11 = MeshFunction("double", mesh, "../unitsquare_32_32_c11.xml.gz")

    c = Expression(cppcode=conductivity_code, degree=0)
    c.c00 = c00
    c.c01 = c01
    c.c11 = c11
    C = as_matrix(((c[0], c[1]), (c[1], c[2])))

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:func:`TrialFunction <dolfin.functions.function.TrialFunction>` and
a :py:func:`TestFunction <dolfin.functions.function.TestFunction>` on
the previously defined :py:class:`FunctionSpace
<dolfin.cpp.function.FunctionSpace>` :math:`V`.

Further, the source :math:`f` is involved in the variational form, and
hence it must be must specified.  Since :math:`f` is given by a simple
mathematical formula, it can easily be declared using the
:py:class:`Expression <dolfin.cpp.function.Expression>` class.  Note
that the string defining :math:`f` uses 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 :math:`a`
and the linear form :math:`L` (using UFL operators).  In summary, this
reads:

.. code-block:: python

    # 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)
    a = inner(C*grad(u), grad(v))*dx
    L = f*v*dx

Now, we have specified the bilinear and linear forms and can consider
the solution of the variational problem.  First, we need to define a
:py:class:`Function <dolfin.cpp.function.Function>` ``u`` to
represent the solution. (Upon initialization, it is simply set to the
zero function.) A :py:class:`Function <dolfin.cpp.function.Function>`
represents a function living in a finite element function space.
Next, we can call the :py:meth:`solve
<dolfin.cpp.fem.GenericAdaptiveVariationalSolver.solve>` function with
the arguments ``a == L``, ``u`` and ``bc`` as follows:

.. code-block:: python

    # 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 <dolfin.cpp.function.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:meth:`plot <dolfin.cpp.io.VTKPlotter.plot>` command:

.. code-block:: python

    # Save solution in VTK format
    file = File("poisson.pvd")
    file << u

    # Plot solution
    plot(u, interactive=True)

Complete code
-------------

demo_tensor-weighted-poisson.py:

.. literalinclude:: demo_tensor-weighted-poisson.py
   :start-after: # Begin demo

generate_data.py:

.. literalinclude:: generate_data.py
   :start-after: # Begin demo