#!/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>

# Weakly imposing a Dirichlet boundary condition
# ==========================
# 
# This tutorial shows how to implement the weak imposition of a Dirichlet boundary condition, as proposed in the paper <a href='https://bempp.com/publications/#weakimposition'>Boundary Element Methods with Weakly Imposed Boundary Conditions (2019)</a>.

# First, we import Bempp and NumPy.

# In[1]:


import bempp.api
import numpy as np


# Next, we define the grid for our problem, and the function spaces that we will use. In this example, we use a sphere with P1 and DUAL0 function spaces.

# In[2]:


h = 0.3
grid = bempp.api.shapes.sphere(h=h)
p1 = bempp.api.function_space(grid, "B-P", 1)
dual0 = bempp.api.function_space(grid, "DUAL", 0)


# Next, we define the blocked operators proposed in the paper:
# $$
# \left(\begin{pmatrix}-\mathsf{K}&\mathsf{V}\\\mathsf{W}&\mathsf{K}'\end{pmatrix}+
# \begin{pmatrix}\tfrac12\mathsf{Id}&0\\\beta\mathsf{Id}&-\tfrac12\mathsf{Id}\end{pmatrix}\right)
# \begin{pmatrix}u\\\lambda\end{pmatrix}
# =
# \begin{pmatrix}g_\text{D}\\\beta g_\text{D}\end{pmatrix},
# $$
# where $\beta>0$ is a parameter of our choice. In this example, we use $\beta=0.1$.

# In[3]:


beta = 0.1
multi = bempp.api.BlockedOperator(2,2)
multi[0,0] = -bempp.api.operators.boundary.laplace.double_layer(p1, p1, dual0)
multi[0,1] = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0)
multi[1,0] = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1)
multi[1,1] = bempp.api.operators.boundary.laplace.adjoint_double_layer(dual0, dual0, p1)

diri = bempp.api.BlockedOperator(2,2)
diri[0,0] = 0.5 * bempp.api.operators.boundary.sparse.identity(p1, p1, dual0)
diri[1,0] = beta * bempp.api.operators.boundary.sparse.identity(p1, dual0, p1)
diri[1,1] = -0.5 * bempp.api.operators.boundary.sparse.identity(dual0, dual0, p1)


# Next, we define the function $g_\text{D}$, and define the right hand side.
# 
# Here, we use $$g_\text{D}=\sin(\pi x)\sin(\pi y)\sinh(\sqrt2\pi z),$$ as in section 5 of the paper.

# In[4]:


def f(x, n, d, res):
    res[0] = np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2])

f_fun = bempp.api.GridFunction(p1, fun=f)

rhs = [2*diri[0,0]*f_fun, diri[1,0]*f_fun]


# Now we solve the system. We set `use_strong_form=True` to activate mass matrix preconditioning.

# In[5]:


sol, info, it_count = bempp.api.linalg.gmres(multi+diri, rhs, return_iteration_count=True, use_strong_form=True)
print("Solution took",it_count,"iterations")


# For this problem, we know the analytic solution. We compute the error in the $\mathcal{B}_\text{D}$ norm:

# In[6]:


def g(x, n, d, res):
    grad = np.array([
            np.cos(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,
            np.sin(np.pi*x[0]) * np.cos(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,
            np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.cosh(np.sqrt(2)*np.pi*x[2]) * np.pi * np.sqrt(2)
        ])
    res[0] = np.dot(grad, n)

g_fun = bempp.api.GridFunction(dual0, fun=g)

e_fun = [sol[0]-f_fun,sol[1]-g_fun]

error = 0
# V norm
slp = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0)
hyp = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1)
error += np.sqrt(np.dot(e_fun[1].coefficients.conjugate(),(slp * e_fun[1]).projections(dual0)))
error += np.sqrt(np.dot(e_fun[0].coefficients.conjugate(),(hyp * e_fun[0]).projections(p1)))
# D part
error += beta**.5 * e_fun[0].l2_norm()

print("Error:",error)