Far-field scattering pattern for FEM-BEM coupling

Applications
1

Hi everyone,
I want to use FEniCS together with Bempp to compute the far-field scattering pattern from a unit sphere. I modified " Simple FEM-BEM coupling for the Helmholtz equation with FEniCS" for a unite sphere and run the code without error. Now I am wondering how I can compute the far-field pattern?

Thanks
-Hassan

2

Hassan,

The FEM-BEM tutorial plots the near field of the solution. If you’d like to calculate the far-field pattern, you can use the same representation formula (see markdown cell in tutorial).
The variables dirichlet_fun and neumann_fun store the Dirichlet and Neumann traces of the solution. Instead of using bempp.api.operators.potential.helmholtz.single_layer, use bempp.api.operators.far_field.helmholtz.single_layer and similar for the double-layer operator. The far-field value of the scattered field is then dlp_pot.evaluate(dirichlet_fun) - slp_pot.evaluate(neumann_fun). See the tutorial “The OSRC preconditioner for high-frequency scattering” for instructions to plot the far field in the correct points.

Best,
Elwin

1 Like
3

Dear Elwin,

Thank you very much for your reply. Based on your explanation, I modified the code. Here is the code for generating far-filed scattering from a unit sphere centered at the orgin:

# Simple FEM-BEM coupling for the Helmholtz equation with FEniCS for a
# unit sphere centered at the origin.

import dolfin
import bempp.api
import numpy as np
import pygalmesh
import meshio

k = 6.2841
d = np.array([1., 1., 1.])
d /= np.linalg.norm(d)

# mesh = dolfin.UnitCubeMesh(10, 10, 10)
s = pygalmesh.Ball([0, 0, 0], 1.0)
mesh = pygalmesh.generate_mesh(s, max_cell_circumradius=0.1)
meshio.write("sphere.xml", mesh)
mesh = dolfin.Mesh("sphere.xml")

from bempp.api.external import fenics

fenics_space = dolfin.FunctionSpace(mesh, "CG", 1)
trace_space, trace_matrix = \
    fenics.fenics_to_bempp_trace_data(fenics_space)
bempp_space = bempp.api.function_space(trace_space.grid, "DP", 0)

print("FEM dofs: {0}".format(mesh.num_vertices()))
print("BEM dofs: {0}".format(bempp_space.global_dof_count))

id_op = bempp.api.operators.boundary.sparse.identity(
    trace_space, bempp_space, bempp_space)
mass = bempp.api.operators.boundary.sparse.identity(
    bempp_space, bempp_space, trace_space)
dlp = bempp.api.operators.boundary.helmholtz.double_layer(
    trace_space, bempp_space, bempp_space, k)
slp = bempp.api.operators.boundary.helmholtz.single_layer(
    bempp_space, bempp_space, bempp_space, k)


u = dolfin.TrialFunction(fenics_space)
v = dolfin.TestFunction(fenics_space)
n = 0.25

@bempp.api.complex_callable
def u_inc(x, n, domain_index, result):
    result[0] = np.exp(1j * k * np.dot(x, d))
u_inc = bempp.api.GridFunction(bempp_space, fun=u_inc)

# The rhs from the FEM
rhs_fem = np.zeros(mesh.num_vertices())
# The rhs from the BEM
rhs_bem = u_inc.projections(bempp_space)
# The combined rhs
rhs = np.concatenate([rhs_fem, rhs_bem])

from bempp.api.assembly.blocked_operator import BlockedDiscreteOperator
from bempp.api.external.fenics import FenicsOperator
from scipy.sparse.linalg.interface import LinearOperator
blocks = [[None,None],[None,None]]

trace_op = LinearOperator(trace_matrix.shape, lambda x:trace_matrix @ x)

A = FenicsOperator((dolfin.inner(dolfin.nabla_grad(u),
                                 dolfin.nabla_grad(v)) \
    - k**2 * n**2 * u * v) * dolfin.dx)

blocks[0][0] = A.weak_form()
blocks[0][1] = -trace_matrix.T * mass.weak_form().to_sparse()
blocks[1][0] = (.5 * id_op - dlp).weak_form() * trace_op
blocks[1][1] = slp.weak_form()

blocked = BlockedDiscreteOperator(np.array(blocks))

from bempp.api.assembly.discrete_boundary_operator import InverseSparseDiscreteBoundaryOperator
from scipy.sparse.linalg import LinearOperator

# Compute the sparse inverse of the Helmholtz operator
# Although it is not a boundary operator we can use
# the SparseInverseDiscreteBoundaryOperator function from
# BEM++ to turn its LU decomposition into a linear operator.
P1 = InverseSparseDiscreteBoundaryOperator(
    blocked[0,0].to_sparse().tocsc())

# For the Laplace slp we use a simple mass matrix preconditioner. 
# This is sufficient for smaller low-frequency problems.
P2 = InverseSparseDiscreteBoundaryOperator(
    bempp.api.operators.boundary.sparse.identity(
        bempp_space, bempp_space, bempp_space).weak_form())

# Create a block diagonal preconditioner object using the Scipy LinearOperator class
def apply_prec(x):
    """Apply the block diagonal preconditioner"""
    m1 = P1.shape[0]
    m2 = P2.shape[0]
    n1 = P1.shape[1]
    n2 = P2.shape[1]
    
    res1 = P1.dot(x[:n1])
    res2 = P2.dot(x[n1:])
    return np.concatenate([res1, res2])

p_shape = (P1.shape[0] + P2.shape[0], P1.shape[1] + P2.shape[1])
P = LinearOperator(p_shape, apply_prec, dtype=np.dtype('complex128'))

# Create a callback function to count the number of iterations
it_count = 0
def count_iterations(x):
    global it_count
    it_count += 1

from scipy.sparse.linalg import gmres
soln, info = gmres(blocked, rhs, M=P, callback=count_iterations)

soln_fem = soln[:mesh.num_vertices()]
soln_bem = soln[mesh.num_vertices():]

print("Number of iterations: {0}".format(it_count))


# Solution function with dirichlet data on the boundary
dirichlet_data = trace_matrix * soln_fem
dirichlet_fun = bempp.api.GridFunction(trace_space, coefficients=dirichlet_data)

# Solution function with Neumann data on the boundary
neumann_fun = bempp.api.GridFunction(bempp_space, coefficients=soln_bem)


# We now want to plot the radar cross section in the z=0 plane. To compute it 
# we use the far-field operators implemented in Bempp.

theta = np.linspace(0, 2 * np.pi, 400)
points = np.array([np.cos(theta), np.sin(theta), np.zeros(len(theta))])

slp_far_field = bempp.api.operators.far_field.helmholtz.single_layer(
    bempp_space, points, k)
dlp_far_field = bempp.api.operators.far_field.helmholtz.double_layer(
    trace_space, points, k)

far_field = dlp_far_field.evaluate(dirichlet_fun) - slp_far_field.evaluate(neumann_fun)

max_incident = np.abs(u_inc.coefficients).max()
radiation_pattern = (np.abs(far_field/max_incident)**2).ravel()
db_pattern = 10 * np.log10(4 * np.pi * radiation_pattern)


from matplotlib import pyplot as plt
plt.polar(theta, db_pattern)
plt.ylim(db_pattern.min(), db_pattern.max())
plt.show()

Best wishes,
-Hassan

1 Like