# Computing the capacity of a cube with a re-entrant corner

### Background

The capacity $\text{cap}(\Omega)$ of an isolated conductor $\Omega\subset\mathbb{R}^3$ with boundary $\Gamma$ measures its ability to store charges. It is defined as the ratio of the total surface equilibrium charge relative to its surface potential value. To compute the capacity, we need to solve the following exterior Laplace problem for the equilibrium potential $u$ with unit surface value:
$$
\begin{align}
-\Delta u &= 0\quad\text{in }\Omega^\text{+},\\
u &= 1\quad\text{on }\Gamma,\\
|u(\mathbf{x})| &=\mathcal{O}\left(|\mathbf{x}|^{-1}\right)\quad\text{as }|\mathbf{x}|\rightarrow\infty.
\end{align}
$$
Here $\Omega^\text{+}$ is the domain exterior to $\Omega$.
The total surface charge of an isolated conductor is given by Gauss law as
$$
\text{cap}(\Omega)=-\epsilon_0\int_{\Gamma}\frac{\partial u}{\partial\nu}(\mathbf{x})\,\mathrm{d}\mathbf{x}.
$$
Here, $\nu(\mathbf{x})$ is the outward pointing normal direction for $\mathbf{x}\in\Gamma$, and $\epsilon_0$ is the electric constant with value $\epsilon_0\approx 8.854\times 10^{-12}\,{\rm F/m}$. For simplicity, in the following we will exclude this constant and use the expression $\text{cap}(\Omega)=-\int_{\Gamma}\partial u/\partial\nu\,d\mathbf{x}$.

Using Green's representation theorem and noting that the exterior Laplace double layer potential is zero for constant densities, we can represent the solution $u$ as
$$
u(x) = -\int_{\Gamma} g(\mathbf{x},\mathbf{y})\phi(\mathbf{y})\,\mathrm{d}\mathbf{y} \quad\text{for all }\mathbf{x}\in\Omega^\text{+},
$$
where $\phi:={\partial u}/{\partial\nu}$ is the normal derivative of the exterior solution $u$ and $g(\mathbf{x},\mathbf{y}):=\frac{1}{4\pi|\mathbf{x}-\mathbf{y}|}$ is the Green's function of the 3D Laplacian. By taking boundary traces, we arrive at the following boundary integral equation of the first kind.
$$
1 = -\int_{\Gamma} g(\mathbf{x},\mathbf{y})\phi(\mathbf{y})\,\mathrm{d}\mathbf{y} =: -\mathsf{V}\phi(\mathrm{x})\quad\text{for all }\mathbf{x}\in\Gamma.
$$
The capacity is now simply given by
$$
\text{cap}(\Omega) = -\int_\Gamma \phi(\mathbf{x}) \,\mathrm{d}\mathbf{x}.
$$

To improve the convergence we will solve the preconditioned equation
$$
\mathsf{V}\tilde{\mathsf{W}}\phi = -1,
$$
where $\tilde{\mathsf{W}}$ is a regularised hypersingular operator defined by the weak form
$$
\langle \tilde{\mathsf{W}}w, v\rangle = \langle \mathsf{W}w, v\rangle + \langle w, 1\rangle\langle v, 1\rangle.
$$

### Implementation

We start with the usual imports.


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

The grid re-entrant cube is predefined in the shapes module. By default it refines towards the singular corner.


```python
grid = bempp.api.shapes.reentrant_cube(h=0.05, refinement_factor=1)
```

Next, we define the right-hand side and the operator. We discretise a pair of hypersingular and single layer operators. Using the corresponding Bempp function guarantees that only one single layer operator needs to be discretised. The parameter `dual` means that we discretise the single-layer on the dual grid using piecewise constant basis functions and the hypersingular operator on the original grid using piecewise linear continuous basis functions. The regularisation is done via a built-in rank one operator that is added to the hypersingular operator.


```python
def one_fun(x, n, domain_index, res):
    res[0] = -1

slp, hyp, base_slp = \
    bempp.api.operators.boundary.laplace.single_layer_and_hypersingular_pair(
        grid, spaces='dual', return_base_slp=True)  

rank_one_op = bempp.api.RankOneBoundaryOperator(
    hyp.domain, hyp.range, hyp.dual_to_range)
hyp_regularized = hyp + rank_one_op
    
rhs_fun = bempp.api.GridFunction(slp.range, fun=one_fun)

lhs = slp * hyp_regularized
```

We use GMRES to solve the system. To improve convergence we use a strong form discretisation that automatically preconditions with the mass matrix.


```python
x, info, it_count = bempp.api.linalg.gmres(
    lhs, rhs_fun, use_strong_form=True,
    return_iteration_count=True)

sol_fun = hyp_regularized * x

print("Number of iterations: {0}".format(it_count))

print("The capacity is {0}.".format(-sol_fun.integrate()[0,0]))
```

    Number of iterations: 6
    The capacity is 8.1104817108.


We now want to compute local residual error estimates for the computed solution. To do this, we use the already calculates single layer operator on the barycentrically refined grid.

First, we compute the map from element indices of the original grid to descendent element indices on the barycentric refinement. ``bary_map[i,:]`` are the indices of the elements in the barycentric refinement that are associated with the element with index ``i`` on the original grid. 


```python
bary_map = grid.barycentric_descendents_map()
```

For the error estimator we want to compute element diameters. This is handled by the following function.


```python
def diameter(element):
    """Compue the diameter of an element."""
    corners = element.geometry.corners
    d0 = np.linalg.norm(corners[:,0] - corners[:,1])
    d1 = np.linalg.norm(corners[:,0] - corners[:, 2])
    d2 = np.linalg.norm(corners[:,1] - corners[:, 2])
    return np.max([d0, d1, d2])
```

We now evaluate the solution on the barycentrically refined grid. The function space on the barycentric refinement is the space of piecewise linear, discontinuous functions. This was used automatically by Bempp to assemble the base single layer operator on the barycentric refinement.

We define an identity mapping the solution into this space. 


```python
map_to_base_space = bempp.api.operators.boundary.sparse.identity(sol_fun.space, base_slp.domain, base_slp.domain)
base_sol_fun = map_to_base_space * sol_fun

base_slp_times_sol = base_slp * base_sol_fun
```

We now evaluate the squared error estimate $\eta_{T}^2 = \text{diam}(T)\|\nabla \mathsf{V}\phi\|_{\mathcal{L}^2(T)}^2$ on each element $T$ of the barycentric refinement, and sum of the contributions of the barycentric triangles associated with each triangle of the original grid. We can ignore the right-hand side function as it is constant and therefore its surface gradient is zero.


```python
bary_grid = grid.barycentric_grid()

local_errors_squared_bary = np.zeros(
    bary_grid.leaf_view.entity_count(0), dtype='float64')
local_errors = np.zeros(
    grid.leaf_view.entity_count(0), dtype='float64')
bary_index_set = bary_grid.leaf_view.index_set()

# Compute the error estimates on the barycentric refinement

for element in bary_grid.leaf_view.entity_iterator(0):
    index = bary_index_set.entity_index(element)
    gradient_norm = base_slp_times_sol.surface_grad_norm(element)
    local_error_squared = diameter(element) * gradient_norm**2
    local_errors_squared_bary[index] = local_error_squared

# Sum up the local contributions for each triangle to obtain an error estimate on the
# original grid.
    
for m in range(bary_map.shape[0]):
    for n in range(bary_map.shape[1]):
        local_errors[m] += local_errors_squared_bary[bary_map[m, n]]

local_errors = np.sqrt(local_errors)
total_error = np.linalg.norm(local_errors)

print("Total H^(-1/2) error estimate for normal derivative: {0}".format(
        np.linalg.norm(local_errors)))
print("Error estimate for capacity: {0}".format(total_error**2))
```

    Total H^(-1/2) error estimate for normal derivative: 0.0
    Error estimate for capacity: 0.0


We now want to plot the errors on each element. We first have to reorder the vector with elementwise errors to fit to the indices of a piecewise constant function space. The numbering of the two is not necessarily identical and depends on the underlying grid manager. We use the ``index_set`` to access the numberings of elements.


```python
const_space = bempp.api.function_space(grid, "DP", 0)

# Resort the error contributions
sorted_local_errors = np.zeros_like(local_errors)
index_set = grid.leaf_view.index_set()
for element in grid.leaf_view.entity_iterator(0):
    index = index_set.entity_index(element)
    global_dof_index = const_space.get_global_dofs(element, dof_weights=False)[0]
    sorted_local_errors[global_dof_index] = local_errors[index]

bempp.api.GridFunction(const_space, coefficients=sorted_local_errors).plot()
```



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