"""
Richard S. Falk, February 24, 2014; updated February 23, 2017

A 1D example with an essential boundary condition at both end points

- u'' + u = e - e^{-1}, 0 < x < 1,   u(0) = 0,  u(1) = 0

The exact solution is given by

u(x) = (e^{-1} -1)*(e^x-1) + (1-e)*(e^{-x}-1)

Study of convergence in the L2 and H1 norms, and max error at the mesh points.

"""
from dolfin import *
from numpy import log, exp

nmeshes = 5  # number of meshes to compute with
deg =1     # Polynomial order of trial/test functions
errors = []  # list into which to store the errors
# Define the exact solution for later use
uexact = Expression('(exp(-1.0)-1.0)*(exp(x[0])-1.0) + (1-exp(1.0))*(exp(-x[0])-1.0)',degree=deg+2)

# Define the boundary
xl = 0 # Left end point
xr = 1 # Right end point

for i in range(nmeshes):
    n = 2**(i+2)
    # Create a uniform mesh of n elements on the interval (xl, xr)
    mesh = IntervalMesh(n,xl,xr)
    # Define function space for this mesh using Continuous Galerkin
    # (Lagrange) functions of order deg on each element
    V = FunctionSpace(mesh, 'CG', deg)
    # Define a boundary value at the left and right end points
    ub = Constant(0.)
    # Define the portion of the boundary at which the Dirichlet BC is imposed
    def Dirichlet_boundary(x,on_boundary):
        return x[0] < DOLFIN_EPS or x[0] > 1 - DOLFIN_EPS and on_boundary
        
    # The following statement enforces the Dirichlet BCs
    # on the portion of the boundary defined by Dirichlet_boundary
    # It is defined using the function space V, the value ub, for points
    # for which Dirichlet_boundary is true.
    bc = DirichletBC(V,ub,Dirichlet_boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(exp(1.0) - exp(-1.0))
    a = inner(grad(u),grad(v))*dx + u*v*dx
    L = f*v*dx

    # Solve the discrete equations for the function uh in V
    uh = Function(V)  # declares uh as a finite element function in the space V
    solve(a == L, uh, bc)

    # compare with exact solution interpolated into high degree space Vex
    Vex = FunctionSpace(mesh, 'CG', deg+2)
    uinterp = interpolate(uexact, Vex)
    err = uinterp - interpolate(uh, Vex)
    # also compare with interpolant of the exact solution in V
    ui = interpolate(uexact, V)
    erri = ui - uh
    # Compute L2 norm of uinterp
    L2norm = sqrt( assemble( uinterp*uinterp*dx ) )
    # Compute H1 seminorm of uinterp (just the derivative part of the norm)
    H1seminorm = sqrt( assemble( inner(grad(uinterp),grad(uinterp))*dx ) )
    # Compute the Max norm of ui
    uiv = ui.vector()
    Maxnorm = norm(uiv,'linf')
    # Compute the L2 norm of the error
    L2normerr = sqrt( assemble( err*err*dx ) )
    # Compute H1 seminorm of the error
    H1seminormerr = sqrt( assemble( inner(grad(err),grad(err))*dx ) )
    # Compute the maximum error at the mesh points
    erriv = ui.vector() - uh.vector()
    Maxerr = norm(erriv,'linf')
    # Append this latest result to the list in which the errors are stored.
    errors.append([1.0/n, L2normerr, H1seminormerr, Maxerr])

print "\n     h      L2 error           H1 error          Max error   \
    L2 rate H1 rate Max rate \n"
print "  {:7.5f}   {:4.2e} ({:5.2f}%)  {:4.2e} ({:5.2f}%) {:4.2e} ({:5.2f}%)" \
    .format(errors[0][0], errors[0][1], 100*errors[0][1]/L2norm, \
    errors[0][2], 100*errors[0][2]/H1seminorm, errors[0][3], \
    100*errors[0][3]/Maxnorm)
for i in range(1,nmeshes):
    print "  {:7.5f}   {:4.2e} ({:5.2f}%)  {:4.2e} ({:5.2f}%)  {:4.2e} ({:5.2f}%) {:5.2f} {:5.2f} {:5.2f}".format(errors[i][0], errors[i][1], \
        100*errors[i][1]/L2norm, errors[i][2], 100*errors[i][2]/H1seminorm, \
        errors[i][3], 100*errors[i][3]/Maxnorm, \
        log(errors[i-1][1]/errors[i][1])/log(2), \
        log(errors[i-1][2]/errors[i][2])/log(2), \
        log(errors[i-1][3]/errors[i][3])/log(2))

# plot solution
plot(uh)
# hold plot 
interactive()  # close window to get program to continue