Hi everyone,
I’m implementing an inverse dispersion solver in DOLFINx that computes complex wavenumbers k(ω) for waves in a periodic elastic 3D structure, using the SLEPc PEP-solver and Bloch periodic boundary conditions. The eigenvalues are correct and the dispersion curves look as expected. However, when plotting the mode shapes from the dispersion curve, the modes don’t look correct. Especially in the first branch, the modes don’t look accurate as there is only rigid body motion and no bending waves.
I have a related forward solver that computes ω(k) with the SLEPc EPS-solver for the same geometry, mesh, and material assignment, which yields a similar dispersion curve (with only real-values k-numbers). The mode shapes from that solver look accurate.
I’ve tried to extract the eigenvectors from the PEP-solver in different ways, plotting both the real part of the displacement as well as the magnitudes, but the problem persists.
The code for solving the inverse problem and generating the dataset with frequencies, eigenvalues (k), and eigenvectors is attached below, along with the code for plotting the dispersion curve and clicking on points to obtain the mode shapes. The STEP-file for the geometry is available at: GitHub - NadineTab02/FEniCSx_Project · GitHub
# 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
# Function for gmsh
def import_CAD_to_FEniCSx(input_file = "object_new_broken_cubic_case", output_file = "from_comsol", lc = 0.25):
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()
volumes = gmsh.model.getEntities(dim=3)
volume_tags = [v[1] for v in volumes]
if len(volume_tags) > 1:
gmsh.model.occ.fragment([(3, tag) for tag in volume_tags], [])
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
if len(volumes) == 1:
gmsh.model.addPhysicalGroup(3, [v[1] for v in volumes], tag=1)
else:
listvols1 = [1, 2, 3, 4, 5, 6, 7, 9]
gmsh.model.addPhysicalGroup(3, listvols1, tag=101)
listvols2 = [8, 10]
gmsh.model.addPhysicalGroup(3, listvols2, tag=201)
surfaces = gmsh.model.getEntities(dim=2)
listsurf1 = [4]
listsurf2 = [5]
listsurf3 = [1, 2, 3, 6]
gmsh.model.addPhysicalGroup(2, listsurf1, tag=1001)
gmsh.model.addPhysicalGroup(2, listsurf2, tag=2001)
gmsh.model.addPhysicalGroup(2, listsurf3, tag=3001)
gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
gmsh.write(f"{output_file}.msh")
mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
final_mesh = mesh_data.mesh
facet_tags = mesh_data.facet_tags
cell_tags = mesh_data.cell_tags
gmsh.option.setNumber("Geometry.Surfaces", 1)
gmsh.option.setNumber("Geometry.SurfaceLabels", 1)
gmsh.fltk.run()
gmsh.finalize()
mesh_info = [final_mesh, cell_tags, facet_tags]
return mesh_info, None
# 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:
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:
slave_local = np.array(mpc.is_slave, dtype=np.int32)
mask = np.ones(num_dofs_local + num_ghosts, dtype=np.int32)
mask[slave_local.astype(bool)] = 0
idx_set = np.flatnonzero(mask).astype(np.int32)
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--------------------
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}")
# Define material properties based on the tags from your import function
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}")
# ------------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)])
# Define weak forms
def epsilon(u):
return 0.5 * (ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def kappa(u):
i, j = ufl.indices(2)
Aij = u[i]*e_d[j] + e_d[i]*u[j]
A = ufl.as_tensor(Aij, (i,j))
return 0.5 * A
# Fourth-order elasticity tensor C
Id = ufl.Identity(dim)
indices = ufl.indices(4)
def delta_product(i, j, k, l):
return ufl.as_tensor(Id[i, j] * Id[k, l], indices)
i, j, k, l = indices
C = lambda_ * delta_product(i, j, k, l) + mu * (delta_product(i, k, j, l) + delta_product(i, l, k, j))
# Bilinear forms with explicit conjugation (complex mode)
# K matrix
a_K = C[i, j, k, l] * ufl.conj(epsilon(v)[i, j]) * epsilon(u)[k, l] * dx
# L matrix
a_linear = (-C[i, j, k, l] * ufl.conj(kappa(v)[i, j]) * epsilon(u)[k, l] +
C[i, j, k, l] * ufl.conj(epsilon(v)[i, j]) * kappa(u)[k, l]) * dx
# H matrix
a_quadratic = C[i, j, k, l] * ufl.conj(kappa(v)[i, j]) * kappa(u)[k, l] * dx
# Mass matrix
m = rho * ufl.conj(v[i]) * u[i] * dx
# For the standard stiffness
a = a_K
# 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_k = 20 # Minimum number of eigenvalues we want
f_name = 'f'
num_f = 20 # Number of frequencies to sweep over
f_values = np.linspace(0, 10000, num_f)
num_f = len(f_values)
dispersion_k_real = []
dispersion_k_imag = []
all_modes_displacements = [] # Store displacement arrays
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
A0 = K0_sub.copy()
A0 += -omega_sq * M_sub
A0.assemble()
A2 = K2_sub.copy()
A2 *= -1.0
pep = SLEPc.PEP().create(mesh.comm)
pep.setOperators([A0, K1_sub, A2])
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}")
# Storage for this frequency (initialize empty lists)
Lambda_values = []
Lambda_values_real = []
Lambda_values_imag = []
modes_at_f = [] # Store mode indices
for i in range(nconv):
# Create vector to hold eigenvector
vr = A0.createVecRight()
Lambda = pep.getEigenpair(i, vr)
Lambda_values.append(Lambda)
Lambda_values_real.append(Lambda.real)
Lambda_values_imag.append(Lambda.imag)
k_real = Lambda.imag
k_imag = -Lambda.real
print(f"Mode {i}: Lambda = {Lambda.real:.6f} + {Lambda.imag:.6f}i "
f"k = {k_real:.6f} + {k_imag:.6f}i")
# Create full vector in function space
eh = Function(V)
full_vec = eh.x.petsc_vec
full_vec.set(0.0)
full_vec.setValues(idx_set_row, vr.getArray())
full_vec.assemble()
dim = V.dofmap.index_map_bs
mpc.backsubstitution(eh)
u_complex = eh.x.array.copy()
u_complex_reshaped = u_complex.reshape(-1, dim)
# Real and imaginary parts
u_real = np.real(u_complex_reshaped)
u_imag = np.imag(u_complex_reshaped)
u_magnitude = np.abs(u_complex_reshaped)
modes_at_f.append({
'mode_index': i,
'k_real': k_real,
'k_imag': k_imag,
'displacement_real': u_real,
'displacement_imag': u_imag,
'displacement_magnitude': u_magnitude,
'displacement_complex': u_complex_reshaped,
})
# Append to dispersion data
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)
# Append mode shapes
all_modes_displacements.append(modes_at_f)
# ----------------Save frequencies and k-values---------------------------------
with open(os.path.join('fenicsx_dump_inverse_cyl_testingJune25.pkl'), 'wb') as f:
pickle.dump([
f_name, # 'f'
num_f, # number of frequencies
f_values, # frequency values
dispersion_k_real, # Re(k) for each frequency
dispersion_k_imag, # Im(k) for each frequency,
all_modes_displacements # mode shapes for later plotting
], f)
print("\nData saved successfully!")
# Minimum Working Example - Interactive plotting with point clicking
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import dolfinx
from dolfinx import mesh, fem, plot
from basix.ufl import element
from mpi4py import MPI
import gmsh
import dolfinx.mesh as msh
from dolfinx.io import gmsh as gmshio
import pickle
import os
#--------------Mesh import function (same as solver)--------------------------------
def import_CAD_to_FEniCSx(input_file = "object_new_broken_cubic_case", output_file = "from_comsol", lc = 0.25):
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()
volumes = gmsh.model.getEntities(dim=3)
volume_tags = [v[1] for v in volumes]
if len(volume_tags) > 1:
gmsh.model.occ.fragment([(3, tag) for tag in volume_tags], [])
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
if len(volumes) == 1:
gmsh.model.addPhysicalGroup(3, [v[1] for v in volumes], tag=1)
else:
listvols1 = [1, 2, 3, 4, 5, 6, 7, 9]
gmsh.model.addPhysicalGroup(3, listvols1, tag=101)
listvols2 = [8, 10]
gmsh.model.addPhysicalGroup(3, listvols2, tag=201)
surfaces = gmsh.model.getEntities(dim=2)
listsurf1 = [4]
listsurf2 = [5]
listsurf3 = [1, 2, 3, 6]
gmsh.model.addPhysicalGroup(2, listsurf1, tag=1001)
gmsh.model.addPhysicalGroup(2, listsurf2, tag=2001)
gmsh.model.addPhysicalGroup(2, listsurf3, tag=3001)
gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
gmsh.write(f"{output_file}.msh")
mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
final_mesh = mesh_data.mesh
facet_tags = mesh_data.facet_tags
cell_tags = mesh_data.cell_tags
gmsh.option.setNumber("Geometry.Surfaces", 1)
gmsh.option.setNumber("Geometry.SurfaceLabels", 1)
gmsh.fltk.run()
gmsh.finalize()
mesh_info = [final_mesh, cell_tags, facet_tags]
return mesh_info, None
#--------------Create mesh--------------------------------
step_filename, mesh_name = 'Test_periodic_rod_geo_r025', 'mesh_rod'
lc = 10
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
# Get length parameter L_f for scaling
z_coords = mesh.geometry.x[:, 2]
L_f = z_coords.max() - z_coords.min()
k_scaling = L_f / np.pi
#--------------Load data--------------------------------
print("Loading data from fenicsx_dump_inverse_cyl_testingJune25.pkl...")
with open('fenicsx_dump_inverse_cyl_testingJune25.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]
all_modes_displacements = obj[5]
print(f"Loaded {len(f_values)} frequencies")
print(f"Frequency range: {f_values[0]:.1f} - {f_values[-1]:.1f} Hz")
# Scale k values
print("Scaling k values...")
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] * k_scaling
dispersion_k_imag[i][j] = dispersion_k_imag[i][j] * k_scaling
#--------------Interactive plot with point clicking--------------------------------
print("\n" + "="*50)
print("INTERACTIVE MODE SHAPE PICKER")
print("Click on points in the dispersion curve to plot mode shapes")
print("="*50)
# Collect data for plotting
all_k = []
all_f = []
points_data = []
for i, f in enumerate(f_values):
for j, k in enumerate(dispersion_k_real[i]):
if k > 0: # Only positive k
all_k.append(k)
all_f.append(f)
points_data.append({
'k': k,
'f': f,
'f_idx': i,
'mode_j': j,
'mode_index': all_modes_displacements[i][j].get('mode_index', j) if j < len(all_modes_displacements[i]) else j
})
# Create interactive plot
fig, ax = plt.subplots(figsize=(10, 8))
sc = ax.scatter(all_k, all_f, c='blue', s=20, picker=True)
ax.set_xlabel(r'Re$(k)$ [$\pi$/a]')
ax.set_ylabel(r'$f$ [Hz]')
ax.set_xlim(0, 1)
ax.set_ylim(0, 70000)
ax.grid(True)
# Store clicked points
clicked_points = []
def on_pick(event):
ind = event.ind[0]
point = points_data[ind]
print(f"\n{'='*50}")
print(f"CLICKED POINT {len(clicked_points)+1}:")
print(f" f = {point['f']:.2f} Hz")
print(f" k = {point['k']:.6f} (π/a)")
print(f" f_idx = {point['f_idx']}")
print(f" mode_j = {point['mode_j']}")
print(f" mode_index = {point['mode_index']}")
print(f"{'='*50}")
clicked_points.append(point)
fig.canvas.mpl_connect('pick_event', on_pick)
plt.show()
#--------------Plot mode shapes for clicked points--------------------------------
print(f"\n{'='*50}")
print(f"Plotting mode shapes for {len(clicked_points)} clicked points...")
print(f"{'='*50}")
VISUAL_SCALE = 10
for idx, point in enumerate(clicked_points):
print(f"\nPlotting {idx+1}/{len(clicked_points)}: f = {point['f']:.2f} Hz, k = {point['k']:.6f}")
f_idx = point['f_idx']
f_actual = f_values[f_idx]
modes_at_f = all_modes_displacements[f_idx]
# Find mode by mode_index
mode_found = None
for mode in modes_at_f:
if mode.get('mode_index') == point['mode_index']:
mode_found = mode
break
# Try using mode_j as fallback
if mode_found is None:
print(f" Trying mode_j = {point['mode_j']} instead...")
for mode in modes_at_f:
if mode.get('mode_index') == point['mode_j']:
mode_found = mode
print(f" Found mode with mode_index = {point['mode_j']}")
break
if mode_found is None:
print(f" No mode found. Skipping...")
continue
u_plot = mode_found['displacement_real']
k_real_scaled = mode_found['k_real'] * k_scaling
k_imag_scaled = mode_found['k_imag'] * k_scaling
# Create grid and attach displacement
topology, cell_types, geometry = plot.vtk_mesh(mesh)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)
grid["u"] = u_plot
# Calculate displacement magnitude
u_mag = np.linalg.norm(u_plot, axis=1)
max_u = np.max(u_mag)
# Normalize for colorbar
if max_u > 0:
u_mag_normalized = u_mag / max_u
else:
u_mag_normalized = u_mag
grid["Normalized |u|"] = u_mag_normalized
# Warp the mesh
warped = grid.warp_by_vector("u", factor=VISUAL_SCALE)
warped["Normalized |u|"] = u_mag_normalized
# Create plotter
p = pv.Plotter(window_size=(1000, 800))
p.set_background('white')
# Show undeformed mesh as wireframe
p.add_mesh(grid, style="wireframe", color="lightgray", line_width=1)
# Show deformed mesh colored by normalized displacement
p.add_mesh(warped, scalars="Normalized |u|", cmap="plasma", show_scalar_bar=True,
scalar_bar_args={
'title': 'Normalized |u|',
'n_labels': 5,
'fmt': '%.2f'
})
# Add text info
p.add_text(f"Point {idx+1}: f = {f_actual:.1f} Hz\nk = {k_real_scaled:.4f} + {k_imag_scaled:.4f}i\nScale: {VISUAL_SCALE:.3f}x",
font_size=12, position='upper_left', color='black')
p.show_axes()
p.view_isometric()
p.show()
print("\nDone! All mode shapes plotted.")
I would be very grateful for any advice on how to solve this problem!
/Nadine
