Challenges Encountered When Implementing Relevant Literature Methods Using BEMPP

Applications
1

Hello everyone!

I have read most of the relevant articles on BEMPP and noticed that Matthew’s doctoral thesis provides a detailed description of BEMPP. However, when I try to implement the direct EFIE, indirect EFIE, MFIE, and CFIE from the thesis using bempp, I consistently encounter issues with either non-convergent results or incorrect outcomes (Only the indirect EFIE format converges.).

I would be grateful if you could advise on whether there might be problems in my code or in the theoretical implementation process. Please find the original text and my implementation process below:

61188ecb0c4249f0fea3bc0b41aad41d
61188ecb0c4249f0fea3bc0b41aad41d1080×3221 426 KB

The code is as follows:

import bempp_cl.api
import numpy as np
from bempp_cl.api.operators.boundary import maxwell

frequency = 3e8
eps0 = 8.854187817e-12
mu0 = 4 * np.pi * 1e-7
zf0 = np.sqrt(mu0/eps0)
k = 2*np.pi*frequency*np.sqrt(eps0*mu0)

polarization = np.array([-1.0, 0, 0])
direction = np.array([0, 0, -1.0])

def incident_field(point):
    return polarization * np.exp(1j * k * np.dot(point, direction))

@bempp_cl.api.complex_callable#(jit=False)
def tangential_trace(point, n, domain_index, result):
    value = polarization * np.exp(1j * k * np.dot(point, direction))
    result[:] = np.cross(value, n)

@bempp_cl.api.complex_callable    
def neumann_trace(point, n, domain_index, result):
    value = np.cross(direction, polarization) * 1j * k * np.exp(1j * k * np.dot(point, direction))
    result[:] = 1.0 / (1j * k) * np.cross(value, n)
grid = bempp_cl.api.shapes.regular_sphere(4)
grid.plot()  

tan_inc = bempp_cl.api.GridFunction(rwg_space, fun = tangential_trace, dual_space = snc_space)
neu_inc = bempp_cl.api.GridFunction(rwg_space, fun = neumann_trace, dual_space = snc_space)

E = maxwell.electric_field(rwg_space, rwg_space, snc_space, k)
E_ik = maxwell.electric_field(rwg_space, rwg_space, snc_space, 1j*k)
H = maxwell.magnetic_field(rwg_space, rwg_space, snc_space, k)
Id = bempp_cl.api.operators.boundary.sparse.identity(rwg_space, rwg_space, snc_space)

rwg_space = bempp_cl.api.function_space(grid, "RWG", 0)
snc_space = bempp_cl.api.function_space(grid, "SNC", 0)

The following is the format of the algorithm in the paper

sol= bempp_cl.api.linalg.lu(E, (0.5 * Id + H) * tan_inc)   #direct EFIE 
#sol= bempp_cl.api.linalg.lu(E, tan_inc)   #indirect EFIE
# sol= bempp_cl.api.linalg.lu((H + 0.5 * Id), -E * tan_inc)   #direct MFIE
# sol= bempp_cl.api.linalg.lu((H - 0.5 * Id), tan_inc)   #indirect MFIE
# sol= bempp_cl.api.linalg.gmres(-E_ik * E + 0.5 * Id + H, -E_ik*(0.5 * Id + H) * tan_inc)   #direct CFIE

# Number of points in the x-direction
nx = 300

# Number of points in the y-direction
ny = 300

# Generate the evaluation points with numpy
# x, y, z = np.mgrid[-3 : 3 : nx * 1j, -3 : 3 : ny * 1j,  0 : 0 : 1j     ] #xoy
x, y, z = np.mgrid[-3 : 3 : nx * 1j,  0 : 0 : 1j     , -3 : 3 : ny * 1j] #xoz
# x, y, z = np.mgrid[ 0 : 0 : 1j     , -3 : 3 : nx * 1j, -3 : 3 : ny * 1j] #yoz面
points = np.vstack((x.ravel(), y.ravel(), z.ravel()))

# Compute interior and exterior indices
all_indices = np.ones(points.shape[1], dtype="uint32")

interior_indices = np.sum((points) ** 2, axis=0) < 1 ** 2
exterior_indices = ~interior_indices

ext_points = points[:, exterior_indices]
int_points = points[:, interior_indices]

ME = bempp_cl.api.operators.potential.maxwell.electric_field(rwg_space, ext_points, k)
MH = bempp_cl.api.operators.potential.maxwell.magnetic_field(rwg_space, ext_points, k)

exterior_values = MH * tan_inc - ME * sol                      #dirct EFIE & direct MFIE & direct CFIE
#exterior_values = -ME * sol                                  #indirct EFIE
# fexterior_values = -MH * sol                                 #indirct MFIE

Thank you for your time and assistance. I look forward to your reply and wish you all the best.