Solving a mixed Neumann-Dirichlet Problem

Download this tutorial as a Jupyter notebook

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 \displaystyle \Gamma_D the part of the boundary that holds the Dirichlet boundary conditions and by \displaystyle \Gamma_N the boundary part that holds the Neumann boundary conditions. We denote by \displaystyle t\in\Gamma_D the unknown Neumann data and by \displaystyle u\in\Gamma_N the unknown Dirichlet data. The given Dirichlet data on \displaystyle \Gamma_D is denoted by \displaystyle g_D and the given Neumann data on \displaystyle \Gamma_N is denoted by \displaystyle g_N.

From Green’s representation theorem it follows that

\displaystyle   \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\\[2mm]  \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

Here (as usual) \displaystyle \mathsf{V}, \displaystyle \mathsf{K}, \displaystyle \mathsf{K}', \displaystyle \mathsf{W} are the single layer, double layer, adjoint double layer and hypersingular boundary operators.

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 \displaystyle \Gamma_N and \displaystyle \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.

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.

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 \displaystyle t on \displaystyle \Gamma_D. For \displaystyle \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 \displaystyle \Gamma_D. This is important for dofs on boundary edges and excludes associated functions that lie just outside \displaystyle \Gamma_D on the other side of the boundary edge.
  • The neumann_space_neumann_segment space is defined on \displaystyle \Gamma_N. \displaystyle \Gamma_N is open: the boundary edges are not part of the space. We again restrict basis functions to \displaystyle \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 \displaystyle \Gamma_D. The space is closed and by default basis functions are allowed to extend into the elements adjacent to \displaystyle \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 \displaystyle \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.
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.

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.

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.

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.

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 \displaystyle \Gamma_N and \displaystyle \Gamma_D into the global Dirichlet and Neumann spaces.

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.