Hello again! I decided to move my questions from the spack installation thread back to this discussion.
The generation of matrix blocks using the Scifem real function spaces looks like the way to go for building these types of matrices. I am now trying to implement this method into a simulation that uses a GMSH generated mesh and I have come up against another roadblock. In the code:
# Dolfin 0.9.0, PETSc 3.24.1
import sys
import numpy as np
from mpi4py import MPI as M
from petsc4py import PETSc as P
import gmsh
import ufl as U
import basix.ufl
from dolfinx.io import gmshio
from dolfinx.fem import petsc
from dolfinx import mesh as msh
from dolfinx import fem
from dolfinx import io
import scifem
import time
eta0 = 377.0 # free space wave impedance
class ModesTE:
def __init__(self, beta, k, a, b, m, n):
self.beta = beta
self.k = k
self.a = a
self.b = b
self.m = m
self.n = n
def eval_H(self, x):
if((self.m != 0) and (self.n != 0)):
B = np.sqrt(self.n * self.m * 4.0 / (self.a * self.b))
elif self.m == 0:
B = np.sqrt(self.n * 2.0 / (self.a * self.b))
else:
B = np.sqrt(self.m * 2.0 / (self.a * self.b))
hx = B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
hy = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
hz = np.full_like(hx, 0.0+0j, dtype=np.complex128)
return(hx, hy, hz)
def eval_E(self, x):
if((self.m != 0) and (self.n != 0)):
B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
elif self.m == 0:
B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
else:
B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
ZTE = eta0 * self.k / self.beta
ey = -ZTE * B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
ex = ZTE * B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
ez = np.full_like(hx, 0.0+0j, dtype=np.complex128)
return(ex, ey, ez)
eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
f = 10.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 6.60
comm = M.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency
if mpiRank == modelRank:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add('Waveguide')
gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(dim=0)
gmsh.model.mesh.setSize(pt, lm)
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 2.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', 30)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber('Mesh.Tetrahedra', 0)
OUT = []
IN = []
PEC = []
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
for bnd in pt:
CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
if np.abs(CoM[2]) < eps:
IN.append(bnd[1])
elif np.abs(CoM[2] - l) < eps:
OUT.append(bnd[1])
else:
PEC.append(bnd[1])
gmsh.model.addPhysicalGroup(3, [1], 1) #WG
gmsh.model.setPhysicalName(3, 1, "Waveguide")
gmsh.model.addPhysicalGroup(2, PEC, 1)
gmsh.model.setPhysicalName(2, 1, "PEC")
gmsh.model.addPhysicalGroup(2, IN, 2)
gmsh.model.setPhysicalName(2, 2, "INPORT")
gmsh.model.addPhysicalGroup(2, OUT, 3)
gmsh.model.setPhysicalName(2, 3, "OUTPORT")
gmsh.model.mesh.generate(3)
# gmsh.fltk.run()
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
gmsh.finalize()
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
#Real function space over boundary tag 2
facetsB = fm.find(2)
Gamma, entity_map1, _, _ = msh.create_submesh(mesh, mesh.topology.dim-1, facetsB)
print(entity_map1)
Z1 = scifem.create_real_functionspace(Gamma, value_shape=(1,))
z1u = U.TrialFunction(Z1)
z1v = U.TestFunction(Z1)
elem = basix.ufl.element('N1curl', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
normal = U.FacetNormal(mesh)
ds = U.Measure("ds", domain=mesh, subdomain_data = fm)
# Dirichlet BCs on boundary tag 1
facets = fm.find(1)
ubc = fem.Function(V)
ubc.interpolate(lambda x: np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
f = ModesTE(kz, k0, a, b, 1, 0)
# Build boundary mode vector
gh = fem.Function(V, dtype=np.complex128)
gh.interpolate(f.eval_H)
t_start = time.perf_counter()
n = U.FacetNormal(mesh)
# Weak form
num_facets_parent = mesh.topology.index_map(mesh.topology.dim-1).size_local + mesh.topology.index_map(mesh.topology.dim-1).num_ghosts
inverse_map = np.full(num_facets_parent, -1, dtype=np.int32)
inverse_map[facetsB] = np.arange(len(facetsB), dtype=np.int32)
L00 = fem.form((U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx)
L10 = fem.form(1j * kz * U.inner(U.dot(U.cross(u, gh), n), z1v) * ds(2), entity_maps={Gamma:inverse_map})
L01 = fem.form(1j * kz * U.inner(U.dot(U.cross(v, gh), n), z1u) * ds(2), entity_maps={Gamma:inverse_map})
A = petsc.create_matrix([[L00, L01], [L10, None]], "mpiaij")
petsc.assemble_matrix(A, [[L00, L01], [L10, None]], bcs=[bc])
A.assemble()
sys.exit(0)
I am generating three sub-matrices. L00 is the usual weak form square matrix. L01 and L10 are two forms that use the scifem real space and should each generate an Nx1 and 1xN matrix, respectively. Unfortunately, there is a shape error in building the real space forms. In the interfacing with the GMSH mesh, I am probably missing something in passing the mesh information to the real space definition. This is the error:
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/WaveguideBC/Rect2.py", line 154, in <module>
L10 = fem.form(1j * kz * U.inner(U.dot(u, n), z1v) * ds(2), entity_maps=[entity_map1])
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/bill/spack/var/spack/environments/Fenics-090-complex/.spack-env/view/lib/python3.12/site-packages/ufl/operators.py", line 213, in inner
return Inner(a, b)
^^^^^^^^^^^
File "/home/bill/spack/var/spack/environments/Fenics-090-complex/.spack-env/view/lib/python3.12/site-packages/ufl/tensoralgebra.py", line 166, in __new__
raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Dot id=125006079896432> and <Argument id=125006079775024>
EDIT 16.55, 2025-11-28: I realized that I posted the wrong code. Updated with the proper entity map lines at the end of the code block.