# Solving a mixed Neumann-Dirichlet Problem

### Background

With Bempp, it is possible to define operators only on segments of a given domain. This makes it possible to solve mixed Neumann-Dirichlet problems. In this tutorial, we solve the Laplace equation inside the unit cube with unit Dirichlet boundary conditions on two sides and unit Neumann boundary conditions on the other four sides.

Denote by $\Gamma_D$ the part of the boundary that holds the Dirichlet boundary conditions and by $\Gamma_N$ the boundary part that holds the Neumann boundary conditions. We denote by $t\in\Gamma_D$ the unknown Neumann data and by $u\in\Gamma_N$ the unknown Dirichlet data. The given Dirichlet data on $\Gamma_D$ is denoted by $g_D$ and the given Neumann data on $\Gamma_N$ is denoted by $g_N$.

From Green's representation theorem it follows that
$$
\begin{align}
\left[\mathsf{V}t\right] (\mathbf{x}) - \left[\mathsf{K}u\right] (\mathbf{x}) &= \left[\tfrac{1}{2}\mathsf{Id} + \mathsf{K}\right]g_D(\mathbf{x}) - \mathsf{V}g_N(\mathbf{x}),\quad \mathbf{x}\in\Gamma_D\\
\left[\mathsf{W}u\right] (\mathbf{x}) + \left[\mathsf{K}'t\right] (\mathbf{x}) &=\left[\tfrac{1}{2}\mathsf{Id} - \mathsf{K}'\right]g_N(\mathbf{x}) - \mathsf{W}g_D(\mathbf{x}),\quad \mathbf{x}\in\Gamma_N
\end{align}
$$
Here (as usual) $\mathsf{V}$, $\mathsf{K}$, $\mathsf{K}'$, $\mathsf{W}$ are the single layer, double layer, adjoint double layer and hypersingular <a href='https://bempp.com/2017/07/11/available_operators/'>boundary operators</a>.

The difficulty in the implementation is the definition of the discrete function spaces and the treatment of degrees of freedom (dofs) that lie on the interface between $\Gamma_N$ and $\Gamma_D$. In the following, we will go through the implementation and point out how to correctly define all spaces involved.

### Implementation

We start with the usual imports. In addition we increase the integration order, as in this example we will be working with spaces of quadratic functions.


```python
import bempp.api
import numpy as np

bempp.api.global_parameters.quadrature.medium.double_order = 4
bempp.api.global_parameters.quadrature.far.double_order = 4
```

We now define the domain. We use a standard unit cube. In the corresponding function all sides of the cube are already associated with different domain indices. We associate the indices 1 and 3 with the Dirichlet boundary and the other indices with the neumann boundary.


```python
grid = bempp.api.shapes.cube()
dirichlet_segments = [1, 3]
neumann_segments = [2, 4, 5, 6]
```

We can now define the spaces. For the Neumann data, we use discontinuous polynomial basis functions of order 1. For the Dirichlet data, we use continuous basis functions of local polynomial order 2.

We need global spaces for the Dirichlet and Neumann data and suitable spaces on the segments. The space definitions are as follows:

* The ``neumann_space_dirichlet_segment`` space holds the unknown Neumann data $t$ on $\Gamma_D$. For $\Gamma_D$ we use the parameter ``closed=True``, meaning that all boundary edges and the associated dofs on the boundary edges are part of the space. The parameter ``element_on_segment=True`` implies that we restrict functions to elements that lie on elements associated with $\Gamma_D$. This is important for dofs on boundary edges and excludes associated functions that lie just outside $\Gamma_D$ on the other side of the boundary edge.

* The ``neumann_space_neumann_segment`` space is defined on $\Gamma_N$. $\Gamma_N$ is open: the boundary edges are not part of the space. We again restrict basis functions to $\Gamma_N$ by the parameter ``element_on_segment=True``. However, we also include all functions which are defined on elements of the space but whose reference points (i.e. the dof positions) are on the excluded boundary. This is achieved by the parameter ``reference_point_on_segment=False``. If it were set to ``True`` (default) it would only include dofs whose reference points lie in the segment and not on the excluded boundary.

* The ``dirichlet_space_dirichlet_segment`` space is a space of continuous basis functions that holds the Dirichlet data on $\Gamma_D$. The space is closed and by default basis functions are allowed to extend into the elements adjacent to $\Gamma_D$. This extension is necessary because of the definition of the underlying Sobolev space on the segment. To control this behavior for continuous spaces the option ``strictly_on_segment`` exists, which is by default set to ``False``.

* The ``dirichlet_space_neumann_segment`` is defined similarly to the ``dirichlet_space_dirichlet_segment`` but on the open segment $\Gamma_N$.

* For the discretisation of the Dirichlet data, we also need the space ``dual_dirichlet_space``. This is the correct dual space for projecting functions into the space of Dirichlet data.


```python
order_neumann = 1
order_dirichlet = 2


global_neumann_space = bempp.api.function_space(grid, "DP", order_neumann)
global_dirichlet_space = bempp.api.function_space(grid, "P", order_dirichlet)

neumann_space_dirichlet_segment = bempp.api.function_space(
    grid, "DP", order_neumann, domains=dirichlet_segments,
    closed=True, element_on_segment=True)

neumann_space_neumann_segment = bempp.api.function_space(
    grid, "DP", order_neumann, domains=neumann_segments,
    closed=False, element_on_segment=True, reference_point_on_segment=False)

dirichlet_space_dirichlet_segment = bempp.api.function_space(
    grid, "P", order_dirichlet, domains=dirichlet_segments, closed=True)

dirichlet_space_neumann_segment = bempp.api.function_space(
    grid, "P", order_dirichlet, domains=neumann_segments, closed=False)

dual_dirichlet_space = bempp.api.function_space(
    grid, "P", order_dirichlet, domains=dirichlet_segments,
    closed=True, strictly_on_segment=True)
```

In the following, we define all operators on the corresponding spaces and the overall blocked operator.


```python
slp_DD = bempp.api.operators.boundary.laplace.single_layer(
    neumann_space_dirichlet_segment,
    dirichlet_space_dirichlet_segment,
    neumann_space_dirichlet_segment)

dlp_DN = bempp.api.operators.boundary.laplace.double_layer(
    dirichlet_space_neumann_segment,
    dirichlet_space_dirichlet_segment,
    neumann_space_dirichlet_segment)

adlp_ND = bempp.api.operators.boundary.laplace.adjoint_double_layer(
    neumann_space_dirichlet_segment,
    neumann_space_neumann_segment,
    dirichlet_space_neumann_segment)

hyp_NN = bempp.api.operators.boundary.laplace.hypersingular(
    dirichlet_space_neumann_segment,
    neumann_space_neumann_segment,
    dirichlet_space_neumann_segment)

slp_DN = bempp.api.operators.boundary.laplace.single_layer(
    neumann_space_neumann_segment,
    dirichlet_space_dirichlet_segment,
    neumann_space_dirichlet_segment)

dlp_DD = bempp.api.operators.boundary.laplace.double_layer(
    dirichlet_space_dirichlet_segment,
    dirichlet_space_dirichlet_segment,
    neumann_space_dirichlet_segment)

id_DD = bempp.api.operators.boundary.sparse.identity(
    dirichlet_space_dirichlet_segment,
    dirichlet_space_dirichlet_segment,
    neumann_space_dirichlet_segment)

adlp_NN = bempp.api.operators.boundary.laplace.adjoint_double_layer(
    neumann_space_neumann_segment,
    neumann_space_neumann_segment,
    dirichlet_space_neumann_segment)

id_NN = bempp.api.operators.boundary.sparse.identity(
    neumann_space_neumann_segment,
    neumann_space_neumann_segment,
    dirichlet_space_neumann_segment)

hyp_ND = bempp.api.operators.boundary.laplace.hypersingular(
    dirichlet_space_dirichlet_segment,
    neumann_space_neumann_segment,
    dirichlet_space_neumann_segment)

blocked = bempp.api.BlockedOperator(2, 2)

blocked[0, 0] = slp_DD
blocked[0, 1] = -dlp_DN
blocked[1, 0] = adlp_ND
blocked[1, 1] = hyp_NN
```

Next, we define the functions of the Dirichlet and Neumann data and their discretisations on the corresponding segments.


```python
def dirichlet_data_fun(x):
    return 1
    
def dirichlet_data(x, n, domain_index, res):
    res[0] = dirichlet_data_fun(x)
    
def neumann_data_fun(x):
    return 1
 
def neumann_data(x, n, domain_index, res):
    res[0] = neumann_data_fun(x)

dirichlet_grid_fun = bempp.api.GridFunction(
    dirichlet_space_dirichlet_segment,
    fun=dirichlet_data, dual_space=dual_dirichlet_space)

neumann_grid_fun = bempp.api.GridFunction(
    neumann_space_neumann_segment,
    fun=neumann_data, dual_space=dirichlet_space_neumann_segment)

rhs_fun1 = (.5 * id_DD + dlp_DD) * dirichlet_grid_fun \
           - slp_DN * neumann_grid_fun
rhs_fun2 = - hyp_ND * dirichlet_grid_fun \
           + (.5 * id_NN - adlp_NN) * neumann_grid_fun
```

We can now discretise and solve the blocked operator system. We solve without preconditioner. This would cause problems if we were to further increase the degree of the basis functions.


```python
lhs = blocked.weak_form()
rhs = np.hstack([rhs_fun1.projections(neumann_space_dirichlet_segment), 
                 rhs_fun2.projections(dirichlet_space_neumann_segment)])

from scipy.sparse.linalg import gmres
x, info = gmres(lhs, rhs)
```

Next, we split up the solution vector and define the grid functions associated with the computed Neumann and Dirichlet data.


```python
nx0 = neumann_space_dirichlet_segment.global_dof_count

neumann_solution = bempp.api.GridFunction(
    neumann_space_dirichlet_segment, coefficients=x[:nx0])
dirichlet_solution = bempp.api.GridFunction(
    dirichlet_space_neumann_segment, coefficients=x[nx0:])
```

We want to recombine the computed Dirichlet and Neumann data with the corresponding known data in order to get Dirichlet and Neumann grid functions defined on the whole grid. To achieve this we define identity operators from $\Gamma_N$ and $\Gamma_D$ into the global Dirichlet and Neumann spaces.


```python
neumann_imbedding_dirichlet_segment = \
    bempp.api.operators.boundary.sparse.identity(
        neumann_space_dirichlet_segment,
        global_neumann_space,
        global_neumann_space)

neumann_imbedding_neumann_segment = \
    bempp.api.operators.boundary.sparse.identity(
        neumann_space_neumann_segment,
        global_neumann_space,
        global_neumann_space)

dirichlet_imbedding_dirichlet_segment = \
    bempp.api.operators.boundary.sparse.identity(
        dirichlet_space_dirichlet_segment,
        global_dirichlet_space,
        global_dirichlet_space)

dirichlet_imbedding_neumann_segment = \
    bempp.api.operators.boundary.sparse.identity(
        dirichlet_space_neumann_segment,
        global_dirichlet_space,
        global_dirichlet_space)

dirichlet = (dirichlet_imbedding_dirichlet_segment * dirichlet_grid_fun +
             dirichlet_imbedding_neumann_segment * dirichlet_solution)

neumann = (neumann_imbedding_neumann_segment * neumann_grid_fun +
           neumann_imbedding_dirichlet_segment * neumann_solution)

dirichlet.plot()
```

We can plot the solution using the command ``dirichlet.plot()``. The solution looks as follows. 

1fa89f6eb9f5c045f185dbe639d23712b1ccd4c1
ipynb
Jupyter Notebook
notebooks/mixed_neumann_dirichlet.ipynb
mscroggs/bempp-tutorials
4ab135f6bad73c56f76c6d93cf0d01e61e683fd9
MIT
notebooks/mixed_neumann_dirichlet.ipynb
mscroggs/bempp-tutorials
4ab135f6bad73c56f76c6d93cf0d01e61e683fd9
MIT
notebooks/mixed_neumann_dirichlet.ipynb
mscroggs/bempp-tutorials
4ab135f6bad73c56f76c6d93cf0d01e61e683fd9
MIT
Qwen/Qwen-72B
1. YES
2. YES