Hello
I am new to fenics, although I have succefully implemented this in sparselizard recently.
I am trying to resolve the eigenmodes of this structure
I have succesully converted the gmsh generated mesh to an .xml using this function
subprocess.run([“dolfin-convert”,“MSL_XY.msh”, “MSL_XY.xml”])
How do I set expressions (formerly known as parameters), such as epsilon_r and mu_r, with respect to their physical region marking?
Thank you
Hi, welcome to the FEniCS community!
It is pretty hard to guess what you would like to achieve without a MWE. I recommend you to read https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21.
Assuming that you import the tags:
mesh = Mesh("MSL_XY.xml")
subdomains = MeshFunction("size_t", mesh, "MSL_XY_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "MSL_XY_facet_region.xml")
and redefine the measures as
dx = Measure('dx', domain=mesh, subdomain_data= subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
Then you can use this to define your expression for each subdomain, i.e.,
def eps1(u):
return 0.5*(grad(u) + grad(u).T)
def eps2(u):
return 0.6*(grad(u) + grad(u).T)
I1 = inner(eps1(u),eps1(v))*dx(1)
I2 = inner(eps2(u),eps2(v))*dx(2)
where dx(1) and dx(2) are the volumetric measures for subdomains tagged as 1 and 2, respectively. Similarly for ds and dS.
You could also define a custom expression:
class Mu(UserExpression):
"""
This class defines a variable mu that takes different values
on the unitsquare domain (0,1)^2 at x<0.5 and x>0.5. The mesh is
assumed to be conforming in x=0.5
"""
def __init__(self, p1, p2, **kwargs):
super().__init__(**kwargs)
self.p1 = p1
self.p2 = p2
def eval(self, value, x):
if x[0] < 0.5. + DOLFIN_EPS:
value[0] = self.p1(x)
else:
value[0] = self.p2(x)
def value_shape(self):
return ()
p1 = Expression('x[1]', degree=1)
p2 = Expression('x[0]+x[1]', degree=1)
mu = Mu(p1, p2, degree=2)
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
Now I have a better view of what you are aiming for. Since you were converting to xml I thought that you were using the Legacy FEniCS. Since your are working in FEniCSx, you can use the GMSHIO to import your mesh tags as follows:
from dolfinx.io import gmshio
from dolfinx import fem, default_scalar_type
import ufl
msh, cell_tags, facet_tags = gmshio.read_from_msh( MSL_XY.msh, MPI.COMM_WORLD, gdim=2)
Please read the post from @dokken on how to save .msh files.
Then, you define your parameters as follows (I will assume same tags as you)
cells_on1 = cell_tags.find(1)
cells_on2 = cell_tags.find(2)
cells_on3 = cell_tags.find(3)
cells_on4 = cell_tags.find(4)
cells_on5 = cell_tags.find(5)
cells_on125= np.concatenate((cells_on1, cells_on2, cells_on5))
DG0 = fem.FunctionSpace(msh, ("DG", 0))
e_r = fem.Function(DG0)
e_r.x.array[cells_on3] = np.full_like(cells_on3, 1.0, dtype=default_scalar_type)
e_r.x.array[cells_on4] = np.full_like(cells_on4, 4.5, dtype=default_scalar_type)
In this way your parameters are defined per cell on the mesh. Since m_r is defined everywhere with single value, then it is enough to take m_r=1.0. If not, follow the above.
Since you want Et=0 and Ez=0 on tags 1,2,5, and they are trial or test functions, you can restrict the integrals by redefining the measures:
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tags)
a_tt = ((1/mr)*ufl.inner(ufl.curl(Et), ufl.curl(Et_tf)) - (k0**2) * e_r * ufl.inner(Et, Et_tf)) * ufl.dx((3,4))
b_tt = (1/mr)*ufl.inner(Et, Et_tf) * ufl.dx((3,4))
b_tz = (1/mr)*ufl.inner(Et, ufl.grad(Ez_tf)) * ufl.dx((3,4))
b_zt = ufl.inner(ufl.grad(Ez), Et_tf) * ufl.dx((3,4))
b_zz = ((1/mr)*ufl.inner(ufl.grad(Ez), ufl.grad(Ez_tf)) - (k0**2) * e_r * ufl.inner(Ez, Ez_tf)) * ufl.dx((3,4))
A fantastic tutorial of what I did can be found in https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html
I hope this is what you’re looking for.
Cheers.
Thank you, @Jesus_Vellojin, for the elaborate answer! I will look into the tutorial, and hopefully will have something up and running, soon enough.
Just one question: Why did you choose dicontinuous Lagrange for the e_r specifically?
Thank you
Because DG-0 is a single constant per cell, just like the tags that you have used for marking domains.