Coefficients of composition of function and gridfunction

Applications
1

When implementing a solver for a nonlinear BIE formulation (by Newton’s method) I stumbled upon the following problem:
Given a nonlinear function a:\mathbb R^3\rightarrow \mathbb R^3 and coefficients c\in \mathbb R^{\text{dof}} corresponding to a gridfunction u in the RT-space, I would like to calculate the coefficients of the projection of the function a(u) in the RT-space.

My current Bempp version is 3.5.5, which I can not update to bempp-cl for unrelated reasons (though to my knowledge RT-spaces are anyhow unavailable in bempp-cl).

I would be grateful for any suggestions, references or ideas on how to approach this problem!

2

The only way I can think of doing this would be to define a new GridFunction using callable that evauates u at each point then applies a to the result. To evaluate u, though, you’d need to find which cell the point is in and its local coordinates. Something like this maybe:

def g(x, n, domain, result):
    for i, x_i in enumerate(x):
        cell_index = ... # Cell that x_i is in
        local_coords = ... # Compute local coordinates of x_i in this cell
        result[0][i] = f.evaluate(cell_index, local_coords)

result = bempp.api.GridFunction(space, fun=g)

I’m not convinced this is the best way, and I can’t think of a nice way to do the missing steps without looping through each cell and checking if the point is in that cell.

It’s also worth adding that RT spaces are now available in Bempp-cl.

3

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.

2 Likes