Hi,
I try to model electromagnetic radiation from a patch antenna. In particular, I am not sure on the considerations for computing correctly the fields when a dielectric is present.
Would appreciate help on the definition of spaces, boundaries, solver and field data, to compute the total field? Code below:
#################################################################################
# Boundary element method solver for electromagnetic fields - Bempp
try:
import bempp.api
except ImportError:
print("installing bempp ...")
!pip install bempp-cl --quiet
!pip install meshio --quiet
!pip install gmsh --quiet
print("installed bempp.")
import bempp.api
LocalComputing=False
if (LocalComputing):
bempp.api.GMSH_PATH = 'C:/WPy64-39100/python-3.9.10.amd64/Scripts/gmsh'#bempp.api.utils.which("gmsh")
import numba
def BemppApiGridUnion(grids, domain_indices=None, swapped_normals=None):
"""
Return the union of a given list of grids.
Parameters
----------
grids: list
A list of grid objects.
domain_indices : list
Attach a list of domain indices to the new grid such that grid[j] received the domain index domain_indices[j]
swapped_normals : list of boolean
A list of the form [False, True, ...], that specifies for each grid if the normals should be swapped (True) or not (False). This is helpful if one grid is defined to be inside another grid.
This method returns a new grid object, which is the union of the input grid objects.
"""
vertex_offset = 0
element_offset = 0
vertex_count = sum([grid.number_of_vertices for grid in grids])
element_count = sum([grid.number_of_elements for grid in grids])
vertices = np.empty((3, vertex_count), dtype="float64")
elements = np.empty((3, element_count), dtype="uint32")
all_domain_indices = np.empty(element_count, dtype="uint32")
if domain_indices is None:
domain_indices = range(len(grids))
if swapped_normals is None:
swapped_normals = len(grids) * [False]
for index, grid in enumerate(grids):
nelements = grid.number_of_elements
nvertices = grid.number_of_vertices
vertices[:, vertex_offset : vertex_offset + nvertices] = grid.vertices
if swapped_normals[index]:
current_elements = grid.elements[[0, 2, 1], :]
else:
current_elements = grid.elements
elements[:, element_offset : element_offset + nelements] = (current_elements + vertex_offset)
all_domain_indices[element_offset : element_offset + nelements] = domain_indices[index]
vertex_offset += nvertices
element_offset += nelements
return bempp.api.grid.grid.Grid(vertices, elements, all_domain_indices)
##############################################################################
# Patch antenna
OperationalFrequency=3e9 # Hz
OperationWavelengthFreeSpace=c_0/OperationalFrequency
wavenumberk=2.0*np.pi/OperationWavelengthFreeSpace
#Patch antenna dielectric
RelativePermittivityPA=3.0
OperationWavelengthPA=c_0/OperationalFrequency/np.sqrt(RelativePermittivityPA)
wavenumberkPA=2.0*np.pi/OperationWavelengthPA
GridSizeh=0.1*np.min([OperationWavelengthFreeSpace,OperationWavelengthPA]) # h : float; A floating point number specifying the grid size. Very critically computation value. Dimensions and positions of interest should be larger than this value.
NxPoints=64 #
NyPoints=NxPoints
NzPoints=128
Xlength=10.0*OperationWavelengthFreeSpace
Ylength=Xlength
Zrange=22.0*OperationWavelengthFreeSpace
[xMatrix, yMatrix, zMatrix] = np.mgrid[-Xlength/2.0:Xlength/2.0:NxPoints*1j,-Ylength/2.0:Ylength/2.0:NyPoints*1j,-2*Zrange/22.0:20*Zrange/22.0:NzPoints*1j]#np.mgrid[-extent:extent:nx * 1j,0:0:1j,-extent:extent:nz * 1j]
xCoord=np.linspace(-Xlength/2.0,Xlength/2.0,NxPoints)
yCoord=np.linspace(-Ylength/2.0,Ylength/2.0,NyPoints)
zCoord=np.linspace(-2*Zrange/22.0,20*Zrange/22.0,NzPoints)
PointsVstack = np.vstack((xMatrix.ravel(), yMatrix.ravel(), zMatrix.ravel()))
lenPointsVstack=PointsVstack.shape[1]
###################################################
# Compute interior and exterior indices
all_indices = np.ones(PointsVstack.shape[1], dtype='uint32')
RectangleMaskX=np.logical_and((PointsVstack[0]-0.0)<0.5*OperationWavelengthPA,(PointsVstack[0]-0.0)>-0.5*OperationWavelengthPA)
RectangleMaskY=np.logical_and((PointsVstack[1]-0.0)<0.5*OperationWavelengthPA,(PointsVstack[1]-0.0)>-0.5*OperationWavelengthPA)
RectangleMaskZ=np.logical_and((PointsVstack[2]+0.1*OperationWavelengthPA)<0.1*OperationWavelengthPA,(PointsVstack[2]+0.1*OperationWavelengthPA)>-0.1*OperationWavelengthPA)
RectangleMask=RectangleMaskX*RectangleMaskY*RectangleMaskZ
interior_indices = np.where(RectangleMask)
exterior_indices = np.where(RectangleMask==False) #~interior_indices0 & ~interior_indices1
int_points = PointsVstack[:, interior_indices]
ext_points = PointsVstack[:, exterior_indices]
##########################
# Geometry
plate1 = np.array([[0.5*OperationWavelengthPA, 0.5*OperationWavelengthPA, -0.2*OperationWavelengthPA],[-0.5*OperationWavelengthPA, 0.5*OperationWavelengthPA, -0.2*OperationWavelengthPA],[-0.5*OperationWavelengthPA, -0.5*OperationWavelengthPA, -0.2*OperationWavelengthPA],[0.5*OperationWavelengthPA, -0.5*OperationWavelengthPA, -0.2*OperationWavelengthPA]])
plate2 = np.array([[0.25*OperationWavelengthPA, 0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA],[-0.25*OperationWavelengthPA, 0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA],[-0.25*OperationWavelengthPA, -0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA],[0.25*OperationWavelengthPA, -0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA]])
gridDielectricRect = bempp.api.shapes.cuboid(length=(1.0*OperationWavelengthPA,1.0*OperationWavelengthPA,0.2*OperationWavelengthPA), origin=(0, 0.0, -0.1*OperationWavelengthPA), h=GridSizeh)
grid1 = bempp.api.shapes.screen(plate1,h=GridSizeh)
grid2 = bempp.api.shapes.screen(plate2,h=GridSizeh)
# BemppApiGridUnion - Not good that elements touch each other
grid = BemppApiGridUnion([grid1,grid2,gridDielectricRect])
# We define the spaces of order 0 RWG div-conforming functions and order 0 scaled Nédélec curl-conforming functions.
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)
### Source driven
#Edges to Identify
edgeV1 = [0.0*OperationWavelengthPA, 0.125*OperationWavelengthPA, -0.2*OperationWavelengthPA]
edgeV2 = [0.0*OperationWavelengthPA, 0.125*OperationWavelengthPA, 0.0*OperationWavelengthPA]
#Get Edges
edgeInd=[] # A list because maybe more than one edge is found
for i, e in enumerate(grid.edges.T):
if ((np.abs(grid.vertices[:, e[0]]-edgeV1)<=GridSizeh*0.5).all() or (np.abs(grid.vertices[:, e[0]]-edgeV2)<=GridSizeh*0.5).all()):
if ((np.abs(grid.vertices[:, e[1]]-edgeV1)<=GridSizeh*0.5).all() or (np.abs(grid.vertices[:, e[1]]-edgeV2)<=GridSizeh*0.5).all()):
edgeInd.append(i)
# Get DOF
dof=edgeInd
print('edgeInd: '+str(edgeInd))
for iIter in range(0,len(edgeInd),1):
for edges, dofs in zip(grid.element_edges.T, div_space.local2global):
for e, d in zip(edges, dofs):
if e == edgeInd[iIter]:
dof[iIter] = d
print('dof: '+str(dof))
#Set-Up Coefficients
coefficients=np.zeros((div_space.global_dof_count),dtype=np.complex64)
coefficients[dof] = 1.0+1j*0.0 # Source phasor. I belive it is a voltage source since the operators are electric field operators
# trace_fun is important to define the feeding:
# - an incident field: probably fun=tangential_trace
# - a source driven element: probably fun=coefficients - https://bempp.discourse.group/t/how-to-use-bempp-for-source-driven-em-problems/25/3
# Boundaries - for electromagnetic fields
# Next, we define the Maxwell electric field boundary operator and the identity operator. For Maxwell problems, the domain and range spaces should be div-conforming, while the dual_to_range space should be curl conforming
A0_int = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, wavenumberkPA)
A0_ext = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, wavenumberk)
identity = bempp.api.operators.boundary.sparse.identity(div_space, div_space, curl_space)
# Operators might be linearly combined - For instance one can dominate over the other?
A = ((np.sqrt(RelativePermittivityPA)-1.0)/np.sqrt(RelativePermittivityPA))*A0_int+(1.0/np.sqrt(RelativePermittivityPA))*A0_ext#np.sqrt(RelativePermittivityPA)*A0_int+A0_ext#bempp.api.GeneralizedBlockedOperator([[A0_int,identity],[A0_ext,identity]])
trace_fun = bempp.api.GridFunction(div_space, coefficients=coefficients, dual_space=curl_space)
rhs = trace_fun
################################
# Solver
lambda_data_sol = bempp.api.linalg.lu(A, rhs)
# Evaluate the fields
# Since there are no magnetic parts, we do not need the magnetic operator
slp_pot_int = bempp.api.operators.potential.maxwell.electric_field(div_space, PointsVstack, wavenumberkPA)
slp_pot_ext = bempp.api.operators.potential.maxwell.electric_field(div_space, PointsVstack, wavenumberk)
# Normalization to the free space perspective
interior_values=-(1.0/np.sqrt(RelativePermittivityPA))*slp_pot_int*lambda_data_sol#[0]
exterior_values=-slp_pot_ext*lambda_data_sol#[1]
# Now compute the total field
field_data = np.empty((3, PointsVstack.shape[1]), dtype='complex128')
field_data[:, interior_indices] = interior_values[:, interior_indices]
field_data[:, exterior_indices] = exterior_values[:, exterior_indices]