Hi everyone, I’m facing problem in performing Bloch analysis for wave propagation; transformation matrix is giving error.
import numpy as np
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction
from mpi4py import MPI
import dolfinx.fem as fem
import dolfinx.fem.petsc as dolfinx_petsc # <--- alias for convenience
from dolfinx.mesh import create_rectangle, CellType
length, height = 1.0, 1.0
Nx, Ny = 5, 5
domain = create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([length, height])],
[Nx, Ny],
cell_type=CellType.triangle,
)
dim = domain.topology.dim
print(f"Mesh topology dimension d={dim}.")
degree = 2
shape = (dim,) # vector field of size `dim`
V = fem.functionspace(domain, ("P", degree, shape)) # using fem.functionspace (lowercase) as in newer dolfinx
u_sol = fem.Function(V, name="Displacement")
# material constants (use dtype-double arrays where necessary)
E = fem.Constant(domain, 210e3)
nu = fem.Constant(domain, 0.3)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
def epsilon(v):
return sym(grad(v))
def sigma(v):
return lmbda * tr(epsilon(v)) * Identity(dim) + 2 * mu * epsilon(v)
# variational forms
u = TrialFunction(V)
v = TestFunction(V)
rho = 2e-3
g = 9.81
f = fem.Constant(domain, np.array([0.0, -rho * g]))
dx = Measure("dx", domain=domain)
a = inner(sigma(u), epsilon(v)) * dx
L = inner(f, v) * dx
m_form = rho * inner(u, v) * dx
# geometric boundary markers (x is shape (dim, Npoints))
def left(x):
return np.isclose(x[0], 0.0)
def right(x):
return np.isclose(x[0], length)
def bottom(x):
return np.isclose(x[1], 0.0)
def top(x):
return np.isclose(x[1], height)
left_dofs = fem.locate_dofs_geometrical(V, left)
right_dofs = fem.locate_dofs_geometrical(V, right)
bottom_dofs = fem.locate_dofs_geometrical(V, bottom)
top_dofs = fem.locate_dofs_geometrical(V, top)
print(left_dofs,'\n',right_dofs)
# Example: create a Dirichlet BC (zero displacement) on the left side
zero = np.zeros((dim,), dtype=np.double)
bc_left = fem.dirichletbc(zero, left_dofs, V)
# bcs must be a list of DirichletBC objects (or an empty list)
bcs = [] # if no BCs wanted, set bcs = []
# IMPORTANT: assemble_matrix expects a dolfinx form object; wrap with fem.form
A = dolfinx_petsc.assemble_matrix(fem.form(a), bcs=bcs)
A.assemble()
M = dolfinx_petsc.assemble_matrix(fem.form(m_form), bcs=bcs)
M.assemble()
from petsc4py import PETSc
import scipy.sparse as sp
# Convert PETSc Mat to SciPy CSR
def petsc_to_scipy(A_petsc: PETSc.Mat):
# Convert to AIJ (CSR) format
A_csr = A_petsc.convert("aij")
ia, ja, a = A_csr.getValuesCSR()
return sp.csr_matrix((a, ja, ia), shape=A_csr.getSize())
# Convert PETSc Mat to dense NumPy
def petsc_to_numpy(A_petsc: PETSc.Mat):
A_dense = A_petsc.convert("dense")
return A_dense.getDenseArray()
# ---- Convert to SciPy CSR ----
A_csr = petsc_to_scipy(A)
M_csr = petsc_to_scipy(M)
# ---- Convert to NumPy dense ----
A_np = petsc_to_numpy(A)
M_np = petsc_to_numpy(M)
# Print matrix sizes
A_size = A.getSize()
M_size = M.getSize()
print("Stiffness matrix A size:", A_size)
print("Mass matrix M size:", M_size)
ndofs = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print("Number of DOFs:", ndofs)
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
# --- assume you have A_np, M_np as dense or CSR matrices ---
def build_TTBLOCH(V, left_dofs, right_dofs, bottom_dofs, top_dofs, kvec, length, height):
dim = V.dofmap.index_map_bs
ndofs = V.dofmap.index_map.size_global * dim
T = np.eye(ndofs, dtype=np.complex128)
dof_coords = V.tabulate_dof_coordinates()
node_ids = np.arange(ndofs) // dim
node_to_dofs = {i: np.where(node_ids == i)[0] for i in range(ndofs // dim)}
left_nodes = np.unique(node_ids[left_dofs])
right_nodes = np.unique(node_ids[right_dofs])
bottom_nodes = np.unique(node_ids[bottom_dofs])
top_nodes = np.unique(node_ids[top_dofs])
# --- X-direction: right -> left ---
phase_x = np.exp(1j * kvec[0] * length)
for rn, ln in zip(right_nodes, left_nodes):
rn_dofs = node_to_dofs[rn]
ln_dofs = node_to_dofs[ln]
for dr, dl in zip(rn_dofs, ln_dofs):
T[dr, :] = 0.0
T[dr, dl] = phase_x
# --- Y-direction: top -> bottom (skip corner nodes already mapped in X) ---
right_set = set(right_nodes)
top_nodes_noncorner = [n for n in top_nodes if n not in right_set]
phase_y = np.exp(1j * kvec[1] * height)
for tn, bn in zip(top_nodes_noncorner, bottom_nodes):
tn_dofs = node_to_dofs[tn]
bn_dofs = node_to_dofs[bn]
for dt, db in zip(tn_dofs, bn_dofs):
T[dt, :] = 0.0
T[dt, db] = phase_y
return T
# High-symmetry path in reciprocal space (2D square lattice)
k_points = [
np.array([0.0, 0.0]), # Gamma
np.array([0.0, np.pi]), # X
np.array([np.pi, np.pi]), # M
np.array([np.pi, 0.0]), # Y
np.array([0.0, 0.0]) # Gamma
]
band_eigs = []
for kvec in k_points:
T = build_TTBLOCH(V, left_dofs, right_dofs, bottom_dofs, top_dofs, kvec, length, height)
# reduced Bloch matrices
Kb = T.conj().T @ A_np @ T
Mb = T.conj().T @ M_np @ T
# Solve generalized eigenproblem (few lowest modes)
nev = 6
vals, vecs = spla.eigsh(Kb, k=nev, M=Mb, sigma=0.0, which='LM')
omegas = np.sqrt(np.real(vals))
band_eigs.append(omegas)
# Convert to array: (n_kpoints, nev)
band_eigs = np.array(band_eigs)
print("Band structure frequencies:\n", band_eigs)