Douglas N. Arnold 2014

The first problem we considered was the PDE

Only slightly more complicated than homogeneous Neumann boundary conditions are the inhomogeneous Neumann boundary conditions

To implement the new term in FEniCS is straightforward. First we have to define an `Expression`

for

In []:

```
g = Expression('pow(x[0], 2)*x[1]')
```

Next we add the extra term to the definition of

In []:

```
F = f*v*dx + g*v*ds
```

The remainder of the program need not be changed.

*Remark:* Note that, in addition to the "measure" `dx`

which is used to signal integration over the domain, here we use the measure `ds`

which means integration over the boundary of the domain. It is also
possible to single out a part of the domain or a part of the boundary
and integrate over that with a measure like `dx(1)`

or `ds(3)`

. We do not need this yet, but will for more complicated boundary conditions (see below).

Instead of the Neumann boundary condition, we can impose the Dirichlet boundary condition *essential boundary condition*, meaning that it is imposed as a constraint on the space

FEniCS imposes these constraints this after assembling the stiffness
matrix and load vector, by modifying entries associated to boundary
DOFs. For this, we need to specify to FEniCS which DOFs are to be
constrained (in this case all the DOFs on the boundary), and what they
are to be constrained to (namely

`DirichletBC(Vh, g, bdry)`

where the three arguments are the function space to constrain, the
function to constrain to, and the portion of the domain where DOFs are
to be constrained. The function `Expression`

just as for the Neumann boundary condition. When the Dirichlet boundary condition is to be imposed on the whole boundary of `bdry = DomainBoundary()`

. The relevant code is then:

In []:

```
g = Expression('...')
bc = DirichletBC(Vh, g, DomainBoundary())
uh = Function(Vh)
solve(b == F, uh, bc)
```

Note that the `solve`

call now takes a third argument, specifying the essential boundary condition.

Sometimes we have different boundary conditions on different, complementary, parts of the boundary. For example, let

`DomainBoundary()`

, but must define the subset `DirichletBC`

. We might choose to call the function `GammaD`

. It takes two arguments, point `x`

, and a Boolean on_boundary. When FEniCS needs to identify whether some point of the domain belongs to `GammaD(x, on_boundary)`

with `x`

set to the point (the array of its coordinates `x[0]`

and `x[1]`

), and `on_boundary`

set to `True`

or `False`

according to whether the point `True`

or `False`

according to whether the point `near(x[1], 0.)`

for this. It returns `True`

if `x[1]`

is within round off error of zero. Thus we can specify
In []:

```
def GammaD(x, on_boundary):
return near(x[1], 0)
```

From this point on the code is similar to that discussed above.

In []:

```
g = Expression('...')
bc = DirichletBC(Vh, g, GammaD)
uh = Function(Vh)
solve(b == F, uh, bc)
```

Note that for this problem, with `F = Constant(0.)*v*dx`

.

*Remark.* People are often confused about what to do at the
interface where the boundary conditions switch from Neumann to
Dirichlet. In our example, this would consist of the two points

It is possible to have multiple Dirichlet boundary conditions, say `DirichletBC`

, each just as above, say `bc1`

and `bc2`

, and then call solve as

`solve(b == F, uh, [bc1, bc2])`

When inhomogeneous Neumann conditions are imposed on part of the boundary, we may need to include an integral like `ds(1)`

which integrates only over `MeshFunction`

which assigns an integer to each edge, then define the relevant portion of the boundary using a subclass of the FEniCS class `SubDomain`

,
then use this subclass to set the values of the mesh function, and then
use the mesh function to make the new measure. Here is some typical
code for this, which redefines `ds`

so that `ds(1)`

integrates over the left half of the boundary (where the first coordinate is `ds(0)`

integrates over the right half of the boundary, and `ds()`

integrates over the whole boundary.

In []:

```
# create a mesh function which assigns an unsigned integer (size_t) to each edge
mf = MeshFunction("size_t", mesh, 1) # 3rd argument is dimension of an edge
mf.set_all(0) # initialize the function to zero
# create a SubDomain subclass, specifying the portion of the boundary with x[0] < 1/2
class Lefthalf(SubDomain):
def inside(self, x, on_boundary):
return x[0] < 0.5 and on_boundary
lefthalf = Lefthalf() # instantiate it
# use this lefthalf object to set values of the mesh function to 1 in the subdomain
lefthalf.mark(mf, 1)
# define a new measure ds based on this mesh function
ds = Measure("ds")[mf]
```

Robin boundary conditions take the form

In the case of periodic boundary conditions, the boundary condition does not constrain individual points on the boundary, but relates the values at two different points. We will not discuss this further except to say that this possibility is available in FEniCS.

In [5]:

```
# this cell is present to set the notebook style
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/dna1.css", "r").read()
return HTML(styles)
css_styling()
```

Out[5]: