# Scaling functions and symmetries in forms

Hello all.
I’ve been scratching my head over this during the last few days. In what seems like a straightforward calculation of the far field approximation for a surface current
\mathbf{E}_{ff} \propto \mathbf{a}_r\times\int\limits_{\partial\Omega}\exp(jk\mathbf{a}_r\cdot \mathbf{r}^\prime )\,\mathbf{n}^\prime\times\mathbf{E}\,dS^\prime ,
it has turned into a perplexing challenge to make this integral work in Dolfinx. I have two problems:

1. There is the complex scaling factor (\exp(\ldots)) that looks straightforward to implement given some of the methods that I have found in this forum.
2. I would like to use symmetries, where I generate solutions over a fraction of the geometry and by adding together scaled versions of the solution that take into account these symmetries, generate the full far-field integral.

I first ran the problem of a quarter sphere and applied symmetries on \mathbf{n}\times\mathbf{E} to form the full result.

import sys
import numpy as np
from dolfinx import default_scalar_type
from mpi4py import MPI as nMPI
from dolfinx import fem
import ufl
import basix.ufl
from dolfinx.io import gmshio
from dolfinx import io
import gmsh

class Fields:
def __init__(self, A, k0):
self.A = A
self.k0 = k0
def E(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Er = 377.0 * self.A * st * np.vstack([ct * cp, ct * sp, -st]) * np.exp(-1j * self.k0 * r) / r
return Er
def H(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Hr = self.A * st * np.vstack([-sp, cp, 0.0 * x[2]]) * np.exp(-1j * self.k0 * r) / r
return Hr

eps = 1.0e-4
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
if mpiRank == 0:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)

gmsh.model.occ.addBox(0, -10, 0, 10, 20, 10, 3)
gmsh.model.occ.intersect([(3,1)],[(3,3)], 4, removeObject=True, removeTool=True)
gmsh.model.occ.cut([(3,4)], [(3,2)], 5, removeObject=True, removeTool=True)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(0)
gmsh.model.mesh.setSize(pt, 0.8)
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
for bx in pt:
M = gmsh.model.occ.getMass(bx[0], bx[1])
CoM = gmsh.model.occ.getCenterOfMass(bx[0], bx[1])
if CoM[2] > 1.0+eps and CoM[0] > 2.0:
else:
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 36)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber("Mesh.Tetrahedra", 0)
gmsh.option.setNumber("Mesh.Smoothing", 10)

gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Gmsh")

#    gmsh.fltk.run()
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, 0, gdim=3)

if mpiRank == 0:
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
E = fem.Function(V)
H = fem.Function(V)
U = fem.Function(V)
k = 2.0 * np.pi / 5
# Generate a generic azimuthally symmetric field function
Fld = Fields(1.0, 2.0 * k)
E.interpolate(Fld.E)
H.interpolate(Fld.H)

SymmType = 1
# Plot field vectors
Vw = basix.ufl.element('DG', mesh.basix_cell(), 0, shape=(mesh.geometry.dim, ))
W = fem.functionspace(mesh, Vw)
Et = fem.Function(W)
Et.interpolate(E)
Et.name = "ElectricField"
with io.XDMFFile(mesh.comm, "Efield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
H_expr = fem.Expression(ufl.curl(E), W.element.interpolation_points())
Et.interpolate(H)
Et.name = "MagneticField"
with io.XDMFFile(mesh.comm, "Hfield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
# Scalar "phase" function space
Q = fem.functionspace(mesh, ('CG', 2))
vv = fem.Function(Q)

n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
a = 0.0
# Field symmetry condition about the x=0 plane for M = -n x E
Mcx = fem.Constant(mesh, default_scalar_type(((1.0, 0.0, 0.0),(0.0, -1.0, 0.0),(0.0, 0.0, -1.0)))) # x = 0 wall is perfect magnetic wall.
Mcz = fem.Constant(mesh, default_scalar_type(((1.0, 0.0, 0.0),(0.0, 1.0, 0.0),(0.0, 0.0, -1.0)))) # z = 0 plane is perfect elec. conductor
# Generate integral over quarter sphere surface, using symmetry to get value over full hemisphere surface for z>0
for p in range(26):
theta = np.pi * p / 50
for q in range(101):
phi = np.pi * q / 50
L= []
# Observation vector at position theta, phi (usual sphetrical coords)
rr = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
# U contains constant a_phi unit vector for observation point
U.interpolate(lambda x: (-np.sin(phi) + 0.0 * x[0], np.cos(phi) + x[1] * 0.0, 0.0 * x[2]))
# This is the scalar phase term
vv.interpolate(lambda x: np.exp(1j * k * (rr[0] * x[0] + rr[1] * x[1] + rr[2] * x[2])))
# Add to integrand original quarter sphere for x>0, z>0
L.append(vv * ufl.cross(n,E))
# The other quarter sphere using symmetry.  Phase factor for x < 0, z > 0. PMC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (-rr[0] * x[0] + rr[1] * x[1] + rr[2] * x[2])))
L.append(Mcx * vv * ufl.cross(n,E))
# The third quarter sphere using symmetry.  Phase factor for x > 0, z < 0. PEC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (rr[0] * x[0] + rr[1] * x[1] - rr[2] * x[2])))
L.append(Mcz * vv * ufl.cross(n,E))
# The last quarter sphere using symmetry.  Phase factor for x < 0, z < 0. PEC and PMC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (-rr[0] * x[0] + rr[1] * x[1] - rr[2] * x[2])))
L.append(Mcz * Mcx * vv * ufl.cross(n,E))
#
# Do integration of the four sphere slices over quarter sphere surface.  Ev should contain the azimuthal component of total field at point theta,phi
Ev = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(sum(L), U) * ds(1))), op=nMPI.SUM)
# Do the same integration for theta component.  Should be almost zero over all theta, phi.
U.interpolate(lambda x: (np.cos(theta)*np.cos(phi)+x[0]*0.0, np.cos(theta)*np.sin(phi)+x[1]*0.0, -np.sin(theta)+x[2]*0.0))
Eh = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(sum(L), U) * ds(1))), op=nMPI.SUM)
# Print out result
if mpiRank == 0:
print(theta, phi, np.absolute(Ev), np.absolute(Eh), Ev, Eh)
# The source current (-n x E) is fully azimuthally symmetric, so the list of Ev's should also exhibit this symmetry,
# and the Eh should be close to zero, but they are not. Why? Problem is with the symmetry matrix multiplication.  How to do this
# correctly?
sys.exit(0)


(Note 2024-08-30: fixed the rad boundary condition so ds(1) integration is only over outer spherical surface. Solution still not correct.)

The results over \theta and \phi exhibit none of the properties that I expected. At this point, I thought the problem was with the invocation of the symmetry conditions or with the phase factor weighting function. To eliminate the symmetry requirement, I ran an example with a full sphere.

import sys
import numpy as np
from dolfinx import default_scalar_type
from mpi4py import MPI as nMPI
from dolfinx import fem
import ufl
import basix.ufl
from dolfinx.io import gmshio
from dolfinx import io
import gmsh

class Fields:
def __init__(self, A, k0):
self.A = A
self.k0 = k0
def E(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Er = 377.0 * self.A * st * np.vstack([ct * cp, ct * sp, -st]) * np.exp(-1j * self.k0 * r) / r
return Er
def H(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Hr = self.A * st * np.vstack([-sp, cp, 0.0 * x[2]]) * np.exp(-1j * self.k0 * r) / r
return Hr

eps = 1.0e-4
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
if mpiRank == 0:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)

gmsh.model.occ.addBox(-10, -10, -10, 20, 20, 20, 3)
gmsh.model.occ.intersect([(3,1)],[(3,3)], 4, removeObject=True, removeTool=True)
gmsh.model.occ.cut([(3,4)], [(3,2)], 5, removeObject=True, removeTool=True)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(0)
gmsh.model.mesh.setSize(pt, 0.8)
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
for bx in pt:
M = gmsh.model.occ.getMass(bx[0], bx[1])
if M > 4.0 * np.pi + eps:
else:
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 36)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber("Mesh.Tetrahedra", 0)
gmsh.option.setNumber("Mesh.Smoothing", 10)

gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Gmsh")

#    gmsh.fltk.run()
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, 0, gdim=3)

if mpiRank == 0:
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
E = fem.Function(V)
H = fem.Function(V)
U = fem.Function(V)
k = 2.0 * np.pi / 5
# Generate a generic azimuthally symmetric field function
Fld = Fields(1.0, 2.0 * k)
E.interpolate(Fld.E)
H.interpolate(Fld.H)

SymmType = 1
# Plot field vectors
Vw = basix.ufl.element('DG', mesh.basix_cell(), 0, shape=(mesh.geometry.dim, ))
W = fem.functionspace(mesh, Vw)
Et = fem.Function(W)
Et.interpolate(E)
Et.name = "ElectricField"
with io.XDMFFile(mesh.comm, "Efield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
H_expr = fem.Expression(ufl.curl(E), W.element.interpolation_points())
Et.interpolate(H)
Et.name = "MagneticField"
with io.XDMFFile(mesh.comm, "Hfield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
# Scalar "phase" function space
Q = fem.functionspace(mesh, ('CG', 2))
vv = fem.Function(Q)

n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
a = 0.0

# Generate integral over quarter sphere surface, using symmetry to get value over full hemisphere surface for z>0
for p in range(26):
theta = np.pi * p / 50
for q in range(101):
phi = np.pi * q / 50
L= []
# Observation vector at position theta, phi (usual sphetrical coords)
rr = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
# U contains constant a_phi unit vector for observation point
U.interpolate(lambda x: (-np.sin(phi) + 0.0 * x[0], np.cos(phi) + x[1] * 0.0, 0.0 * x[2]))
# This is the scalar phase term
vv.interpolate(lambda x: np.exp(1j * k * (rr[0] * x[0] + rr[1] * x[1] + rr[2] * x[2])))
L.append(vv * ufl.cross(n,E))
# Do integration over sphere surface.  Ev should contain the azimuthal component of total field at point theta,phi
Ev = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(sum(L), U) * ds(1))), op=nMPI.SUM)
# Do the same for theta component (should be close to zero over whole theta/phi range)
U.interpolate(lambda x: (np.cos(theta)*np.cos(phi)+x[0]*0.0, np.cos(theta)*np.sin(phi)+x[1]*0.0, -np.sin(theta)+x[2]*0.0))
Eh = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(sum(L), U) * ds(1))), op=nMPI.SUM)
# Print out results
if mpiRank == 0:
print(theta, phi, np.absolute(Ev), np.absolute(Eh), Ev, Eh)
# Here we see the desired azimuthal symmetry.   Scaling by the phase factor is working correctly.
sys.exit(0)


I now see the azimuthal symmetry and the solution behavior is correct. The conclusion is there seems to be a problem with my treatment of the symmetry conditions of the integrand (the Mcx and Mcz factors in the first code). Am doing the symmetry matrix multiplication incorrectly? The numerical result for the symmetry-based solution should be exactly the same as (or very close to) the full sphere case.

For the interest of the community, I managed to get the symmetry problem to work. It turns out that I should have been more careful in the application of the rectangular-to-spherical vector systems and the symmetry conditions. Integration needs to be carried out over each of the components with cartesian unit vectors before applying spherical transformation - a bit of a silly error on my part. Here is an example illustrative code for reader reference, if this problem should crop up for someone else. The component-by-component integration is a bit ugly, but I need to think about how to do this better. Cheers!

import sys
import numpy as np
from dolfinx import default_scalar_type
from mpi4py import MPI as nMPI
from dolfinx import fem
import ufl
import basix.ufl
from dolfinx.io import gmshio
from dolfinx import io
import gmsh

# Prototypical field for point radiator
class Fields:
def __init__(self, A, k0):
self.A = A
self.k0 = k0
def E(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Er = self.A * st * np.array([ct * cp, ct * sp, -st]) * np.exp(-1j * self.k0 * r) / r
return Er
def H(self, x):
r = np.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
ct = x[2] / r
st = np.sqrt(1.0 - ct * ct)
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
cp = x[0] / rr
sp = x[1] / rr
Hr = self.A * st * np.array([-sp, cp, 0.0 * x[2]]) * np.exp(-1j * self.k0 * r) / r
return Hr

# Spherical geometry
eps = 1.0e-4
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
if mpiRank == 0:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)

gmsh.model.occ.addBox(0, -10, 0, 10, 20, 10, 3)
gmsh.model.occ.intersect([(3,1)],[(3,3)], 4, removeObject=True, removeTool=True)
gmsh.model.occ.cut([(3,4)], [(3,2)], 5, removeObject=True, removeTool=True)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(0)
gmsh.model.mesh.setSize(pt, 0.8)
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
for bx in pt:
CoM = gmsh.model.occ.getCenterOfMass(bx[0], bx[1])
if CoM[2] > 1.0+eps and CoM[0] > 2.0:
else:
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 36)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber("Mesh.Tetrahedra", 0)
gmsh.option.setNumber("Mesh.Smoothing", 10)

gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Gmsh")

#    gmsh.fltk.run()
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, 0, gdim=3)

if mpiRank == 0:
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=3)
V = fem.functionspace(mesh, elem)
E = fem.Function(V)
H = fem.Function(V)
S = fem.Function(V)
T = fem.Function(V)
U = fem.Function(V)
k = 2.0 * np.pi / 5
# Generate a generic azimuthally symmetric field function
Fld = Fields(1.0, 2.0 * k)
E.interpolate(Fld.E)
H.interpolate(Fld.H)

SymmType = 1
# Plot field vectors - see what they look like
Vw = basix.ufl.element('DG', mesh.basix_cell(), 0, shape=(mesh.geometry.dim, ))
W = fem.functionspace(mesh, Vw)
Et = fem.Function(W)
Et.interpolate(E)
Et.name = "ElectricField"
with io.XDMFFile(mesh.comm, "Efield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
H_expr = fem.Expression(ufl.curl(E), W.element.interpolation_points())
Et.interpolate(H)
Et.name = "MagneticField"
with io.XDMFFile(mesh.comm, "Hfield_{0}_{1}.xdmf".format(0, SymmType), "w") as xx:
xx.write_mesh(mesh)
xx.write_function(Et)
# Scalar "phase" function space
Q = fem.functionspace(mesh, ('CG', 2))
vv = fem.Function(Q)

n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
a = 0.0

# Generate integral over quarter sphere surface, using symmetry to get value over full hemisphere surface for z>0
Ex = np.zeros(4, dtype=np.complex128)
Ey = np.zeros(4, dtype=np.complex128)
Ez = np.zeros(4, dtype=np.complex128)
Hx = np.zeros(4, dtype=np.complex128)
Hy = np.zeros(4, dtype=np.complex128)
Hz = np.zeros(4, dtype=np.complex128)
# The field symmetry conditions for E and H on the PEC ground plane and PMC symmetry wall
PMCsymmE = np.array([[-1, 0, 0],[0, 1, 0],[0, 0, 1]])
PECsymmE = np.array([[-1, 0, 0],[0, -1 ,0],[0 ,0 ,1]])
PMCsymmH = np.array([[1, 0, 0],[0, -1, 0],[0, 0, -1]])
PECsymmH = np.array([[1, 0, 0],[0, 1 ,0],[0 ,0 ,-1]])

#Loop over theta-phi angles
for p in range(26):
theta = np.pi * p / 50
for q in range(21):
phi = np.pi * q / 10

# Transformation from rectangular to spherical coordinates
Rot = np.array([[np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)],\
[np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)],\
[-np.sin(phi), np.cos(phi), 0.0]])
# Observation vector at position theta, phi (usual sphetrical coords)
rr = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
# S, T, U contains constant unit vectors along x-y-z axes for integration
S.interpolate(lambda x: (1.0 + 0.0 * x[0], x[1] * 0.0, 0.0 * x[2]))
T.interpolate(lambda x: (0.0 * x[0], 1.0 + x[1] * 0.0, 0.0 * x[2]))
U.interpolate(lambda x: (0.0 * x[0], 0.0 * x[1], 1.0 + 0.0 * x[2]))
# This is the scalar phase term
vv.interpolate(lambda x: np.exp(1j * k * (rr[0] * x[0] + rr[1] * x[1] + rr[2] * x[2])))
# Add to integrand original quarter sphere for x>0, z>0 for all three components of current sources
Hx[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), S) * ds(1))), op = nMPI.SUM)
Hy[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), T) * ds(1))), op = nMPI.SUM)
Hz[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), U) * ds(1))), op = nMPI.SUM)
Ex[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), S) * ds(1))), op = nMPI.SUM)
Ey[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), T) * ds(1))), op = nMPI.SUM)
Ez[0] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), U) * ds(1))), op = nMPI.SUM)
# The other quarter sphere using symmetry.  Phase factor for x < 0, z > 0. PMC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (-rr[0] * x[0] + rr[1] * x[1] + rr[2] * x[2])))
Hx[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), S) * ds(1))), op = nMPI.SUM)
Hy[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), T) * ds(1))), op = nMPI.SUM)
Hz[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), U) * ds(1))), op = nMPI.SUM)
Ex[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), S) * ds(1))), op = nMPI.SUM)
Ey[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), T) * ds(1))), op = nMPI.SUM)
Ez[1] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), U) * ds(1))), op = nMPI.SUM)
# The third quarter sphere using symmetry.  Phase factor for x > 0, z < 0. PEC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (rr[0] * x[0] + rr[1] * x[1] - rr[2] * x[2])))
Hx[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), S) * ds(1))), op = nMPI.SUM)
Hy[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), T) * ds(1))), op = nMPI.SUM)
Hz[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), U) * ds(1))), op = nMPI.SUM)
Ex[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), S) * ds(1))), op = nMPI.SUM)
Ey[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), T) * ds(1))), op = nMPI.SUM)
Ez[2] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), U) * ds(1))), op = nMPI.SUM)
# The last quarter sphere using symmetry.  Phase factor for x < 0, z < 0. PEC and PMC symmetry
vv.interpolate(lambda x: np.exp(1j * k * (-rr[0] * x[0] + rr[1] * x[1] - rr[2] * x[2])))
Hx[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), S) * ds(1))), op = nMPI.SUM)
Hy[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), T) * ds(1))), op = nMPI.SUM)
Hz[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,E), U) * ds(1))), op = nMPI.SUM)
Ex[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), S) * ds(1))), op = nMPI.SUM)
Ey[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), T) * ds(1))), op = nMPI.SUM)
Ez[3] = mesh.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(vv * ufl.cross(n,ufl.curl(E)), U) * ds(1))), op = nMPI.SUM)
#
# Do integration of the four sphere slices over quarter sphere surface.  Ev should contain the azimuthal component of total field at point theta,phi