# Combining Dirichlet and Neumann conditions#

Author: Jørgen S. Dokken

Let’s return to the Poisson problem from the Fundamentals chapter and see how to extend the mathematics and the implementation to handle Dirichlet condition in combination with a Neumann condition. The domain is still the unit square, but now we set the Dirichlet condition \(u=u_D\) at the left and right sides, while the Neumann condition

is applied to the remaining sides \(y=0\) and \(y=1\).

## The PDE problem#

Let \(\Lambda_D\) and \(\Lambda_N\) denote parts of the boundary \(\partial \Omega\) where the Dirichlet and Neumann conditions apply, respectively. The complete boundary-value problem can be written as

Again, we choose \(u=1+x^2+2y^2\) as the exact solution and adjust \(f, g,\) and \(u_D\) accordingly

For the ease of programming, we define \(g\) as a function over the whole domain \(\Omega\) such that \(g\) takes on the correct values at \(y=0\) and \(y=1\). One possible extension is

## The variational formulation#

The first task is to derive the variational formulatin. This time we cannot omit the boundary term arising from integration by parts, because \(v\) is only zero on \(\Lambda_D\). We have

and since \(v=0\) on \(\Lambda_D\),

by applying the boundary condition on \(\Lambda_N\). The resulting weak from reads

Expressing this equation in the standard notation \(a(u,v)=L(v)\) is straight-forward with

## Implementation#

As in the previous example, we define our mesh,function space and bilinear form \(a(u,v)\).

```
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad
import numpy as np
import pyvista
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
```

Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fulfill this condition.

```
def u_exact(x):
return 1 + x[0]**2 + 2 * x[1]**2
def boundary_D(x):
return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
dofs_D = locate_dofs_geometrical(V, boundary_D)
u_bc = Function(V)
u_bc.interpolate(u_exact)
bc = dirichletbc(u_bc, dofs_D)
```

The next step is to define the Neumann condition. We first define \(g\) uses `UFL`

s `SpatialCoordinate`

-function, and then in turn create a boundary integration measure `ds`

. As the test function \(v\) is zero on the boundary integrals over the Dirichlet boundary disappears, and we can integrate `g*v*ds`

over the entire boundary.

```
x = SpatialCoordinate(mesh)
g = -4 * x[1]
f = Constant(mesh, default_scalar_type(-6))
L = f * v * dx - g * v * ds
```

We can now assemble and solve the linear system of equations

```
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V2 = functionspace(mesh, (
```