Support multiple Neumann boundary conditions per domain and provide a tutorial on this feature

[BenjaminRodenberg]

Our adapter only supports a single Neumann boundary condition on the whole domain (see FEniCS tutorial book, p.85). We add g*ds in the weak form to implement Neumann boundary conditions.

To be able to support multiple Neumann boundary conditions (for example, one coupling boundary condition and one on the original boundary), we have to use markers to restrict the boundary integral to the respective subdomain. Pseudocode:

coupling_domain.mark(boundary_markers, 1)
F += g*ds(1)

For more information see FEniCS tutorial book, section 4.4