#!/usr/bin/env python
# coding: utf-8

# <a href='http://bempp.com' style='line-height:50px'><img src='bempp-logo.jpg' style='float:left;height:50px;margin-right:10px'> This tutorial forms part of the documentation for Bempp. Find out more at bempp.com</a>

# # 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.

# In[1]:


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.

# In[2]:


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.

# In[3]:


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.

# In[4]:


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.

# In[5]:


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.

# In[6]:


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.

# In[7]:


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.

# In[8]:


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. <img src="cube_mixed_solution.png">