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:
- 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.
- 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.add("Prototypical field profile sphere")
gmsh.model.occ.addSphere(0, 0, 0, 10.0, 1)
gmsh.model.occ.addSphere(0, 0, 0, 1.0, 2)
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)
rad = []
norad = []
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:
rad.append(bx[1])
else:
norad.append(bx[1])
gmsh.model.addPhysicalGroup(3, [5], 1)
gmsh.model.setPhysicalName(3, 1, "Rad Region")
gmsh.model.addPhysicalGroup(2, rad, 1)
gmsh.model.setPhysicalName(2, 1, "Rad")
gmsh.model.addPhysicalGroup(2, norad, 2)
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])))
# Add to integrand
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])))
# Add to integrand
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])))
# Add to integrand
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.add("Prototypical field profile sphere")
gmsh.model.occ.addSphere(0, 0, 0, 10.0, 1)
gmsh.model.occ.addSphere(0, 0, 0, 1.0, 2)
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)
rad = []
norad = []
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:
rad.append(bx[1])
else:
norad.append(bx[1])
gmsh.model.addPhysicalGroup(3, [5], 1)
gmsh.model.setPhysicalName(3, 1, "Rad Region")
gmsh.model.addPhysicalGroup(2, rad, 1)
gmsh.model.setPhysicalName(2, 1, "NoRad")
gmsh.model.addPhysicalGroup(2, norad, 2)
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])))
# Add to integrand
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.