Interactions missing with _curl_curl_product for multishape geometries.

[pescap]

Dear all,

I recently found out that the _curl_curl_product operator was not properly evaluated when using multishape geometries (i.e., in the presence of segments and swapped_normals).

I have attached a notebook reproducing the behavior. I could not find a similar issue with other operators (EFIE, MFIE, identities, Laplace-Beltrami, and the other operators involved in the OSRC preconditioner for Maxwell) as they behave as expected.

Any ideas?

import bempp.api
import numpy as np

def describe(c):
    cr, ci = np.real(c), np.imag(c)
    print(cr.max(), cr.min(), cr.mean())
    print(ci.max(), ci.min(), ci.mean())

grid = bempp.api.shapes.multitrace_cube(h=1)

segments = [[1, 2, 3, 4, 5, 6]]
swapped_normals = [[]]

grid1=bempp.api.shapes.cuboid(length = (1,1,0.5), origin=(0,0,0), h=1)

#grid.plot()
#grid1.plot()

k0 = 1
rwg = [bempp.api.function_space(grid, "RWG", 0, segments=seg, swapped_normals=normals, include_boundary_dofs=True)
              for seg, normals in zip(segments, swapped_normals)]

snc = [bempp.api.function_space(grid, "SNC", 0, segments=seg, swapped_normals=normals,include_boundary_dofs=True)
              for seg, normals in zip(segments, swapped_normals)]

p1dA= [bempp.api.function_space(grid, "DP", 1, segments=seg, swapped_normals=normals,include_boundary_dofs=True)
              for seg, normals in zip(segments, swapped_normals)]

rwg1 = bempp.api.function_space(grid1, "RWG", 0,include_boundary_dofs=True)
snc1 = bempp.api.function_space(grid1, "SNC", 0,include_boundary_dofs=True)
p1d1 = bempp.api.function_space(grid1, "DP", 1,include_boundary_dofs=True)

efiem = bempp.api.as_matrix(bempp.api.operators.boundary.maxwell.electric_field(rwg[0], rwg[0], snc[0], k0).weak_form())
efie1m = bempp.api.as_matrix(bempp.api.operators.boundary.maxwell.electric_field(rwg1, rwg1, snc1, k0).weak_form())

print('describe the EFIE')
describe(efiem)
print('---')
describe(efie1m)
print('---')

I = bempp.api.operators.boundary.sparse.identity(snc[0], snc[0], snc[0])
I1 = bempp.api.operators.boundary.sparse.identity(snc1, snc1, snc1)

Im = bempp.api.as_matrix(I.weak_form())
I1m = bempp.api.as_matrix(I1.weak_form())

print('describe the Identity for snc/snc')
describe(Im)
print('---')
describe(I1m)
print('---')

CC = bempp.api.operators.boundary.sparse._curl_curl_product(snc[0], snc[0], snc[0])
CC1 = bempp.api.operators.boundary.sparse._curl_curl_product(snc1, snc1, snc1)

CCm = bempp.api.as_matrix(CC.weak_form())
CC1m = bempp.api.as_matrix(CC1.weak_form())

print('describe the curl_curl_product')
describe(CCm)
print('---')
describe(CC1m)
print('---')


from matplotlib import pyplot as plt
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12,6))

ax1.spy(CCm)
ax2.spy(CC1m)

[pescap]

The plt.spy shows that some terms for the multishape case (LEFT) are equal to zero (refer, e.g., to the diagonal terms in the left-hand figure):