Here’s my code and there is no problem to solve it. I tested assemble_scalar(inner(dot(grad(u)(’+’), n(’+’)), v(’+’)) * dS(7)) and it can give a value. While I calculated assemble_scalar(inner(dot(grad(u)(’+’), n(’+’)), v(’+’)) * dS(7) +
Constant(mesh, 0)*dx(domain=mesh, subdomain_data=mf_triangle)), an type error happened. If I executed assemble_scalar(inner(dot(grad(u)(’-’), n(’-’)), v(’-’)) * dS(7), the python kernel crash.
mesh file:
import pygmsh
import meshio
import numpy as np
points_list = []
unit_coeff = 0.1
#l = 11.17113576837
l = 11.17113576837
w = 11.43
shift_len = l*0.1
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)
points_list = [[x+shift_len, y, z] for x, y, z in points_list]
geom = pygmsh.built_in.Geometry()
# Draw a cross.
lcar = 0.08
points = [[-0.35+shift_len, 0, 0], [0.35+shift_len, 0, 0], [0+shift_len, 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')
geom.add_physical(edges, 'Interface')
geo_file = 'test.geo'
with open(geo_file, 'w') as f:
f.write(geom.get_code())
mesh = pygmsh.generate_mesh(geom, prune_z_0=True, verbose=True, geo_filename="test.geo")
points, cells, cell_data, field_data = mesh.points, mesh.cells, mesh.cell_data, mesh.field_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]}
))
print(cells, cell_data, field_data)
main file
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 = 0.0176
gm = 2.8e-3
f0 = gm * 200
fm = gm * 1317
f = 9
omega = 2*np.pi*f
#omega = 10
c0 = c * 100
fa = gm * 135 / 2
w = 11.43
c1 = 1 + fm * (f0 + 1j * fa) / ((f0 + 1j * fa) ** 2 - f ** 2)
c2 = -f * fm / ((f0 + 1j * fa) ** 2 - f ** 2)
k0 = omega / c0 * 1e9
wavenumber = (-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)
W = FunctionSpace(mesh, ('P', 4))
g = Function(W)
g.interpolate(lambda x: np.cos(np.pi * x[1] / (2 * w * 0.1)))
m_i = c1
m_o = 1
k_i = c2
k_o = 0
shift_len = 1.117113576837
def m(x):
inner_condition = (x[0]-shift_len) ** 2 + x[1] ** 2 < 0.35 ** 2
inner_value = c1
outer_value = 1
return inner_condition * inner_value + outer_value * (~inner_condition)
def kappa(x):
inner_condition = (x[0]-shift_len) ** 2 + x[1] ** 2 < 0.35 ** 2
inner_value = c2
outer_value = 0
return inner_condition * inner_value + outer_value * (~inner_condition)
def epsilon(x):
inner_condition = (x[0]-shift_len) ** 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 == 4)[0]
bdofs = locate_dofs_topological(V, 1, bndry_facets)
bc = DirichletBC(u0, bdofs)
ds = Measure('ds', subdomain_data=mf_line)
dS = Measure('dS', subdomain_data=mf_line)
dx = Measure('dx', subdomain_data=mf_triangle)
n = FacetNormal(mesh)
t = perp(n)
a = mu / (mu ** 2 - k ** 2) * inner(grad(u), grad(v)) * dx - \
(m_i / (m_i ** 2 - k_i ** 2) - m_o / (m_o ** 2 - k_o ** 2)) * \
inner(dot(grad(u)('+'), n('+')), v('+')) * dS(7) -\
1j * k / (mu ** 2 - k ** 2) * (inner(u.dx(0), v.dx(1)) - inner(u.dx(1), v.dx(0))) * dx -\
1j * k_i / (m_i ** 2 - k_i ** 2) * inner(dot(grad(u)('+'), t('+')), v('+')) * dS(7) -\
e_r * k0 ** 2 * inner(u, v) * dx + wavenumber * inner(u, v) * (ds(1) + ds(2) + ds(3))
L = 2 * wavenumber *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)
