I have touched on this topic before, without much satisfaction. Lately, I have taken a look at BEMPP python API and how it interfaces with Fenics. Unfortunately, for the Maxwell problem, only the lowest order polynomial interpolation (Nedelec or RT/RWG) is present. I would like to use the higher orders of hierachical basis functions that Fenics permits to do the boundary operator.
The boundary operator (which augments the Fenics variational form) is B=S+D:
S_{mn} = -k_0^2\iint \left ( \iint g(\mathbf{r},\mathbf{r}^\prime)(\mathbf{v}_m\times\mathbf{n})\cdot (\mathbf{u}_n\times\mathbf{n}^\prime) dS^\prime \right )dS ,
D_{mn} = \iint \left ( \iint g(\mathbf{r},\mathbf{r}^\prime)\nabla_\Gamma\cdot(\mathbf{v}_m\times\mathbf{n})\nabla_{\Gamma^\prime}\cdot (\mathbf{u}_n\times\mathbf{n}^\prime) dS^\prime \right )dS ,
where
g(\mathbf{r, r^\prime})=\frac{\exp(jk_0\vert\mathbf{ r-r}^\prime\vert)}{4\pi\vert\mathbf{r-r}^\prime\vert} . The primed coordinates are the source (trial function) and unprimed coordinates are tied to the testing functions. The \Gamma subscript indicates differentiation with respect to the surface coordinates.
I would like to insert the matrix elements from these integrals into the Fenics matrix data structure. Is there a way to access the finite element matrix elements to permit this augmentation? I realize that all this means getting into the nitty-gritty of low-level stuff, but it is quite difficult to find complete descriptions on how to do these manipulations. Can someone in the community point me in a suitable direction?
I should say that the number of surface DoFs is likely to be small w.r.t. the overall problem size - perhaps 1000-5000 DoFs in total out of several million, so sparsity is not affected inordinately.
Perhaps I am being unrealistic - feel free to tell me so!
Cheers and happy holidays!