Finally I found the solution according to the dolfinx demos. I’ll also attach the code. My code is not perfect and hope anyone can refine it.
points_list = []
unit_coeff = 0.1
l = 11.17113576837
w = 11.43
initial_points = [[-l, w], [-l, -w], [-6.599113576837424, -w]]
initial_points = (np.array(initial_points)*unit_coeff).tolist()
[points_list.append(x + [0]) for x in initial_points]
for phi in [2*np.pi/3, 4*np.pi/3]:
rot_matrix = np.array([[np.cos(phi), -np.sin(phi)], [np.sin(phi), np.cos(phi)]])
for i, coord in enumerate(initial_points):
point = np.dot(rot_matrix, coord)
point_extend = point.tolist() + [0]
points_list.append(point_extend)
geom = pygmsh.built_in.Geometry()
# Draw a cross.
lcar = 0.1
points = [[-0.35, 0, 0], [0.35, 0, 0], [0, 0, 0]]
circle_points = [geom.add_point(p, lcar=lcar) for p in points]
circle_edge1 = geom.add_circle_arc(circle_points[0], circle_points[-1], circle_points[1])
circle_edge2 = geom.add_circle_arc(circle_points[1], circle_points[-1], circle_points[0])
edges = [circle_edge1, circle_edge2]
line_loop = [geom.add_line_loop(edges)]
circle = geom.add_plane_surface(line_loop[0])
poly = geom.add_polygon(points_list, lcar=0.1, make_surface=False)
domain = geom.add_plane_surface(poly.line_loop, line_loop)
geom.add_physical(poly.lines[0], 'PortOne')
geom.add_physical(poly.lines[3], 'PortTwo')
geom.add_physical(poly.lines[6], 'PortThree')
geom.add_physical([poly.lines[i] for i in [1, 2, 4, 5, 7, 8]], 'PEC')
geom.add_physical(domain, 'OutDomain')
geom.add_physical(circle, 'InnerDomain')
geo_file = 'test.geo'
with open(geo_file, 'w') as f:
f.write(geom.get_code())
mesh = pygmsh.generate_mesh(geom, dim=2)
points, cells, cell_data, field_data = mesh.points, mesh.cells, mesh.cell_data, mesh.field_data
print(cell_data)
meshio.xdmf.write('tag_triangle.xdmf', meshio.Mesh(
points=points,
cells={'triangle': cells[1].data},
cell_data={
'triangle': [cell_data['gmsh:physical'][1]]
},
field_data=field_data
))
line_tags = cell_data['gmsh:physical'][0]
meshio.xdmf.write("tag_all.xdmf", meshio.Mesh(
points=points,
cells={"line": cells[0].data},
cell_data={"line": [line_tags]}
))
with XDMFFile(MPI.comm_world, "tag_triangle.xdmf") as xdmf_infile:
mesh = xdmf_infile.read_mesh(cpp.mesh.GhostMode.none)
mvc_subdomain = xdmf_infile.read_mvc_size_t(mesh)
# tag_info = xdmf_infile.read_information_int()
with XDMFFile(MPI.comm_world, "tag_all.xdmf") as xdmf_infile:
mvc_boundaries = xdmf_infile.read_mvc_size_t(mesh)
mf_triangle = cpp.mesh.MeshFunctionSizet(mesh, mvc_subdomain, 0)
mf_line = cpp.mesh.MeshFunctionSizet(mesh, mvc_boundaries, 0)
gm = 2.8e-3
f0 = gm * 200
fm = gm * 1317
f = 10
c = c * 100
fa = gm * 135
c1 = 1 + fm * (f0 + 1j * fa) / ((f0 + 1j * fa) ** 2 - f ** 2)
c2 = -f * fm / ((f0 + 1j * fa) ** 2 - f ** 2)
k0 = 2 * np.pi * f / c * 1e9
gm = (k0 ** 2 - (np.pi / (2 * w * 0.1)) ** 2) ** 0.5
V = FunctionSpace(mesh, ('P', 2))
u = TrialFunction(V)
v = TestFunction(V)
S = FunctionSpace(mesh, ('DG', 0))
mu = Function(S)
k = Function(S)
e_r = Function(S)
g = Function(V)
g.interpolate(lambda x: np.cos(np.pi * x[1] / (2 * w * 0.1)))
def m(x):
inner_condition = x[0] ** 2 + x[1] ** 2 < 0.35 ** 2
inner_value = 1 + (f0 + 1j * fa) * fm / ((f0 + 1j * fa) ** 2 - f ** 2)
outer_value = 1
return inner_condition * inner_value + outer_value * (~inner_condition)
def kappa(x):
inner_condition = x[0] ** 2 + x[1] ** 2 < 0.35 ** 2
inner_value = -f * fm / ((f0 + 1j * fa) ** 2 - f ** 2)
outer_value = 0
return inner_condition * inner_value + outer_value * (~inner_condition)
def epsilon(x):
inner_condition = x[0] ** 2 + x[1] ** 2 < 0.35 ** 2
inner_value = 11.7
outer_value = 1
return inner_condition * inner_value + outer_value * (~inner_condition)
mu.interpolate(lambda x: m(x))
k.interpolate(lambda x: kappa(x))
e_r.interpolate(lambda x: epsilon(x))
u0 = Function(V)
with u0.vector.localForm() as bc_local:
bc_local.set(0.0)
bndry_facets = np.where(mf_line.values == 3)[0]
bdofs = locate_dofs_topological(V, 1, bndry_facets)
bc = DirichletBC(u0, bdofs)
ds = ds(subdomain_data=mf_line)
a = mu / (mu ** 2 - k ** 2) * inner(grad(u), grad(v)) * dx - \
1j * k / (mu ** 2 - k ** 2) * (inner(u.dx(0), v.dx(1)) - inner(u.dx(1), v.dx(0))) * dx - \
e_r * k0 ** 2 * inner(u, v) * dx + 1j * gm * inner(u, v) * (ds(1) + ds(2))
L = 2j * gm *inner(g, v) * ds(1)
u = Function(V)
solve(a == L, u, bc)
with XDMFFile(MPI.comm_world, "plane_wave.xdmf",
encoding=XDMFFile.Encoding.HDF5) as file:
file.write(u)