Hey @Jesus_Vellojin, thank you for the swift answer.
Sadly, I do not have a MWE, as I can’t find this basic tool for me to work with. I did go over the demos, available in the Docker image, as well as the 2011 version of Automated Solution of Differential Equations by the Finite Element Method fenics guide found on the web and the documented example here:
https://docs.fenicsproject.org/dolfinx/v0.7.2/python/demos/demo_half_loaded_waveguide.html
I’ll add what I can, here, though.
When I generate the mesh for the structure above (a microstrip line), I assign physical tags for boundaries and faces. Here is the mesh generation script:
# Import modules:
import gmsh
import sys
import os
# Initialize gmsh:
gmsh.initialize()
port_h_fact = 5.5
port_w_fact = 6
unit_conv = 1e-3
Line_W = 1.875*unit_conv
Sub_W = Line_W*(1.0 + port_w_fact)
Sub_Thick = 1*unit_conv;
Cu_Thick = 0.1*unit_conv;
meshSize = 0.5*unit_conv
tag_gnd = 1
tag_line = 2
tag_air = 3
tag_fr4 = 4
tag_bounds = 5
# cube points:
# Substrate rectangle
# XxZ = -Y
point1 = gmsh.model.geo.add_point(-Sub_W*0.5, 0, 0, meshSize)
point2 = gmsh.model.geo.add_point(Sub_W*0.5, 0, 0, meshSize)
point3 = gmsh.model.geo.add_point(Sub_W*0.5, Sub_Thick, 0, meshSize)
point4 = gmsh.model.geo.add_point(-Sub_W*0.5, Sub_Thick, 0, meshSize)
# line rectangle
point5 = gmsh.model.geo.add_point(-Line_W*0.5, Sub_Thick, 0.0, meshSize)
point6 = gmsh.model.geo.add_point(Line_W*0.5, Sub_Thick, 0.0, meshSize)
point7 = gmsh.model.geo.add_point(Line_W*0.5, Sub_Thick + Cu_Thick, 0.0, meshSize)
point8 = gmsh.model.geo.add_point(-Line_W*0.5, Sub_Thick + Cu_Thick, 0.0, meshSize)
# additional air edges
point9 = gmsh.model.geo.add_point(-Sub_W*0.5, Sub_Thick*(1.0 + port_h_fact), 0.0, meshSize)
point10 = gmsh.model.geo.add_point(Sub_W*0.5, Sub_Thick*(1.0 + port_h_fact), 0.0, meshSize)
# Edges of substrate
l1 = gmsh.model.geo.add_line(point1, point2,1)
l2 = gmsh.model.geo.add_line(point2, point3,2)
l3 = gmsh.model.geo.add_line(point3, point6,3)
l4 = gmsh.model.geo.add_line(point6, point5,4)
l5 = gmsh.model.geo.add_line(point5, point4,5)
l6 = gmsh.model.geo.add_line(point4, point1,6)
# # Remaining lines of microstrip line
l7 = gmsh.model.geo.add_line(point6, point7,7)
l8 = gmsh.model.geo.add_line(point7, point8,8)
l9 = gmsh.model.geo.add_line(point8, point5,9)
#
# # Remaining lines for airs
l10 = gmsh.model.geo.add_line(point3, point10,10)
l11 = gmsh.model.geo.add_line(point10, point9,11)
l12 = gmsh.model.geo.add_line(point9, point4,12)
# faces of cube:
curve_sub = gmsh.model.geo.add_curve_loop([1,2,3,4,5,6])
curve_line = gmsh.model.geo.add_curve_loop([-4,7,8,9])
curve_air = gmsh.model.geo.add_curve_loop([-5,-9,-8,-7,-3,10,11,12])
# surfaces of cube:
surf_sub = gmsh.model.geo.add_plane_surface([curve_sub])
surf_air = gmsh.model.geo.add_plane_surface([curve_air])
# Create the relevant Gmsh data structures
# from Gmsh model.
gmsh.model.geo.synchronize()
# Add metals to physical group
gmsh.model.addPhysicalGroup(1, [-1], tag_gnd)
gmsh.model.addPhysicalGroup(1,[-2,-10,-11,-12,-6],tag_bounds)
gmsh.model.addPhysicalGroup(1, [-4,7,8,9], tag_line)
# Add metals to physical group
gmsh.model.addPhysicalGroup(2, [surf_air], tag_air)
gmsh.model.addPhysicalGroup(2, [surf_sub], tag_fr4)
gmsh.write("MSL_XY.geo_unrolled")
os.rename("MSL_XY.geo_unrolled","MSL_XY.geo")
# gmsh.write("MSL.geo")
# gmsh.open("MSL.geo")
# Generate mesh:
gmsh.model.mesh.generate()
gmsh.write("MSL_XY.msh2")
os.rename("MSL_XY.msh2","MSL_XY.msh")
# Creates graphical user interface
if 'close' not in sys.argv:
gmsh.fltk.run()
# It finalize the Gmsh API
gmsh.finalize()
Now I am trying to run the eigenmode analysis in the documentation
import sys
from mpi4py import MPI
import numpy as np
from analytical_modes import verify_mode
import subprocess
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type, fem, io, plot
from dolfinx.fem.petsc import assemble_matrix
try:
import pyvista
have_pyvista = True
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
try:
from slepc4py import SLEPc
except ModuleNotFoundError:
print("slepc4py is required for this demo")
sys.exit(0)
vector_basisFunction_order = 2
scalar_basisFunction_order = 2
I am doing the conversion in a pretty old fashion I found in another thread. This can probably be done differently, as well.
# Convert mesh
subprocess.run(["dolfin-convert","MSL_XY.msh", "MSL_XY.xml"])
mesh = Mesh("MSL_XY.xml")
And from here I am mostly following the documentation example
# Scalar and vector function spaces
FS_V = element("Nedelec 1st kind H(curl)", msh.basix_cell(), vector_basisFunction_order)
FS_S = element("Lagrange", msh.basix_cell(), scalar_basisFunction_order)
combined_space = fem.functionspace(mesh, mixed_element([FS_V, FS_S]))
Et_tf, Ez_tf = ufl.TrialFunctions(combined_space)
Et, Ez = ufl.TestFunctions(combined_space)
c0 = 299792458.0;
f0 = 2.5e9
k0 = 2.0*np.pi*f0/c0
Here is the part I am struggling with:
<Somehow define a parameter 'e_r' for relative permittivity>
<Somehow define a parameter 'm_r' for relative permeability>
<Somehow define e_r = 1.0 in physical tag 3>
<Somehow define e_r = 4.5 in physical tag 4>
<Somehow define m_r = 1.0 everywhere>
<Somehow define Et=0 and Ez=0 on metalic boundaries, tags 1,2 and 5>
After that the example is pretty clear
a_tt = ((1/mr)*ufl.inner(ufl.curl(Et), ufl.curl(Et_tf)) - (k0**2) * e_r * ufl.inner(Et, Et_tf)) * ufl.dx
b_tt = (1/mr)*ufl.inner(Et, Et_tf) * ufl.dx
b_tz = (1/mr)*ufl.inner(Et, ufl.grad(Ez_tf)) * ufl.dx
b_zt = ufl.inner(ufl.grad(Ez), Et_tf) * ufl.dx
b_zz = ((1/mr)*ufl.inner(ufl.grad(Ez), ufl.grad(Ez_tf)) - (k0**2) * e_r * ufl.inner(Ez, Ez_tf)) * ufl.dx
A = fem.form(a_tt)
B = fem.form(b_tt + b_tz + b_zt + b_zz)
# Now we can solve the problem with SLEPc. First of all, we need to
# assemble our and matrices with PETSc in this way:
A = assemble_matrix(a, bcs=[bc])
A.assemble()
B = assemble_matrix(b, bcs=[bc])
B.assemble()
# Now, we need to create the eigenvalue problem in SLEPc. Our problem is a
# linear eigenvalue problem, that in SLEPc is solved with the EPS module.
# We can initialize this solver in the following way:
eps = SLEPc.EPS().create(mesh.comm)
# We can pass to EPS our matrices by using the setOperators routine:
eps.setOperators(A, B)
I probably did a million things wrong here, but there are too many degrees of freedom for me to grasp at one time.
Since I have already succesfully implemented this in a different FEM library, I am diving straight into the water…
Again, appreciate any help.
Cheers