Thank you for your help and of course generally for your great code! I have found a similar solution which is a bit hacky but works for my setup, which I will share here for the case that anyone stumbles upon the same problem.

The code implements the projection exactly as one expects it: One approximates the tested integrals b_j =( \varphi_j,a(u))_\Gamma for {j=1,..., \text{dof}} with the quadrature module and then finds the desired coefficient vector by solving the system Mc =b (M denoting the mass matrix with all spaces set to the RT-space).

The following function has to be called once to calculate a dictionary, which maps elements to elementary gridfunctions (gridfunctions where only one coefficient is nonzero), which are nonzero on that element.

```
def precomp(space):
nontrivialEntries = []
dof = space.global_dof_count
import numpy as np
coeffs = np.zeros(dof)
element_list = list(space.grid.leaf_view.entity_iterator(0))
gridfunList = []
domainList = [[] for _ in element_list]
domainDict = dict(zip(element_list,[[] for _ in element_list]))
for j in range(dof):
coeffs[j] = 1
gridfun = bempp.api.GridFunction(space,coefficients = coeffs)
gridfunList.append(gridfun)
indices = np.nonzero(np.sum(gridfun.evaluate_on_element_centers()**2,axis = 0))
domainList[indices[0][0]].append(j)
domainList[indices[0][1]].append(j)
domainDict[element_list[indices[0][0]]].append(j)
domainDict[element_list[indices[0][1]]].append(j)
coeffs[j] = 0
identity = bempp.api.operators.boundary.sparse.identity(space,space,space)
id_weak = identity.weak_form()
id_sparse = aslinearoperator(id_weak).sparse_operator
neighborlist = id_sparse.tolil().rows
return gridfunList,neighborlist,domainDict
```

Then, passing the resulting arguments together with a Gridfunction and a nonlinearity into the following function will yield the desired GridFunction object.

```
def applyNonlinearity(gridFun,nonlinearity,gridfunList,domainDict):
space = gridFun.space
dof = space.global_dof_count
coeff = gridFun.coefficients
from bempp.api.integration import gauss_triangle_points_and_weights
accuracy_order = gridFun.parameters.quadrature.far.single_order
points, weights = gauss_triangle_points_and_weights(accuracy_order)
element_list = list(gridfunList[0].grid.leaf_view.entity_iterator(0))
weightIntegrals = np.zeros(dof)
for element in element_list:
integration_elements = element.geometry.integration_elements(points)
for i in domainDict[element]:
testFuncEvals = gridfunList[i].evaluate(element, points)
gFevals = gridFun.evaluate(element,points)
weightIntegrals[i] += sum([testFuncEvals[:,k].dot(nonlinearity(gFevals[:,k]))*weights[k]*integration_elements[k] for k in range(len(testFuncEvals[0,:]))] )
identity = bempp.api.operators.boundary.sparse.identity(RT_space, RT_space, RT_space)
id_weak = identity.weak_form()
from scipy.sparse.linalg import gmres
coeffsol,info = gmres(id_weak,weightIntegrals,tol=10**(-5))
return bempp.api.GridFunction(gridFun.space,coefficients = coeffsol)
```

A minimal example running the code is now given by:

```
RT_space = bempp.api.function_space(grid,"RT",0)
gridfunList,neighborlist,domainDict = precompMM(RT_space)
def a(x):
return np.linalg.norm(x)**2*x
def u(x, n, domain_index, result):
result[:] = np.cross(n,np.cross(np.array([np.exp(-1*(x[2])), 0. * x[2], 0. * x[2]]), n))
u_gridFun = bempp.api.GridFunction(RT_space, fun=tangential_trace,dual_space=RT_space)
a_of_u_gridFun = applyNonlinearity(trace_fun,a,gridfunList,domainDict)
```

Again thanks for your help.