A more general solver function¶

We consider the flat program developed in the section FEniCS implementation. Some of the code in this program is needed to solve any Poisson problem \(-\nabla^2 u=f\) on \([0,1]\times [0,1]\) with \(u=u_{_\mathrm{D}}\) on the boundary, while other statements arise from our simple test problem. Let us collect the general, reusable code in a function called solver. Our special test problem will then just be an application of solver with some additional statements. We limit the solver function to just compute the numerical solution. Plotting and comparing the solution with the exact solution are considered to be problem-specific activities to be performed elsewhere.

We parameterize solver by \(f\), \(u_{_\mathrm{D}}\), and the resolution of the mesh. Since it is so trivial to use higher-order finite element functions by changing the third argument to FunctionSpace, we also add the polynomial degree of the finite element function space as an argument to solver.

from fenics import *

def solver(f, u_D, Nx, Ny, degree=1):
    """
    Solve -Laplace(u) = f on [0,1] x [0,1] with 2*Nx*Ny Lagrange
    elements of specified degree and u=u_D (Expresssion) on
    the boundary.
    """

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)

    return u

The remaining tasks of our initial program, such as calling the solve function with problem-specific parameters and plotting, can be placed in a separate function. Here we choose to put this code in a function named run_solver:

def run_solver():
    "Run solver to compute and post-process solution"

    # Set up problem parameters and call solver
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')
    f = Constant(-6.0)
    u = solver(f, u_D, 6, 4, 1)

    # Plot solution
    u.rename('u', 'u')
    plot(u)

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

The solution can now be computed, plotted, and saved to file by simply calling the run_solver function.

Writing the solver as a Python module¶

The refactored code is put in a file ft13_poisson_solver.py. We should make sure that such a file can be imported (and hence reused) in other programs. This means that all statements in the main program that are not inside functions should appear within a test if __name__ == '__main__':. This test is true if the file is executed as a program, but false if the file is imported. If we want to run this file in the same way as we can run ft13_poisson_solver.py, the main program is simply a call to run_solver() followed by a call interactive() to hold the plot:

if __name__ == '__main__':
    run_solver()
    interactive()

The complete code can be found in the file ft13_poisson_solver.py.

Verification and unit tests¶

The remaining part of our first program is to compare the numerical and the exact solutions. Every time we edit the code we must rerun the test and examine that max_error is sufficiently small so we know that the code still works. To this end, we shall adopt unit testing, meaning that we create a mathematical test and corresponding software that can run all our tests automatically and check that all tests pass. Python has several tools for unit testing. Two very popular ones are pytest and nose. These are almost identical and very easy to use. More classical unit testing with test classes is offered by the built-in module unittest, but here we are going to use pytest (or nose) since that will result in shorter and clearer code.

Mathematically, our unit test is that the finite element solution of our problem when \(f=-6\) equals the exact solution \(u=u_{_\mathrm{D}}=1+x^2+2y^2\) at the vertices of the mesh. We have already created a code that finds the error at the vertices for our numerical solution. Because of rounding errors, we cannot demand this error to be zero, but we have to use a tolerance, which depends to the number of elements and the degrees of the polynomials in the finite element basis functions. If we want to test that the solver function works for meshes up to \(2\times(20\times 20)\) elements and cubic Lagrange elements, \(10^{-10}\) is an appropriate tolerance for testing that the maximum error vanishes (see the section Dissection of the program).

To make our test case work together with pytest and nose, we have to make a couple of small adjustments to our program. The simple rule is that each test must be placed in a function that

• has a name starting with test_,
• has no arguments, and
• implements a test expressed as assert success, msg.

Regarding the last point, success is a boolean expression that is False if the test fails, and in that case the string msg is written to the screen. When the test fails, assert raises an AssertionError exception in Python, and otherwise runs silently. The msg string is optional, so assert success is the minimal test. In our case, we will write assert max_error < tol, where tol is the tolerance mentioned above.

A proper test function for implementing this unit test in the pytest or nose testing frameworks has the following form. Note that we perform the test for different mesh resolutions and degrees of finite elements.

def test_solver():
    "Test solver by reproducing u = 1 + x^2 + 2y^2"

    # Set up parameters for testing
    tol = 1E-10
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')
    f = Constant(-6.0)

    # Iterate over mesh sizes and degrees
    for Nx, Ny in [(3, 3), (3, 5), (5, 3), (20, 20)]:
        for degree in 1, 2, 3:
            print('Solving on a 2 x (%d x %d) mesh with P%d elements.'
                  % (Nx, Ny, degree))

            # Compute solution
            u = solver(f, u_D, Nx, Ny, degree)

            # Extract the mesh
            mesh = u.function_space().mesh()

            # Compute maximum error at vertices
            vertex_values_u_D = u_D.compute_vertex_values(mesh)
            vertex_values_u  = u.compute_vertex_values(mesh)
            import numpy as np
            error_max = np.max(np.abs(vertex_values_u_D - \
                                      vertex_values_u))

            # Check maximum error
            msg = 'error_max = %g' % error_max
            assert error_max < tol, msg

if __name__ == '__main__':
    run_solver()
    interactive()

To run the test, we type the following command:

Terminal> py.test ft13_poisson_solver.py

This will run all functions test_*() (currently only the test_solver function) found in the file and report the results. For more verbose output, add the flags -s -v.

We shall make it a habit in the following test to encapsulate numerical test problems in unit tests as done above, and we strongly encourage the reader to create similar unit tests whenever a FEniCS solver is implemented.

Parameterizing the number of space dimensions¶

FEniCS makes it is easy to write a unified simulation code that can operate in 1D, 2D, and 3D. We will conveniently make use of this feature in forthcoming examples. As an appetizer, go back to the introductory programs ft01_poisson.py or ft13_poisson_solver.py and change the mesh construction from UnitSquareMesh(16, 16) to UnitCubeMesh(16, 16, 16). Now the domain is the unit cube partitioned into \(16\times 16\times16\) boxes, and each box is divided into six tetrahedron-shaped finite elements for computations. Run the program and observe that we can solve a 3D problem without any other modifications (!). The visualization allows you to rotate the cube and observe the function values as colors on the boundary.

If we want to parameterize the creation of unit interval, unit square, or unit cube over dimension, we can do so by encapsulating this part of the code in a function. Given a list or tuple with the divisions into cells in the various spatial coordinates, the following function returns the mesh for a \(d\)-dimensional cube:

def UnitHyperCube(divisions):
    mesh_classes = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
    d = len(divisions)
    mesh = mesh_classes[d-1](*divisions)
    return mesh

The construction mesh_class[d-1] will pick the right name of the object used to define the domain and generate the mesh. Moreover, the argument *divisions sends all the components of the list divisions as separate arguments to the constructor of the mesh construction class picked out by mesh_class[d-1]. For example, in a 2D problem where divisions has two elements, the statement

mesh = mesh_classes[d-1](*divisions)

is equivalent to

mesh = UnitSquareMesh(divisions[0], divisions[1])

The solver function from ft13_poisson_solver.py may be modified to solve \(d\)-dimensional problems by replacing the Nx and Ny parameters by divisions, and calling the function UnitHyperCube to create the mesh. Note that UnitHyperCube is a function and not a class, but we have named it with CamelCase to make it look like a class:

mesh = UnitHyperCube(divisions)

Exercise 4: Solve a Poisson problem¶

Solve the following problem

The exact solution is given by

Compute the maximum numerical approximation error in a mesh with \(2(N_x\times N_y)\) elements and in a mesh with double resolution: \(4(N_x\times N_y)\) elements. Show that the doubling the resolution reduces the error by a factor 4 when using Lagrange elements of degree one. Make an illustrative plot of the solution too.

a) Base your implementation on editing the program ft01_poisson.py.

Hint 1. In the string for an Expression object, pi is the value of \(\pi\). Also note that \(\pi^2\) must be expressed with syntax pow(pi,2) and not (the common Python syntax) pi**2.

FEniCS will abort with a compilation error if you type the expressions in a wrong way syntax-wise. Search for error: in the /very/long/path/compile.log file mentioned in the error message to see what the C++ compiler reported as error in the expressions.

Hint 2. The result that with P1 elements, doubling the resolution reduces the error with a factor of four, is an asymptotic result so it requires a sufficiently fine mesh. Here one may start with \(N_x=N_y=20\).

Filename: poisson_fsin_flat.

Solution. Looking at the ft01_poisson.py code, we realize that the following edits are required:

• Modify the mesh computation.
• Modify u_b and f.
• Add expression for the exact solution.
• Modify the computation of the numerical error.
• Insert a loop to enable solving the problem twice.
• Put the error reduction computation and the plot statements after the loop.

Here is the modified code:

from fenics import *

Nx = Ny = 20
error = []
for i in range(2):
    Nx *= (i+1)
    Ny *= (i+1)

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'Lagrange', 1)

    # Define boundary conditions
    u0 = Constant(0)

    def u0_boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u0, u0_boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Expression('-2*exp(-2*x[0])*sin(pi*x[1])*('
                   '(4-5*pow(pi,2))*sin(2*pi*x[0]) '
                   ' - 8*pi*cos(2*pi*x[0]))')
    # Note: no need for pi=DOLFIN_PI in f, pi is valid variable
    a = inner(nabla_grad(u), nabla_grad(v))*dx
    L = f*v*dx

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)

    u_e = Expression(
        '2*exp(-2*x[0])*sin(2*pi*x[0])*sin(pi*x[1])')

    u_e_Function = interpolate(u_e, V)         # exact solution
    u_e_array = u_e_Function.vector().array()  # dof values
    max_error = (u_e_array - u.vector().array()).max()
    print('max error:', max_error, '%dx%d mesh' % (Nx, Ny))
    error.append(max_error)

print('Error reduction:', error[1]/error[0])

# Plot solution and mesh
plot(u)

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

# Hold plot
interactive()

The number \(\pi\) has the symbol M_PI in C and C++, but in C++ strings in Expression objects, the symbol pi can be used directly (or one can use the less readable DOLFIN_PI).

b) Base your implementation on a new file that imports functionality from the module ft13_poisson_solver.py. Embed the check of the reduction of the numerical approximation error in a unit test. Filename: poisson_fsin_func.

Solution. Solving the two problems is a matter of calling solver with different sets of arguments. To compute the numerical error, we need code that is close to what we have in test_solver.

from poisson_solver import (
    solver, Expression, Constant, interpolate, File, plot,
    interactive)

def data():
    """Return data for this Poisson problem."""
    u0 = Constant(0)
    u_e = Expression(
        '2*exp(-2*x[0])*sin(2*pi*x[0])*sin(pi*x[1])')
    f = Expression('-2*exp(-2*x[0])*sin(pi*x[1])*('
                   '(4-5*pow(pi,2))*sin(2*pi*x[0]) '
                   ' - 8*pi*cos(2*pi*x[0]))')
    return u0, f, u_e

def test_solver():
    """Check convergence rate of solver."""
    u0, f, u_e = data()
    Nx = 20
    Ny = Nx
    error = []
    # Loop over refined meshes
    for i in range(2):
        Nx *= i+1
        Ny *= i+1
        print('solving on 2(%dx%d) mesh' % (Nx, Ny))
        u = solver(f, u0, Nx, Ny, degree=1)
        # Make a finite element function of the exact u_e
        V = u.function_space()
        u_e_array = interpolate(u_e, V).vector().array()
        max_error = (u_e_array - u.vector().array()).max()  # Linf norm
        error.append(max_error)
        print('max error:', max_error)
    for i in range(1, len(error)):
        error_reduction = error[i]/error[i-1]
        print('error reduction:', error_reduction)
        assert abs(error_reduction - 0.25) < 0.1

def application():
    """Plot the solution."""
    u0, f, u_e = data()
    Nx = 40
    Ny = Nx
    u = solver(f, u0, Nx, Ny, 1)
    # Dump solution to file in VTK format
    file = File("poisson.pvd")
    file << u
    # Plot solution and mesh
    plot(u)

if __name__ == '__main__':
    test_solver()
    application()
    # Hold plot
    interactive()

The unit test is embedded in a proper test function test_solver for the pytest or nose testing frameworks. Visualization of the solution is encapsulated in the application function. Since we need u_e, u_b, and f in two functions, we place the definitions in a function data to avoid copies of these expressions.

Exercise 5: Refactor the code for membrane deflection¶

The ft14_membrane.py program simulates the deflection of a membrane. Refactor this code such that we have a solver function as in the program with name ft13_poisson_solver.py. Let the user have the option to choose a direct or iterative solver for the linear system. Also implement a unit test where you have \(p=4\) (constant) and use P2 and P3 elements. In this case, the exact solution is quadratic in \(x\) and \(y\) and will be “exactly” reproduced by P2 and higher-order elements.

Solution. We can use the solver function from ft13_poisson_solver.py right away. The major difference is that the domain is now a circle and not a square. We change the solver function by letting the mesh be an argument mesh (instead of Nx and Ny):

def solver(
    f, u_b, mesh, degree=1,
    linear_solver='Krylov', # Alt: 'direct'
    ...):
    V = FunctionSpace(mesh, 'P', degree)
    # code as before

The complete code becomes

def application(beta, R0, num_elements_radial_dir):
    # Scaled pressure function
    p = Expression(
        '4*exp(-pow(beta,2)*(pow(x[0], 2) + pow(x[1]-R0, 2)))',
        beta=beta, R0=R0)

    # Generate mesh over the unit circle
    domain = Circle(Point(0.0, 0.0), 1.0)
    mesh = generate_mesh(domain, num_elements_radial_dir)

    w = solver(p, Constant(0), mesh, degree=1,
               linear_solver='direct')
    w.rename('w', 'deflection')  # set name and label (description)

    # Plot scaled solution, mesh and pressure
    plot(mesh, title='Mesh over scaled domain')
    plot(w, title='Scaled ' + w.label())
    V = w.function_space()
    p = interpolate(p, V)
    p.rename('p', 'pressure')
    plot(p, title='Scaled ' + p.label())

    # Dump p and w to file in VTK format
    vtkfile1 = File('membrane_deflection.pvd')
    vtkfile1 << w
    vtkfile2 = File('membrane_load.pvd')
    vtkfile2 << p

The key function to simulate membrane deflection is named application.

For \(p=4\), we have \(w=1-x^2-y^2\) as exact solution. The unit test for P2 and P3 goes as follows:

def test_membrane():
    """Verification for constant pressure."""
    p = Constant(4)
    # Generate mesh over the unit circle
    domain = Circle(Point(0.0, 0.0), 1.0)
    for degree in 2, 3:
        print('********* P%d elements:' % degree)
        n = 5
        for i in range(4):  # Run some resolutions
            n *= (i+1)
            mesh = generate_mesh(domain, n)
            #info(mesh)
            w = solver(p, Constant(0), mesh, degree=degree,
                       linear_solver='direct')
            print('max w: %g, w(0,0)=%g, h=%.3E, dofs=%d' %
                  (w.vector().array().max(), w((0,0)),
                   1/np.sqrt(mesh.num_vertices()),
                   w.function_space().dim()))

            w_exact = Expression('1 - x[0]*x[0] - x[1]*x[1]')
            w_e = interpolate(w_exact, w.function_space())
            error = np.abs(w_e.vector().array() -
                           w.vector().array()).max()
            print('error: %.3E' % error)
            assert error < 9.61E-03

def application2(
    beta, R0, num_elements_radial_dir):
    """Explore more built-in visulization features."""
    # Scaled pressure function
    p = Expression(
        '4*exp(-pow(beta,2)*(pow(x[0], 2) + pow(x[1]-R0, 2)))',
        beta=beta, R0=R0)

    # Generate mesh over the unit circle
    domain = Circle(Point(0.0, 0.0), 1.0)
    mesh = generate_mesh(domain, num_elements_radial_dir)

    w = solver(p, Constant(0), mesh, degree=1,
               linear_solver='direct')
    w.rename('w', 'deflection')

    # Plot scaled solution, mesh and pressure
    plot(mesh, title='Mesh over scaled domain')
    viz_w = plot(w,
                 wireframe=False,
                 title='Scaled membrane deflection',
                 axes=False,
                 interactive=False,
                 )
    viz_w.elevate(-10) # adjust (lift) camera from default view
    viz_w.plot(w)      # bring new settings into action
    viz_w.write_png('deflection')
    viz_w.write_pdf('deflection')

    V = w.function_space()
    p = interpolate(p, V)
    p.rename('p', 'pressure')
    viz_p = plot(p, title='Scaled pressure', interactive=False)
    viz_p.elevate(-10)
    viz_p.plot(p)
    viz_p.write_png('pressure')
    viz_p.write_pdf('pressure')

    # Dump w and p to file in VTK format
    vtkfile1 = File('membrane_deflection.pvd')
    vtkfile1 << w
    vtkfile2 = File('membrane_load.pvd')
    vtkfile2 << p

The striking feature is that the solver does not reproduce the solution to an accuracy more than about 0.01 (!), regardless of the resolution and type of element.

Filename: membrane_func.

Controlling the solution process¶

solve(a == L, u, bc)
      solver_parameters={'linear_solver': 'gmres',
                         'preconditioner': 'ilu'})
# Alternative syntax
solve(a == L, u, bc,
      solver_parameters=dict(linear_solver='gmres',
                             preconditioner='ilu'))

parameters['linear_algebra_backend'] = backendname

print parameters['linear_algebra_backend']

# Alternative syntax for Python 3
print(parameters['linear_algebra_backend'])

info(parameters, verbose=True)

prm = parameters['krylov_solver'] # short form
prm['absolute_tolerance'] = 1E-10
prm['relative_tolerance'] = 1E-6
prm['maximum_iterations'] = 1000

File('fenics_parameters.xml') << parameters

File('fenics_parameters.xml') >> parameters

def solver(
    f, u_b, Nx, Ny,
    degree=1,               # Polynomial degree of function space
    linear_solver='Krylov', # Linear solver method (alt: 'direct')
    abs_tol=1E-5,           # Absolute tolerance in Krylov solver
    rel_tol=1E-3,           # Relative tolerance in Krylov solver
    max_iter=1000,          # Max iterations in Krylov solver
    log_level=PROGRESS,     # Amount of solver output
    print_parameters=False, # Print solver parameters to screen?
    ):
    ...
    # Set up variational problem: a, L, declare u, etc.

    if linear_solver == 'Krylov':
        prm = parameters['krylov_solver']
        prm['absolute_tolerance'] = abs_tol
        prm['relative_tolerance'] = rel_tol
        prm['maximum_iterations'] = max_iter
        print(parameters['linear_algebra_backend'])
        set_log_level(log_level)
        if print_parameters:
            info(parameters, True)
        solver_parameters = {'linear_solver': 'gmres',
                             'preconditioner': 'ilu'}
    else:
        solver_parameters = {'linear_solver': 'lu'}

    solve(a == L, u, bc, solver_parameters=solver_parameters)
    return u

List of linear solver methods and preconditioners¶

Which linear solvers and preconditioners that are available in FEniCS depends on how FEniCS has been configured and which linear algebra backend is currently active. The following table shows an example of which linear solvers that can be available through FEniCS when the PETSc backend is active:

The set of available preconditioners also depends on configuration and linear algebra backend. The following table shows an example of which preconditioners may be available:

An up-to-date list of the available solvers and preconditioners for your FEniCS installation can be produced by

list_linear_solver_methods()
list_krylov_solver_preconditioners()

Linear variational problem and solver objects¶

The FEniCS interface allows different ways to access the core functionality, ranging from very high-level to low-level access. So far, we have mostly used the high-level call solve(a == L, u, bc) to solve a variational problem a == L with a certain boundary condition bc. However, sometimes you may need more fine-grained control over the solution process. In particular, the call to solve will create certain objects that are thrown away after the solution has been computed, and it may be practical or efficient to reuse those objects.

In this section, we will look at an alternative interface to solving linear variational problems in FEniCS, which may be preferable in many situations compared to the high-level solve function interface. This interface uses the two classes LinearVariationalProblem and LinearVariationalSolver. Using this interface, the equivalent of solve(a == L, u, bc) looks as follows:

u = Function(V)
problem = LinearVariationalProblem(a, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.solve()

Many FEniCS objects have an attribute parameters corresponding to a parameter set in the global parameters database, but local to the object. Here, solver.parameters play that role. Setting the CG method with ILU preconditioning as the solution method and specifying solver-specific parameters can be done like this:

solver.parameters['linear_solver'] = 'gmres'
solver.parameters['preconditioner'] = 'ilu'
prm = solver.parameters['krylov_solver'] # short form
prm['absolute_tolerance'] = 1E-7
prm['relative_tolerance'] = 1E-4
prm['maximum_iterations'] = 1000

Settings in the global parameters database are propagated to parameter sets in individual objects, with the possibility of being overwritten as above. Note that global parameter values can only affect local parameter values if set before the time of creation of the local object. Thus, changing the value of the tolerance in the global parameter database will not affect the parameters for already created solvers.

The linear variational problem and solver objects as outlined above are incorporated in an alternative solver function, named solver_objects, in ft15_poisson_iterative.py. Otherwise, this function is similar to the previously shown solver function.

Explicit assembly and solve¶

As we saw already in the section The Navier - Stokes equations, linear variational problems can be assembled explicitly in FEniCS into matrices and vectors using the assemble function. This allows even more fine-grained control of the solution process compared to using the high-level solve function or using the classes LinearVariationalProblem and LinearVariationalSolver. We will now look more closely into how to use the assemble function and how to combine this with low-level calls for solving the assembled linear systems.

Given a variational problem \(a(u,v)=L(v)\), the discrete solution \(u\) is computed by inserting \(u=\sum_{j=1}^N U_j \phi_j\) into \(a(u,v)\) and demanding \(a(u,v)=L(v)\) to be fulfilled for \(N\) test functions \(\hat\phi_1,\ldots,\hat\phi_N\). This implies

which is nothing but a linear system,

where the entries of \(A\) and \(b\) are given by

The examples so far have specified the left- and right-hand sides of the variational formulation and then asked FEniCS to assemble the linear system and solve it. An alternative is to explicitly call functions for assembling the coefficient matrix \(A\) and the right-hand side vector \(b\), and then solve the linear system \(AU=b\) for the vector \(U\). Instead of solve(a == L, U, b) we now write

A = assemble(a)
b = assemble(L)
bc.apply(A, b)
u = Function(V)
U = u.vector()
solve(A, U, b)

The variables a and L are the same as before; that is, a refers to the bilinear form involving a TrialFunction object u and a TestFunction object v, and L involves the same TestFunction object v. From a and L, the assemble function can compute \(A\) and \(b\).

Creating the linear system explicitly in a program can have some advantages in more advanced problem settings. For example, \(A\) may be constant throughout a time-dependent simulation, so we can avoid recalculating \(A\) at every time level and save a significant amount of simulation time.

The matrix \(A\) and vector \(b\) are first assembled without incorporating essential (Dirichlet) boundary conditions. Thereafter, the call bc.apply(A, b) performs the necessary modifications of the linear system such that u is guaranteed to equal the prescribed boundary values. When we have multiple Dirichlet conditions stored in a list bcs, we must apply each condition in bcs to the system:

for bc in bcs:
    bc.apply(A, b)

# Alternative syntax using list comprehension
[bc.apply(A, b) for bc in bcs]

Alternatively, we can use the function assemble_system, which takes the boundary conditions into account during the assembly of the matrix and vector:

A, b = assemble_system(a, L, bcs)

The assemble_system function is preferable to the combination of assemble and bc.apply when the linear system is symmetric, since assemble_system will incorporate the boundary conditions in a symmetric way. Even if the matrix A that comes out of the call to assemble is symmetric (for a symmetric bilinear form a), the call to bc.apply will break the symmetry.

Once the linear system has been assembled, we need to compute the solution \(U=A^{-1}b\) and store the solution \(U\) in the vector U = u.vector(). In the same way as linear variational problems can be programmed using different interfaces in FEniCS—the high-level solve function, the class LinearVariationalSolve, and the low-level assemble function—linear systems can also be programmed using different interfaces in FEniCS. The high-level interface to solving a linear system in FEniCS is also named solve:

solve(A, U, b)

By default, solve(A, U, b) uses sparse LU decomposition to compute the solution. Specification of an iterative solver and preconditioner can be made through two optional arguments:

solve(A, U, b, 'cg', 'ilu')

Appropriate names of solvers and preconditioners are found in the section List of linear solver methods and preconditioners.

This high-level interface is useful for many applications, but sometimes more fine-grained control is needed. One can then create one or more KrylovSolver objects that are then used to solve linear systems. Each different solver object can have its own set of parameters and selection of iterative method and preconditioner. Here is an example:

solver = KrylovSolver('cg', 'ilu')
prm = solver.parameters
prm['absolute_tolerance'] = 1E-7
prm['relative_tolerance'] = 1E-4
prm['maximum_iterations'] = 1000
u = Function(V)
U = u.vector()
solver.solve(A, U, b)

The function solver_linalg in the program file ft15_poisson_iterative.py implements a solver function where the user can choose between different types of assembly. The function demo_linalg runs a test problem on a sequence of meshes and solves the problem with symmetric and non-symmetric modification of the system matrix. One can monitor the number of Krylov method iterations and realize that with a symmetric coefficient matrix, the Conjugate Gradient method requires slightly fewer iterations than GMRES in the non-symmetric case. Taking into account that the Conjugate Gradient method has less work per iteration, there is some efficiency to be gained by using assemble_system for this problem.

The choice of start vector for the iterations in a linear solver is often important. By default, the values of u and thus the vector U = u.vector() will be initialized to zero. If we instead wanted to initialize U with random numbers in the interval \([-100,100]\) this can be done as follows:

n = u.vector().array().size
U = u.vector()
U[:] = numpy.random.uniform(-100, 100, n)
solver.parameters['nonzero_initial_guess'] = True
solver.solve(A, U, b)

Note that we must both turn off the default behavior of setting the start vector (“initial guess”) to zero, and also set the values of the vector U to nonzero values.

Using a nonzero initial guess can be particularly important for time-dependent problems or when solving a linear system as part of a nonlinear iteration, since then the previous solution vector U will often be a good initial guess for the solution in the next time step or iteration. In this case, the values in the vector U will naturally be initialized with the previous solution vector (if we just used it to solve a linear system), so the only extra step necessary is to set the parameter nonzero_initial_guess to True.

Examining matrix and vector values¶

When calling A = assemble(a) and b = assemble(L), the object A will be of type Matrix, while b and u.vector() are of type Vector. To examine the values, we may convert the matrix and vector data to numpy arrays by calling the array() method as shown before. For example, if you wonder how essential boundary conditions are incorporated into linear systems, you can print out A and b before and after the bc.apply(A, b) call:

A = assemble(a)
b = assemble(L)
if mesh.num_cells() < 16: # print for small meshes only
    print(A.array())
    print(b.array())
bc.apply(A, b)
if mesh.num_cells() < 16:
    print(A.array())
    print(b.array())

With access to the elements in A through a numpy array, we can easily perform computations on this matrix, such as computing the eigenvalues (using the eig function in numpy.linalg). We can alternatively dump A.array() and b.array() to file in MATLAB format and invoke MATLAB or Octave to analyze the linear system. Dumping the arrays to MATLAB format is done by

import scipy.io
scipy.io.savemat('Ab.mat', {'A': A.array(), 'b': b.array()})

Writing load Ab.mat in MATLAB or Octave will then make the array variables A and b available for computations.

Matrix processing in Python or MATLAB/Octave is only feasible for small PDE problems since the numpy arrays or matrices in MATLAB file format are dense matrices. FEniCS also has an interface to the eigensolver package SLEPc, which is the preferred tool for computing the eigenvalues of large, sparse matrices of the type encountered in PDE problems (see demo/documented/eigenvalue in the FEniCS source code tree for a demo).

Examining the degrees of freedom¶

We have seen before how to grab the degrees of freedom array from a finite element function u:

nodal_values = u.vector().array()

For a finite element function from a standard continuous piecewise linear function space (\(\mathsf{P}_1\) Lagrange elements), these values will be the same as the values we get by the following statement:

vertex_values = u.compute_vertex_values(mesh)

Both nodal_values and vertex_values will be numpy arrays and they will be of the same length and contain the same values, but with possibly different ordering. The array vertex_values will have the same ordering as the vertices of the mesh, while nodal_values will be ordered in a way that (nearly) minimizes the bandwidth of the system matrix and thus improves the efficiency of linear solvers.

A fundamental question is: what are the coordinates of the vertex whose value is nodal_values[i]? To answer this question, we need to understand how to get our hands on the coordinates, and in particular, the numbering of degrees of freedom and the numbering of vertices in the mesh.

The function mesh.coordinates() returns the coordinates of the vertices as a numpy array with shape \((M,d\)), \(M\) being the number of vertices in the mesh and \(d\) being the number of space dimensions:

>>> from fenics import *
>>>
>>> mesh = UnitSquareMesh(2, 2)
>>> coordinates = mesh.coordinates()
>>> coordinates
array([[ 0. ,  0. ],
       [ 0.5,  0. ],
       [ 1. ,  0. ],
       [ 0. ,  0.5],
       [ 0.5,  0.5],
       [ 1. ,  0.5],
       [ 0. ,  1. ],
       [ 0.5,  1. ],
       [ 1. ,  1. ]])

We see from this output that for this particular mesh, the vertices are first numbered along \(y=0\) with increasing \(x\) coordinate, then along \(y=0.5\), and so on.

Next we compute a function u on this mesh. Let’s take \(u=x+y\):

>>> V = FunctionSpace(mesh, 'P', 1)
>>> u = interpolate(Expression('x[0] + x[1]'), V)
>>> plot(u, interactive=True)
>>> nodal_values = u.vector().array()
>>> nodal_values
array([ 1. ,  0.5,  1.5,  0. ,  1. ,  2. ,  0.5,  1.5,  1. ])

We observe that nodal_values[0] is not the value of \(x+y\) at vertex number 0, since this vertex has coordinates \(x=y=0\). The numbering of the nodal values (degrees of freedom) \(U_1,\ldots,U_{N}\) is obviously not the same as the numbering of the vertices.

Note that we may examine the vertex numbering in the plot window. We type w to turn on wireframe instead of a fully colored surface, m to show the mesh, and then v to show the numbering of the vertices.

Let’s instead examine the values we get by calling u.compute_vertex_values():

>>> vertex_values = u.compute_vertex_values()
>>> for i, x in enumerate(coordinates):
...     print('vertex %d: vertex_values[%d] = %g\tu(%s) = %g' %
...           (i, i, vertex_values[i], x, u(x)))
vertex 0: vertex_values[0] = 0      u([ 0.  0.]) = 8.46545e-16
vertex 1: vertex_values[1] = 0.5    u([ 0.5  0. ]) = 0.5
vertex 2: vertex_values[2] = 1          u([ 1.  0.]) = 1
vertex 3: vertex_values[3] = 0.5    u([ 0.   0.5]) = 0.5
vertex 4: vertex_values[4] = 1          u([ 0.5  0.5]) = 1
vertex 5: vertex_values[5] = 1.5    u([ 1.   0.5]) = 1.5
vertex 6: vertex_values[6] = 1          u([ 0.  1.]) = 1
vertex 7: vertex_values[7] = 1.5    u([ 0.5  1. ]) = 1.5
vertex 8: vertex_values[8] = 2          u([ 1.  1.]) = 2

We can ask FEniCS to give us the mapping from vertices to degrees of freedom for a certain function space \(V\):

v2d = vertex_to_dof_map(V)

Now, nodal_values[v2d[i]] will give us the value of the degree of freedom in u corresponding to vertex i (v2d[i]). In particular, nodal_values[v2d] is an array with all the elements in the same (vertex numbered) order as coordinates. The inverse map, from degrees of freedom number to vertex number is given by dof_to_vertex_map(V). This means that we may call coordinates[dof_to_vertex_map(V)] to get an array of all the coordinates in the same order as the degrees of freedom. Note that these mappings are only available in FEniCS for \(\mathsf{P}_1\) elements.

For Lagrange elements of degree larger than 1, there are degrees of freedom (nodes) that do not correspond to vertices. For these elements, we may get the vertex values by calling u.compute_vertex_values(mesh), and we can get the degrees of freedom by calling u.vector().array(). To get the coordinates associated with all degrees of freedom, we need to iterate over the elements of the mesh and ask FEniCS to return the coordinates and dofs associated with each cell. This information is stored in the FiniteElement and DofMap object of a FunctionSpace. The following code illustrates how to iterate over all elements of a mesh and print the degrees of freedom and coordinates associated with the element.

element = V.element()
dofmap = V.dofmap()
for cell in cells(mesh):
    print(element.tabulate_dof_coordinates(cell))
    print(dofmap.cell_dofs(cell.index()))

To demonstrate the use of point evaluation of Function objects, we write out the computed u at the center point of the domain and compare it with the exact solution:

center = (0.5, 0.5)
error = u_D(center) - u(center)
print('Error at %s: %g' % (center, error))

For a \(2\times(3\times 3)\) mesh, the output from the previous snippet becomes

Error at (0.5, 0.5): -0.0833333

The discrepancy is due to the fact that the center point is not a node in this particular mesh, but a point in the interior of a cell, and u varies linearly over the cell while u_D is a quadratic function. When the center point is a node, as in a \(2\times(2\times 2)\) or \(2\times(4\times 4)\) mesh, the error is of the order \(10^{-15}\).

We have seen how to extract the nodal values in a numpy array. If desired, we can adjust the nodal values too. Say we want to normalize the solution such that \(\max_j |U_j| = 1\). Then we must divide all \(U_j\) values by \(\max_j |U_j|\). The following function performs the task:

def normalize_solution(u):
    """Normalize solution by dividing by max(|u|)."""
    nodal_values = u.vector().array()
    u_max = np.abs(nodal_values).max()
    nodal_values /= u_max
    u.vector()[:] = nodal_values
    #u.vector().set_local(dofs) # alternative

The /= operator implies an in-place modification of the object on the left-hand side: all elements of the array nodal_values are divided by the value u_max. Alternatively, we could do nodal_values = nodal_values / u_max, which implies creating a new array on the right-hand side and assigning this array to the name nodal_values.

All the code in this can be found in the file ft15_poisson_iterative.py (Poission solver with use of iterative methods).

A variable-coefficient Poisson problem¶

As a test problem, we will extend the Poisson problem from the chapter Fundamentals: Solving the Poisson equation with a variable coefficient \(\kappa(x,y)\) in the Laplace operator:

Let us continue to use our favorite solution \(u(x,y)=1+x^2+2y^2\) and then prescribe \(\kappa(x,y)=x+y\). It follows that \(u_{_\mathrm{D}}(x,y) = 1 + x^2 + 2y^2\) and \(f(x,y)=-8x-10y\).

We shall quickly demonstrate that this simple extension of our model problem only requires an equally simple extension of the FEniCS program. The following simple changes must be made to the previously shown codes:

• the solver function must take k (\(\kappa\)) as an argument,
• a new Expression k must be defined for the variable coefficient,
• the right-hand side f must be an Expression since it is no longer a constant,
• the formula for \(a(u,v)\) in the variational problem must be updated.

We first address the modified variational problem. Multiplying the PDE by a test function \(v\) and integrating by parts now results in

The function spaces for \(u\) and \(v\) are the same as in the problem with \(\kappa=1\), implying that the boundary integral vanishes since \(v=0\) on \(\partial\Omega\) where we have Dirichlet conditions. The variational forms \(a\) and \(L\) in the variational problem \(a(u,v)=L(v)\) then become

We must thus replace

a = dot(grad(u), grad(v))*dx

by

a = k*dot(grad(u), grad(v))*dx

Moreover, the definitions of k and f in the test problem read

k = Expression('x[0] + x[1]')
f = Expression('-8*x[0] - 10*x[1]')

No additional modifications are necessary. The file ft11_poisson_bcs.py (Poisson problem, variable coefficients) is a copy of ft15_poisson_iterative.py with the mentioned changes incorporated. Observe that \(\kappa=1\) recovers the original problem in ft15_poisson_iterative.py.

You can execute the file and confirm that it recovers the exact solution \(u\) at the nodes.

Flux computations¶

It is often of interest to compute the flux \(Q = -\kappa\nabla u\). Since \(u = \sum_{j=1}^N U_j \phi_j\), it follows that

However, the gradient of a piecewise continuous finite element scalar field is a discontinuous vector field since the basis functions \(\{\phi_j\}\) have discontinuous derivatives at the boundaries of the cells. For example, using Lagrange elements of degree 1, \(u\) is linear over each cell, and the gradient becomes a piecewise constant vector field. On the contrary, the exact gradient is continuous. For visualization and data analysis purposes, we often want the computed gradient to be a continuous vector field. Typically, we want each component of \(\nabla u\) to be represented in the same way as \(u\) itself. To this end, we can project the components of \(\nabla u\) onto the same function space as we used for \(u\). This means that we solve \(w = \nabla u\) approximately by a finite element method, using the same elements for the components of \(w\) as we used for \(u\). This process is known as projection.

Projection is a common operation in finite element analysis and FEniCS has a function for easily performing the projection: project(expression, W), which returns the projection of some expression into the space W.

In our case, the flux \(Q = -\kappa\nabla u\) is vector-valued and we need to pick W as the vector-valued function space of the same degree as the space V where u resides:

V = u.function_space()
mesh = V.mesh()
degree = V.ufl_element().degree()
W = VectorFunctionSpace(mesh(), 'P', degree)

grad_u = project(grad(u), W)
flux_u = project(-k*grad(u), W)

The applications of projection are many, including turning discontinuous gradient fields into continuous ones, comparing higher- and lower-order function approximations, and transforming a higher-order finite element solution down to a piecewise linear field, which is required by many visualization packages.

Plotting the flux vector field is naturally as easy as plotting anything else:

plot(flux_u, title='flux field')

flux_x, flux_y = flux_u.split(deepcopy=True) # extract components
plot(flux_x, title='x-component of flux (-kappa*grad(u))')
plot(flux_y, title='y-component of flux (-kappa*grad(u))')

The deepcopy=True argument signifies a deep copy, which is a general term in computer science implying that a copy of the data is returned. (The opposite, deepcopy=False, means a shallow copy, where the returned objects are just pointers to the original data.)

For data analysis of the nodal values of the flux field, we can grab the underlying numpy arrays (which demands a deepcopy=True in the split of flux):

flux_x_nodal_values = flux_x.vector().dofs()
flux_y_nodal_values = flux_y.vector().dofs()

The degrees of freedom of the flux_u vector field can also be reached by

flux_u_nodal_values = flux_u.vector().array()

However, this is a flat numpy array containing the degrees of freedom for both the \(x\) and \(y\) components of the flux and the ordering of the components may be mixed up by FEniCS in order to improve computational efficiency.

The function demo_test_flux in the program ft11_poisson_bcs.py demonstrates the computations described above.

Computing functionals¶

After the solution \(u\) of a PDE is computed, we occasionally want to compute functionals of \(u\), for example,

which often reflects some energy quantity. Another frequently occurring functional is the error

where \({u_{\small\mbox{e}}}\) is the exact solution. The error is of particular interest when studying convergence properties of finite element methods. Other times, we may instead be interested in computing the flux out through a part \(\Gamma\) of the boundary \(\partial\Omega\),

where \(n\) is an outward unit normal at \(\Gamma\) and \(\kappa\) is a coefficient (see the problem in the section Postprocessing computations for a specific example).

All these functionals are easy to compute with FEniCS, as we shall see in the examples below.

energy = 0.5*dot(grad(u), grad(u))*dx
E = assemble(energy)

error = (u_e - u)**2*dx
E = sqrt(abs(assemble(error)))

E = errornorm(u_e, u)

n = FacetNormal(mesh)
flux = -k*dot(grad(u), n)*ds
total_flux = assemble(flux)

Computing convergence rates¶

A central question for any numerical method is its convergence rate: how fast does the error approach zero when the resolution is increased? For finite element methods, this typically corresponds to proving, theoretically or empirically, that the error \(e = {u_{\small\mbox{e}}} - u\) is bounded by the mesh size \(h\) to some power \(r\); that is, \(\|e\| \leq C h^r\) for some constant \(C\). The number \(r\) is called the convergence rate of the method. Note that different norms, like the \(L^2\)-norm \(\|e\|\) or \(H^1_0\)-norm \(\|\nabla e\|\) typically have different convergence rates.

To illustrate how to compute errors and convergence rates in FEniCS, we have included the function convergence_rate in ft11_poisson_bcs.py. This is a tool that is very handy when verifying finite element codes and will therefore be explained in detail here.

error = (u_e - u)**2*dx
E = sqrt(abs(assemble(error)))

u_e = Expression('sin(x[0])', degree=1)

u_e = Expression('sin(x[0])', degree=3)

E = errornorm(u_e, u, normtype='L2')

def errornorm(u_e, u):
    V = u.function_space()
    mesh = V.mesh()
    degree = V.ufl_element().degree()
    W = FunctionSpace(mesh, 'P', degree + 3)
    u_e_W = interpolate(u_e, W)
    u_W = interpolate(u, W)
    e_W = Function(W)
    e_W.vector()[:] = u_e_W.vector().array() - u_W.vector().array()
    error = e_W**2*dx
    return sqrt(abs(assemble(error)))

E = errornorm(u_e, u, norm_type='H10')

def convergence_rate(u_exact, f, u_D, kappa):
    """
    Compute convergence rates for various error norms for a
    sequence of meshes and elements.
    """

    h = {}  # discretization parameter: h[degree][level]
    E = {}  # error measure(s): E[degree][level][error_type]
    degrees = 1, 2, 3, 4
    num_levels = 5

    # Iterate over degrees and mesh refinement levels
    for degree in degrees:
        n = 4  # coarsest mesh division
        h[degree] = []
        E[degree] = []
        for i in range(num_levels):
            n *= 2
            h[degree].append(1.0 / n)
            u = solver(kappa, f, u_D, n, n, degree,
                       linear_solver='direct')
            errors = compute_errors(u_exact, u)
            E[degree].append(errors)
            print('2 x (%d x %d) P%d mesh, %d unknowns, E1=%g' %
                  (n, n, degree, u.function_space().dim(),
                   errors['u - u_exact']))

    # Compute convergence rates
    from math import log as ln  # log is a fenics name too
    error_types = list(E[1][0].keys())
    rates = {}
    for degree in degrees:
        rates[degree] = {}
        for error_type in sorted(error_types):
            rates[degree][error_type] = []
            for i in range(num_meshes):
                Ei   = E[degree][i][error_type]
                Eim1 = E[degree][i-1][error_type]
                r = ln(Ei/Eim1)/ln(h[degree][i]/h[degree][i-1])
                rates[degree][error_type].append(round(r,2))

    return rates

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

omega = 1.0
u_e = Expression('sin(omega*pi*x[0])*sin(omega*pi*x[1])', omega=omega)

f = 2*pi**2*omega**2*u_e

Taking advantage of structured mesh data¶

Many readers have extensive experience with visualization and data analysis of 1D, 2D, and 3D scalar and vector fields on uniform, structured meshes, while FEniCS solvers exclusively work with unstructured meshes. Since it can many times be practical to work with structured data, we discuss in this section how to extract structured data for finite element solutions computed with FEniCS.

A necessary first step is to transform our Mesh object to an object representing a rectangle with equally-shaped rectangular cells. The second step is to transform the one-dimensional array of nodal values to a two-dimensional array holding the values at the corners of the cells in the structured mesh. We want to access a value by its \(i\) and \(j\) indices, \(i\) counting cells in the \(x\) direction, and \(j\) counting cells in the \(y\) direction. This transformation is in principle straightforward, yet it frequently leads to obscure indexing errors, so using software tools to ease the work is advantageous.

In the directory src/modules, associated with this book, we have included the Python module BoxField that provides utilities for working with structured mesh data in FEniCS. Given a finite element function u, the following function returns a BoxField object that represents u on a structured mesh:

def structured_mesh(u, divisions):
    """Represent u on a structured mesh."""
    # u must have P1 elements, otherwise interpolate to P1 elements
    u2 = u if u.ufl_element().degree() == 1 else \
         interpolate(u, FunctionSpace(mesh, 'P', 1))
    mesh = u.function_space().mesh()
    from BoxField import fenics_function2BoxField
    u_box = fenics_function2BoxField(
        u2, mesh, divisions, uniform_mesh=True)
    return u_box

[AL 8: Comment below needs to be modified if interpolation to P1 is included in BoxField.]

Note that we can only turn functions on meshes with \(\mathsf{P}_1\) elements into BoxField objects, so if u is based on another element type, we first interpolate the scalar field onto a mesh with \(\mathsf{P}_1\) elements. Also note that to use the function, we need to know the divisions into cells in the various spatial directions (divisions).

The u_box object contains several useful data structures:

• u_box.grid: object for the structured mesh
• u_box.grid.coor[X]: grid coordinates in X=0 direction
• u_box.grid.coor[Y]: grid coordinates in Y=1 direction
• u_box.grid.coor[Z]: grid coordinates in Z=2 direction
• u_box.grid.coorv[X]: vectorized version of u_box.grid.coor[X] (for vectorized computations or surface plotting)
• u_box.grid.coorv[Y]: vectorized version of u_box.grid.coor[Y]
• u_box.grid.coorv[Z]: vectorized version of u_box.grid.coor[Z]
• u_box.values: numpy array holding the u values; u_box.values[i,j] holds u at the mesh point with coordinates

[AL 9: Is something missing in the last line above? The i and j indices?]

u = solver(p, f, u_b, nx, ny, 1, linear_solver='direct')
u_box = structured_mesh(u, (nx, ny))
u_ = u_box.values       # numpy array

print('u_ is defined on a structured mesh with %s points' %
      str(u_.shape))

# Iterate over 2D mesh points (i,j)
for j in range(u_.shape[1]):
    for i in range(u_.shape[0]):
        print('u[%d, %d] = u(%g, %g) = %g' %
              (i, j,
               u_box.grid.coor[X][i], u_box.grid.coor[X][j],
               u_[i,j]))

x = u_box.grid.coor[X]
dx = x[1] - x[0]
u_xx = (u_[i - 1, j] - 2*u_[i, j] + u_[i + 1, j]) / dx**2

Making surface plots¶

The ability to access a finite element field as structured data is handy in many occasions, e.g., for visualization and data analysis. Using Matplotlib, we can create a surface plot, as shown in Figure Various plots of the solution on a structured mesh (upper left):

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = fig.gca(projection='3d')
cv = u_box.grid.coorv  # vectorized mesh coordinates
ax.plot_surface(cv[X], cv[Y], u_, cmap=cm.coolwarm,
                rstride=1, cstride=1)
plt.title('Surface plot of solution')

The key issue is to know that the coordinates needed for the surface plot is in u_box.grid.coorv and that the values are in u_.

Various plots of the solution on a structured mesh

fig = plt.figure()
ax = fig.gca()
levels = [1.5, 2.0, 2.5, 3.5]
cs = ax.contour(cv[X], cv[Y], u_, levels=levels)
plt.clabel(cs)  # add labels to contour lines
plt.axis('equal')
plt.title('Contour plot of solution')

start = (0, 0.4)
x, u_val, y_fixed, snapped = u_box.gridline(start, direction=X)

start = (0, 0.4)
x, u_val, y_fixed, snapped = u_box.gridline(start, direction=X)
u_e_val = [u_b((x_, y_fixed)) for x_ in x]

plt.figure()
plt.plot(x, u_val, 'r-')
plt.plot(x, u_e_val, 'bo')
plt.legend(['P1 elements', 'exact'], loc='upper left')
plt.title('Solution along line y=%g' % y_fixed)
plt.xlabel('x');  plt.ylabel('u')

flux_u = flux(u, p)
flux_u_x, flux_u_y = flux_u.split(deepcopy=True)

# Plot the numerical and exact flux along the same line
if flux_u_x.ufl_element().degree() == 1:
   flux2_x = flux_u_x
else:
   V1 = FunctionSpace(u.function_space().mesh(), 'P', 1)
   flux2_x = interpolate(flux_x, V)
flux_u_x_box = structured_mesh(flux_u_x, (nx, ny))
x, flux_u_val, y_fixed, snapped = \
   flux_u_x_box.gridline(start, direction=X)
y = y_fixed

plt.figure()
plt.plot(x, flux_u_val, 'r-')
plt.plot(x, flux_u_x_exact(x, y_fixed), 'bo')
plt.legend(['P1 elements', 'exact'], loc='upper right')
plt.title('Flux along line y=%g' % y_fixed)
plt.xlabel('x')
plt.ylabel('u')

from sympy import exp, sin, pi  # for use in math formulas
import sympy as sym

H = lambda x: exp(-16*(x-0.5)**2)*sin(3*pi*x)
x, y = sym.symbols('x[0], x[1]')
u = H(x)*H(y)

u_c = sym.printing.ccode(u)

-exp(-16*pow(x[0] - 0.5, 2) - 16*pow(x[1] - 0.5, 2))*
sin(3*M_PI*x[0])*sin(3*M_PI*x[1])

u_c = u_c.replace('M_PI', 'pi')
u_b = Expression(u_c)

kappa = 1
f = sym.diff(-kappa*sym.diff(u, x), x) + \
    sym.diff(-kappa*sym.diff(u, y), y)
f = sym.simplify(f)
f_c = sym.printing.ccode(f)
f_c = f_c.replace('M_PI', 'pi')
f = Expression(f_c)

flux_u_x_exact = sym.lambdify([x, y], -kappa*sym.diff(u, x),
                              modules='numpy')