I would like to implement a natural boundary conditions (i. e. such that result in right side contributions via boundary integrals) that are different on different parts on the boundary. The whole contribution by the boundary conditions is given by sum like \int_{\Gamma_1} ( \vec f_1, \vec n) \mathrm{d}s + \int_{\Gamma_2} ( \vec f_2, \vec n) \mathrm{d}s

Originally, I used a piece of code like

```
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
```

taken from the older fenics example to construct measures corresponding to \Gamma_i. This approach is unfortunately not working in my installation of 2018.1.0 as `FacetFunction`

is no longer available. I would be grateful for any advice on the preffered reimplementation in the current version.