Hello everyone,
I’m implementing an inverse dispersion solver in DOLFINx that computes complex wavenumbers k(\omega) for waves in a periodic elastic 3D structure. Instead of solving the forward eigenvalue problem \omega(k) (which works well), I want to solve the quadratic eigenvalue problem
[A_0 + \Lambda A_1 + \Lambda ^2 A_2]\textbf{U} = 0,
where \Lambda = ik and the matrices A_0, A_1 and A_2 are defined in the strong form as follows:
A_0 = K_0 - \rho \omega ^2\textbf{U},
K_0 = (\lambda + \mu) \nabla (\nabla \cdot \textbf{U}) + \mu \nabla ^2 \textbf{U}
A_1 = [(\lambda + \mu)(\nabla \textbf{U}_z + \textbf{e}_z (\nabla \cdot \textbf{U}) ) + 2\mu \partial_z \textbf{U}]
A_2 = -[\mu \textbf{U} + (\lambda + \mu)\textbf{e}_z \textbf{U}_z]
These equation are based on eq. A.13, page 30 in the paper given in the following link: https://arxiv.org/pdf/2601.18911.
I’m struggling with the correct implementation of periodic boundary conditions and the PEP solver configuration. The geometry is a cylindrical rod with periodicity along the axial (z) direction, imported as a STEP-file from COMSOL and the meshing is done in Gmsh. My forward solver for the \omega(k) problem uses dolfinx_mpc with the Bloch periodicity imposed at the left and right boundaries of the rod
\textbf{u}_{right} = \textbf{u}_{left} e ^{ik (\textbf{r}_{left} - \textbf{r}_{right})}
where the phase factor e^{ikL} is set as a scaling factor for the mpc.create_periodic_constraint_geometrical module and works correctly.
For the inverse problem, I’ve reformulated the weak formulations using the Bloch substitution \nabla \rightarrow \nabla + ike_z as proposed in the paper, where e_z is the unit vector in z-direction, which leads to the quadratic eigenvalue problem with the coefficient matrices given above. I’ve also modified the mpc formulation to set u_{right} = u_{left} e^{ikL} with e^{ikL} being set to 1 to impose the continuity between the both ends as suggested in the section A.2.3 of the paper.
The weak forms I’ve derived are expressed in the code as:
a: standard elasticity (K_0 term) which is composed bya = a_standard1 + a_standard2corresponding to the two terms in K_0 above.a_linear: linear in k term (coefficient matrix A_1)a_quadratic: quadratic in k term (coefficient matrix A_2)m: mass matrix corresponding to the inertia term
When I solve the quadratic eigenvalue problem using SLEPC with the pep-solver object for each frequency and extract Re(k) = Im(\Lambda) to plot the dispersion curves, I get eigenvalues that are not converging to the expected dispersion curves. I’ve tried various solver types (TOAR , CISS ) and problem types (HERMITIAN , GENERAL ). I’ve also tried different sorting methods (smallest magnitude, smallest imaginary magnitude etc) with no result.
Below is a minimum working example that prints out the eigenvalues \Lambda and plots the corresponding dispersion curves for real and imaginary part of k.
# Minimum Working Example - Solving for complex k (quadratic eigenvalue problem)
#--------------Import packages--------------------------
import numpy as np
import matplotlib.pyplot as plt
from petsc4py import PETSc
from mpi4py import MPI
import dolfinx
import gmsh
import dolfinx.mesh as msh
from dolfinx import mesh, fem, plot
from ufl import dx, inner
from basix.ufl import element
from dolfinx import default_scalar_type
from dolfinx.fem import form, functionspace, assemble_scalar, Function
from dolfinx.fem.petsc import assemble_matrix, LinearProblem
from slepc4py import SLEPc
from dolfinx.io import gmsh as gmshio
import pickle
import os
import dolfinx_mpc.utils
import ufl
import pyvista as pv
from dolfinx_mpc import MultiPointConstraint, assemble_matrix
from dolfinx.fem import form
from slepc4py import SLEPc
from petsc4py import PETSc
import cmath
import os
import pickle
import numpy as np
import matplotlib.pyplot as plt
#--------------Defining functions--------------------------------
def import_CAD_to_FEniCSx(input_file = "object_new_broken_cubic_case", output_file = "from_comsol", lc = 0.25):
"""
Docstring for import_CAD_to_FEniCSx --- This function covert a STEP CAD file into a216 readable FEniCSx mesh with physical groups.
This functions aims at being adapted for you code regarding your need.
You can edit the physcal groups, the mesh size caracteristics, the order of the elements, etc.
Take great care of the listsurf1 variable.
You can adapt the code if you don't need submesh.
:param input_file: str - Name of the input STEP file without the .step extension
:param output_file: str - Name of the output files without the .msh extension
:param lc: float - Typical mesh size
:return: mesh_info: list - List containing the FEniCSx mesh, cell tags and facet tags
submesh_info: list - List containing a submesh and its entity map for a physical group
"""
gmsh.initialize()
comm = MPI.COMM_WORLD
model_rank = 0
model = gmsh.model()
gmsh.model.add(output_file)
gmsh.model.occ.importShapes(f"{input_file}.step")
gmsh.model.occ.synchronize()
# Get all volumes
volumes = gmsh.model.getEntities(dim=3)
print(f'volumes: {volumes}')
volume_tags = [v[1] for v in volumes]
print(f'volume_tags: {volume_tags}')
# Fuse all volumes together to create shared interfaces
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
print(f'volumes after fragment: {volumes}')
gmsh.model.addPhysicalGroup(3, [v[1] for v in volumes], tag=1)
surfaces = gmsh.model.getEntities(dim=2)
print(f'surfaces: {surfaces}')
###### Define surface groups cylinder
listsurf1 = [4]
listsurf2 = [5] # To be editted according to the ge2 elements with jac. < 0ometry -> find the tags thanks to gmsh.fltk.run()
listsurf3= [1, 2, 3, 6]
gmsh.model.addPhysicalGroup(2, listsurf1, tag=1001) # Remember that 1001 will be the tag of integration over this surface in FEniCSx with ds(tag)
gmsh.model.addPhysicalGroup(2, listsurf2, tag=2001)
gmsh.model.addPhysicalGroup(2, listsurf3, tag=3001)
### Mesh characteristics
##Mesh size
gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.option.setNumber("Mesh.HighOrderOptimize", 2) # Better to use it with curved geometries
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc) #2*lc
gmsh.model.mesh.generate(3)
gmsh.write(f"{output_file}.msh") # We are done with Gmsh, we can save the mesh if needed
mesh_data = gmshio.model_to_mesh(
model, comm, model_rank
) # Import the mesh into FEniCSx
final_mesh = mesh_data.mesh # The dolfinx mesh object -> The actual mesh
facet_tags = mesh_data.facet_tags # The tags of the physical surfaces
cell_tags = mesh_data.cell_tags # The tags of the physical volumes
# Comment the following lines to stop visualizing the mesh in Gmsh
gmsh.option.setNumber("Geometry.Surfaces", 1) # To visualize the surfaces ids in Gmsh
gmsh.option.setNumber("Geometry.SurfaceLabels", 1)
#gmsh.fltk.run() # To visualize the mesh in Gmsh
gmsh.finalize()
tdim = final_mesh.topology.dim
fdim = tdim - 1
mesh_info = [final_mesh, cell_tags, facet_tags]
# If you need to create a submesh corresponding to a physical group
submesh, entity_map = msh.create_submesh(final_mesh, fdim, facet_tags.find(1001))[0:2]
entity_maps_mesh = [entity_map]
submesh_info = [submesh, entity_maps_mesh]
return mesh_info, submesh_info
# Function for removing rows and columns to handle rigid body modes
def extract_sub_matrix(A: PETSc.Mat,
V: dolfinx.fem.FunctionSpace,
bcs: list[dolfinx.fem.DirichletBC] = None,
mpc=None) -> PETSc.Mat:
"""
Extract a submatrix corresponding to unconstrained DOFs:
• Dirichlet BCs (bcs != None)
• Periodic/MPC constraints (mpc != None)
"""
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_ghosts = V.dofmap.index_map.num_ghosts * V.dofmap.index_map_bs
if mpc is not None:
# Identify slave DOFs locally
slave_local = np.array(mpc.is_slave, dtype=np.int32) #returns a vector of integers where the ith entry idnicate if a dof is a "slave"
# Create mask: 1 if not slave, 0 if slave
mask = np.ones(num_dofs_local + num_ghosts, dtype=np.int32)
mask[slave_local.astype(bool)] = 0
# Now the unconstrained DOF indices:
idx_set = np.flatnonzero(mask).astype(np.int32)
# Only keep local (non-ghost) rows/columns
idx_set = idx_set[idx_set < num_dofs_local]
elif bcs is not None:
not_dirichlet = np.full(num_dofs_local + num_ghosts, 1, dtype=np.int32)
for bc in bcs:
not_dirichlet[bc.dof_indices()[0]] = 0
idx_set = np.flatnonzero(not_dirichlet).astype(np.int32)
idx_set = idx_set[idx_set < num_dofs_local]
else:
raise ValueError("Either bcs or mpc must be provided.")
# Build PETSc index sets
isx = PETSc.IS(A.getComm()).createGeneral(idx_set)
lgm_row = A.getLGMap()[0].applyIS(isx)
lgm_col = A.getLGMap()[1].applyIS(isx)
# Extract the submatrix
A_sub = A.createSubMatrix(lgm_row, lgm_col)
A_sub.assemble()
return A_sub, (idx_set, idx_set)
# Functions to handle periodicity in the system
def periodic_indicator(x, direction='z', toli=1e-8):
"""Returns True for points on the periodic boundary (max coordinate side)"""
coord_map = {'x': 0, 'y': 1, 'z': 2}
idx = coord_map[direction]
x_max = mesh.geometry.x[:, idx].max()
return np.isclose(x[idx], x_max, atol=toli)
def periodic_relation(x, direction='z', a_basis=None, toli=1e-8):
"""Maps points from max boundary to min boundary"""
out = x.copy()
coord_map = {'x': 0, 'y': 1, 'z': 2}
idx = coord_map[direction]
out[idx] = x[idx] - a_basis[idx]
return out
#--------------Create mesh from the step file--------------------
#SAVE_PATH = '/home/nadine/Documents/Msc Thesis/Code_conda/results'
step_filename, mesh_name = 'Test_periodic_rod_geo_r025', 'mesh_rod'
lc = 10 # characteristic mesh element size
mesh_info, submesh_info = import_CAD_to_FEniCSx(step_filename, mesh_name, lc)
mesh, cell_tags, facet_tags = mesh_info
dim = mesh.geometry.dim
mesh.geometry.x[:] *= 0.001 # Need to downscale the geometry to original size 1 m
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)
#-------------------Geometry properties---------------------------
z_coords = mesh.geometry.x[:, 2]
L_f = z_coords.max() - z_coords.min()
# ---------------Define space function and element----------------------------
deg = 2
Ve = element("Lagrange", mesh.basix_cell(), deg, shape=(dim,))
V = functionspace(mesh, Ve)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#-------------Material parameters-----------------------------------------------------------
material_tags = np.unique(cell_tags.values) if cell_tags is not None else np.array([])
print(f"Material tags found in mesh: {material_tags}")
material_properties = {}
if len(material_tags) == 1:
# Single material case
print("Single material detected")
# Use the actual tag from the mesh
tag = material_tags[0]
material_properties[tag] = {
"name": "main_structure",
"E": 2.10e14,
"rho": 7850,
"nu": 0.3
}
print("Material configuration:")
for tag, props in material_properties.items():
print(f" Tag {tag} ({props['name']}): E={props['E']:.2e}, rho={props['rho']}, nu={props['nu']}")
Q = functionspace(mesh, ("DG", 0))
E = Function(Q)
rho = Function(Q)
nu = Function(Q)
mu = Function(Q)
lambda_ = Function(Q)
# Get cell tags values
if hasattr(cell_tags, 'values'):
cell_values = cell_tags.values
else:
cell_values = cell_tags
for tag, props in material_properties.items():
cells = np.where(cell_values == tag)[0]
if len(cells) > 0:
E_val = props['E']
rho_val = props['rho']
nu_val = props['nu']
E.x.array[cells] = np.full_like(cells, E_val, dtype=default_scalar_type)
rho.x.array[cells] = np.full_like(cells, rho_val, dtype=default_scalar_type)
nu.x.array[cells] = np.full_like(cells, nu_val, dtype=default_scalar_type)
# Compute mu and lambda for these cells
mu_val = E_val / (2 * (1 + nu_val))
lambda_val = E_val * nu_val / ((1 + nu_val) * (1 - 2 * nu_val))
# Assign to mu and lambda functions
mu.x.array[cells] = np.full_like(cells, mu_val, dtype=default_scalar_type)
lambda_.x.array[cells] = np.full_like(cells, lambda_val, dtype=default_scalar_type)
print(f" Assigned to {len(cells)} cells for tag {tag}")
# Verify assignment
print(f"\nMaterial function summary:")
print(f" E_func: min={E.x.array.min():.2e}, max={E.x.array.max():.2e}")
print(f" rho_func: min={rho.x.array.min():.2e}, max={rho.x.array.max():.2e}")
print(f" nu_func: min={nu.x.array.min():.2e}, max={nu.x.array.max():.2e}")
# ------------Boundary tagging-----------------------------------------------------
dx = ufl.Measure("dx", domain=mesh)
# ------------Define weak form for linear elasticity -----------
# Create coordinate component
direction = 'z'
direction_map = {'x': 0, 'y': 1, 'z': 2}
d_idx = direction_map[direction]
# Unit vector in periodicity direction
e_d = ufl.as_vector([float(i == d_idx) for i in range(3)])
u_z = u[2]
partial_z_u = u.dx(2)
# Define weak forms
first_order_stress = ((lambda_ + mu) *(ufl.nabla_grad(u_z) + e_d * (ufl.nabla_div(u)) ) + 2*mu * partial_z_u)
a_linear = ufl.inner(first_order_stress, v) * dx
quadratic_stress = mu * u + (lambda_ + mu)* e_d * u_z
a_quadratic = ufl.inner(-quadratic_stress, v) * dx
# The two following terms belong to the K0 matrix, an intgration by part has been performed thanks to the divergence theorem and using grad(div(U)) = div(grad(U).T)
a_standard1 = -(lambda_ + mu) * ufl.nabla_grad(u).T
a_standard2 = -mu * ufl.nabla_grad(u)
a = (ufl.inner(a_standard1, ufl.nabla_grad(v)) + ufl.inner(a_standard2, ufl.nabla_grad(v))) * dx
m = rho * inner(u, v) * dx
# Set periodicity parameters
direction = 'z' # Periodicity along z-axis
direction_map = {'x': 0, 'y': 1, 'z': 2}
d_idx = direction_map[direction]
# Set basis vector
a_basis = np.zeros(3)
a_basis[d_idx] = L_f # L_f is the length in periodic direction
# Create indicator and relation functions for MPC
indicator = lambda x: periodic_indicator(x, direction=direction)
relation = lambda x: periodic_relation(x, direction=direction, a_basis=a_basis)
#-----------------MPC, assembling matrices and handling rigid body modes -------------------------------------------------------------------------------------
# Create MPC with phase=1 for periodic boundary conditions (no phase shift for the base matrices)
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, indicator=indicator, relation=relation, bcs=[], scale=1.0)
mpc.finalize()
print("Assembling matrices with periodic boundary conditions (phase=1)...")
# Assemble ALL matrices ONCE with periodic BCs (phase=1)
K0 = dolfinx_mpc.assemble_matrix(form(a), mpc)
K1 = dolfinx_mpc.assemble_matrix(form(a_linear), mpc)
K2 = dolfinx_mpc.assemble_matrix(form(a_quadratic), mpc)
M = dolfinx_mpc.assemble_matrix(form(m), mpc)
K0.assemble()
K1.assemble()
K2.assemble()
M.assemble()
print("Matrix assembly complete!")
rows_K0, cols_K0= K0.getSize()
rows_K1, cols_K1= K1.getSize()
rows_K2, cols_K2= K2.getSize()
rows_M, cols_M= M.getSize()
K0_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(K0, V, mpc=mpc)
K1_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(K1, V, mpc=mpc)
K2_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(K2, V, mpc=mpc)
M_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(M, V, mpc=mpc)
print('Done with extracting dof')
#-------------------Frequency sweep to solve for k-------------------------------------
num_f = 10 # Number of frequencies to sweep over
num_k = 5 # Minimum number of eigenvalues we want
f_name = 'f'
f_values = np.linspace(0, 10000, num_f)
dispersion_k_real = []
dispersion_k_imag = []
for idx, f in enumerate(f_values):
print(f"\n{'='*60}")
print(f"Frequency {idx+1}/{len(f_values)}: f = {f:.1f} Hz")
print(f"{'='*60}")
omega = 2 * np.pi * f
omega_sq = omega**2
# Build A0 = K0 - ω²M (using the matrices with periodic BCs)
A0 = K0_sub.copy()
A0 += -omega_sq * M_sub
A0.assemble()
pep = SLEPc.PEP().create()
pep.setOperators([A0, K1_sub, K2_sub]) #important: add bracket aroudn these to form "one" array, and it should go in increasing order
pep.setProblemType(SLEPc.PEP.ProblemType.HERMITIAN)
pep.setDimensions(nev=num_k)
k0 = 1
target = 1j * k0
pep.setTarget(target)
pep.setWhichEigenpairs(SLEPc.PEP.Which.TARGET_MAGNITUDE)
st = pep.getST()
st.setType(SLEPc.ST.Type.SINVERT)
ksp = st.getKSP()
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pep.setFromOptions()
pep.solve()
# Get converged eigenvalues
nconv = pep.getConverged()
print(f"Number of converged eigenvalues for f = {f} Hz is {nconv}")
Lambda_values = []
Lambda_values_real = []
Lambda_values_imag = []
for i in range(nconv):
Lambda = pep.getEigenpair(i)
Lambda_values.append(Lambda)
Lambda_values_real.append(Lambda.real)
Lambda_values_imag.append(Lambda.imag)
print(f"Mode {i}: Lambda = {Lambda.real:.6f} + {Lambda.imag:.6f}i")
dispersion_k_real.append(Lambda_values_imag)
minus_Lambda_values_real = []
for i in range(len(Lambda_values_real)):
minus_Lambda_values_real.append(-Lambda_values_real[i])
dispersion_k_imag.append(minus_Lambda_values_real)
# ----------------Save frequencies and k-values---------------------------------
with open(os.path.join('fenicsx_dump_inverse_cyl_v23.pkl'), 'wb') as f:
pickle.dump([f_name, num_f, f_values, dispersion_k_real, dispersion_k_imag], f)
# ------------------Extract values to plot------------------------------------
with open(os.path.join('fenicsx_dump_inverse_cyl_v23.pkl'),'rb') as f:
obj = pickle.load(f)
f_name = obj[0]
num_f = obj[1]
f_values = obj[2]
dispersion_k_real = obj[3]
dispersion_k_imag = obj[4]
L_f= 1
#----------------------Post-processing by scaling with L_f and pi-----------------------------
for i in range(len(dispersion_k_real)):
for j in range(len(dispersion_k_real[i])):
dispersion_k_real[i][j] = dispersion_k_real[i][j]* L_f/np.pi
for i in range(len(dispersion_k_imag)):
for j in range(len(dispersion_k_imag[i])):
dispersion_k_imag[i][j] = dispersion_k_imag[i][j]* L_f/np.pi
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# --------------------------Plotting dispersion curves----------------------------------
#Plot 1: Real wavenumber vs Frequency (frequency on y-axis)
for i, f in enumerate(f_values):
for k_real in dispersion_k_real[i]:
ax1.scatter(k_real, f, c='blue', s=15, alpha=0.6)
ax1.set_xlabel('Re(k) (1/m)')
ax1.set_ylabel('Frequency (Hz)')
ax1.set_title('Real Wavenumber')
ax1.grid(True, alpha=0.3)
# Plot 2: Imaginary wavenumber vs Frequency (frequency on y-axis)
for i, f in enumerate(f_values):
for k_imag in dispersion_k_imag[i]:
if abs(k_imag) < 10: # Filter extreme values
ax2.scatter(k_imag, f, c='red', s=15, alpha=0.6)
ax2.set_xlabel('Im(k) (1/m)')
ax2.set_ylabel('Frequency (Hz)')
ax2.set_title('Imaginary Wavenumber')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join('plots', 'inverse_dispersion_vertical.png'), dpi=300, bbox_inches='tight')
plt.show()
fig = plt.figure()
The STEP-file for the geometry named “Test_periodic_rod_geo_r025” can be found in the following Github repo: GitHub - NadineTab02/FEniCSx_Project · GitHub
My questions are:
- Can there be a problem with the formulation of the equations?
- Is the mpc-periodicity correctly implemented?
- Is there another solver method that should be used for these kinds of problems?
Any insights on correctly formulating and solving this quadratic eigenvalue problem in FEniCSx would be greatly appreciated!