Hi dokken, it looks like the mesh generated by gmsh caused the problem. However I don’t know why. Here is my mesh and some relevant code. Sorry to trouble you again. Since the mesh cannot be uploaded due to my slow internet, I upload the geo file insteady, you can generate the mesh from strip.geo.
strip.geo file:
//+
d_w = DefineNumber[ 0.0127, Name "Parameters/d_w" ];
//+
d_h = DefineNumber[ 1.27e-3, Name "Parameters/d_h" ];
//+
s_w = DefineNumber[ 1.27e-3, Name "Parameters/s_w" ];
//+
s_h = DefineNumber[ 1.27e-4, Name "Parameters/s_h" ];
//+
Point(1) = {0, 0, 0, 0.0005};
//+
Point(2) = {d_w, 0, 0, 0.0005};
//+
Point(3) = {d_w, d_h, 0, 0.0005};
//+
Point(4) = {d_w, d_w, 0, 0.0005};
//+
Point(5) = {0, d_w, 0, 0.0005};
//+
Point(6) = {0, d_h, 0, 0.0005};
//+
Point(7) = {d_w/2-s_w/2, d_h, 0, 0.00005};
//+
Point(8) = {d_w/2+s_w/2, d_h, 0, 0.00005};
//+
Point(9) = {d_w/2+s_w/2, d_h+s_h, 0, 0.00005};
//+
Point(10) = {d_w/2-s_w/2, d_h+s_h, 0, 0.00005};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {6, 1};
//+
Line(7) = {6, 7};
//+
Line(8) = {7, 8};
//+
Line(9) = {8, 3};
//+
Line(10) = {8, 9};
//+
Line(11) = {9, 10};
//+
Line(12) = {10, 7};
//+
Curve Loop(1) = {5, 7, -12, -11, -10, 9, 3, 4};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {7, 8, 9, -2, -1, -6};
//+
Plane Surface(2) = {2};
//+
Physical Curve("PEC") = {5, 6, 1, 2, 3, 4, 11, 12, 8, 10};
//+
Physical Surface("Air") = {1};
//+
Physical Surface("Dielectric") = {2};
mesh->xdmf:
import pygmsh
import meshio
import numpy as np
filename = 'strip.msh'
mesh = meshio.read(
filename, # string, os.PathLike, or a buffer/open file
)
points, cells, cell_data, field_data = mesh.points, mesh.cells, mesh.cell_data, mesh.field_data
triangle_data = None
line_data = None
tri_cell_idx = []
line_cell_idx = []
for i, cell in enumerate(cells):
if cell.type == 'triangle':
print(i, 'triangle')
tri_cell_idx.append(i)
if triangle_data is None:
triangle_data = cell.data
else:
triangle_data = np.vstack([triangle_data, cell.data])
else:
print(i, 'line')
line_cell_idx.append(i)
if line_data is None:
line_data = cell.data
else:
line_data = np.vstack([line_data, cell.data])
meshio.xdmf.write('tag_triangle.xdmf', meshio.Mesh(
points=points,
cells={'triangle': triangle_data},
cell_data={
'triangle': [np.concatenate([cell_data['gmsh:physical'][i] for i in tri_cell_idx])]
},
field_data=field_data
))
meshio.xdmf.write("tag_all.xdmf", meshio.Mesh(
points=points,
cells={"line": line_data},
cell_data={
"line": [np.concatenate([cell_data['gmsh:physical'][i] for i in line_cell_idx])]
}
))
xdmf->dolfinx:
import numpy as np
import pygmsh
import meshio
from dolfinx.io import XDMFFile
from dolfinx import (MPI, FacetNormal, Function, FunctionSpace,
has_petsc_complex, DirichletBC, cpp)
from dolfinx.fem.assemble import assemble_scalar, assemble_matrix
from dolfinx.fem import locate_dofs_topological
from ufl import (TestFunctions, TrialFunctions, grad, inner, curl, split,
Measure, FiniteElement, dx)
import dolfinx
from scipy.constants import c
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)
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if mf_triangle.values[i] == 2:
e_r.vector.setValueLocal(i, 8.875)
else:
e_r.vector.setValueLocal(i, 1)
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = 10e9
k0_sqr = (2*np.pi*f / c)**2
electric_wall = Function(W)
with electric_wall.vector.localForm() as bc_local:
bc_local.set(0)
bndry_facets = np.where(mf_line.values == 1)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc= DirichletBC(electric_wall, bdofs)
a = inner(curl(u), curl(v))*dx - k0_sqr*e_r*inner(u, v)*dx
b = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
A = assemble_matrix(a, [bc])
B = assemble_matrix(b, [bc])
A.assemble()
C = A.convert('dense')
