Modal analysis using dolfin-x

dolfinx
,
1

Hello, I am trying to do modal analysis using dolfin-x. I have assembled mass and stiffness matrix. but not able to use eigen solver in it. I have to compute natural frequencies and mode shapes of a system. The below code was used for computing both parameters using dolfin. I want to replicate these steps in dolfin-x.

eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.

N_eig = 6   # number of eigenvalues
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)



# Set up file for exporting results
file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

# Extraction
eig_v = []
eig_l = []
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    # 3D eigenfrequency
    freq_3D = sqrt(r) /2/pi
    print("{0:02d}. Solid FE: {1:8.5f} [Hz] ".format(i+1,freq_3D))
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector {0:02d}".format((i+1)))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode, 0)
    eig_l.append(freq_3D)
    eig_v.append(rx)
2

Please supply your attempt of implementing it in DOLFINx, and explain what doesn’t work with your current adaptation.

3

I am using this for computing first 10 natural frequencies. This was used for computing eigen values for a eigen value problem, where there was only single matrix. Here we have generalized eigen value problem, where there are K and M matrices.

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(K,M)
eigensolver.setDimensions(nev=10)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.getVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    for i in range (evs):
            l = eigensolver.getEigenpair(i ,vr, vi)
            print(l.real)
4

Please add the full minimal code (of the DOLFINx code) that you are trying to use, not just code snippets.

5 
import dolfinx
import numpy as np
from mpi4py import MPI

import dolfinx.fem as fem, dolfinx.generation as generation
from petsc4py import PETSc
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
                 div, dx, ds, grad, inner)

from dolfinx.generation import BoxMesh
from dolfinx.mesh import MeshTags, locate_entities
from dolfinx.mesh import CellType, locate_entities_boundary
from dolfinx.fem import (Constant, DirichletBC, Function, LinearProblem, FunctionSpace, VectorFunctionSpace, 
                         locate_dofs_topological)
import ufl
mesh = BoxMesh(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([5, 0.6, 0.3])], [25,3,2], cell_type=CellType.hexahedron)
E, nu = (2e11), (0.3)  
rho = (7850) 
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)
def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
V = VectorFunctionSpace(mesh, ("CG", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_D = Function(V)
with u_D.vector.localForm() as loc:
    loc.set(0)
bc = DirichletBC(u_D, locate_dofs_topological(V, fdim, boundary_facets))

T = Constant(mesh, (0, 0, 0))
import ufl
ds = ufl.Measure("ds", domain=mesh)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
k_form = inner(sigma(u),epsilon(v))*dx
m_form = rho*ufl.dot(u,v)*dx


K = dolfinx.fem.assemble_matrix(k_form, bcs=[bc])
K.assemble()

M = dolfinx.fem.assemble_matrix(m_form, bcs=[bc])
M.assemble()

import sys, slepc4py
slepc4py.init(sys.argv)

from petsc4py import PETSc
from slepc4py import SLEPc

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(K,M)
eigensolver.setDimensions(nev=10)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.getVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    for i in range (evs):
            l = eigensolver.getEigenpair(i ,vr, vi)
            print(l.real)
6

Hi Shubham,

Did you ever get this working?

I’m interested in doing some (EM) eigenmode calcs, and am looking for (any) example code of how to do an eigenmode simulation to start to get my teeth into it.

Thanks in advance!
Gary

7

Hi Gary,
I could not figure it out.

Thanks

8

Here is a code to solve the eigenmodes of a clamped-free beam in DOLFINx:

# Imports
import numpy as np
import ufl
import sys, slepc4py
slepc4py.init(sys.argv)

from dolfinx import fem
from dolfinx.fem import (Constant, dirichletbc, Function, VectorFunctionSpace,
        locate_dofs_topological)
from dolfinx.mesh import (create_box, meshtags, locate_entities, CellType,
        locate_entities_boundary)
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

# Define computational domain
L = np.array([5, 0.6, 0.4])
N = [25, 3, 2]
mesh = create_box(MPI.COMM_WORLD, [np.array([0,0,0]), L], N,
        cell_type=CellType.hexahedron)
        
# Material constants
E, nu = (2e11), (0.3)  
rho = (7850) 
mu = Constant(mesh, E/2./(1+nu))
lambda_ = Constant(mesh, E*nu/(1+nu)/(1-2*nu))

# Convenience functions
def epsilon(u):
    return ufl.sym(ufl.grad(u))
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

# Define function space, trial and test functions
V = VectorFunctionSpace(mesh, ("CG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define Dirichlet boundary condition
def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_D = Function(V)
u_D.interpolate(lambda x: np.zeros_like(x))
bc = dirichletbc(u_D, locate_dofs_topological(V, fdim, boundary_facets))

# Define variational form
k_form = ufl.inner(sigma(u),epsilon(v))*ufl.dx
m_form = rho*ufl.inner(u,v)*ufl.dx

# Assemble stiffness and mass matrices
#
# Using the "diagonal" kwarg ensures that Dirichlet BC modes will not be among
# the lowest-frequency modes of the beam. 
K = fem.petsc.assemble_matrix(fem.form(k_form), bcs=[bc], diagonal=62831)
M = fem.petsc.assemble_matrix(fem.form(m_form), bcs=[bc], diagonal=1/62831)
K.assemble()
M.assemble()

# Create and configure eigenvalue solver
N_eig = 6
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(N_eig)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
st = SLEPc.ST().create(MPI.COMM_WORLD)
st.setType(SLEPc.ST.Type.SINVERT)
st.setShift(0.1)
st.setFromOptions()
eigensolver.setST(st)
eigensolver.setOperators(K, M)
eigensolver.setFromOptions()

# Compute eigenvalue-eigenvector pairs
eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.getVecs()
u_output = Function(V)
u_output.name = "Eigenvector"
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    with XDMFFile(MPI.COMM_WORLD, "eigenvectors.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        for i in range (min(N_eig, evs)):
            l = eigensolver.getEigenpair(i, vr, vi)
            freq = np.sqrt(l.real)/2/np.pi
            print(f"Mode {i}: {freq} Hz")
            u_output.x.array[:] = vr
            xdmf.write_function(u_output, i)

Tested with the latest DOLFINx Docker container, which is currently:

python3 -c 'import dolfinx;print(dolfinx.__version__);print(dolfinx.git_commit_hash)'

0.4.2.0
040fa6553c0cf746040a79adc5ca886df8ec0111
9

For future visitors to this thread, Here is the script from conpierce8 updated to run on fenicsx circa 2025:

# Imports
import numpy as np
import ufl
import sys, slepc4py
slepc4py.init(sys.argv)

from dolfinx import fem
from dolfinx.fem import (Constant, dirichletbc, Function,
        locate_dofs_topological)
from dolfinx.mesh import (create_box, CellType,
        locate_entities_boundary)
from dolfinx.io import XDMFFile
from mpi4py import MPI
from slepc4py import SLEPc
from dolfinx.fem.petsc import assemble_vector, assemble_matrix


# Define computational domain
L = np.array([5, 0.6, 0.4])
N = [25, 3, 2]
mesh = create_box(MPI.COMM_WORLD, [np.array([0,0,0]), L], N,
        cell_type=CellType.hexahedron)
        
# Material constants
E, nu = (2e11), (0.3)  
rho = (7850) 
mu = Constant(mesh, E/2./(1+nu))
lambda_ = Constant(mesh, E*nu/(1+nu)/(1-2*nu))

# Convenience functions
def epsilon(u):
    return ufl.sym(ufl.grad(u))
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(3) + 2*mu*epsilon(u)

# Define function space, trial and test functions
V = fem.functionspace(mesh, ("Lagrange", 1, (3, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define Dirichlet boundary condition
def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_D = Function(V)
u_D.interpolate(lambda x: x * 0)
bc = dirichletbc(u_D, locate_dofs_topological(V, fdim, boundary_facets))

# Define variational form
k_form = ufl.inner(sigma(u),epsilon(v))*ufl.dx
m_form = rho*ufl.inner(u,v)*ufl.dx

# Assemble stiffness and mass matrices
#
# Using the "diagonal" kwarg ensures that Dirichlet BC modes will not be among
# the lowest-frequency modes of the beam. 
M = assemble_matrix(fem.form(m_form), bcs=[bc], diagonal=1/62831)
K = assemble_matrix(fem.form(k_form), bcs=[bc], diagonal=62831)

K.assemble()
M.assemble()

# Create and configure eigenvalue solver
N_eig = 6
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(N_eig)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
st = SLEPc.ST().create(MPI.COMM_WORLD)
st.setType(SLEPc.ST.Type.SINVERT)
st.setShift(0.1)
st.setFromOptions()
eigensolver.setST(st)
eigensolver.setOperators(K, M)
eigensolver.setFromOptions()

# Compute eigenvalue-eigenvector pairs
eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.getVecs()
u_output = Function(V)
u_output.name = "Eigenvector"
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    with XDMFFile(MPI.COMM_WORLD, "eigenvectors.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        for i in range (min(N_eig, evs)):
            l = eigensolver.getEigenpair(i, vr, vi)
            freq = np.sqrt(l.real)/2/np.pi
            print(f"Mode {i}: {freq} Hz")
            u_output.x.array[:] = vr
            xdmf.write_function(u_output, i)
10

Hi, I was trying to use this code and unfortunately I ran into an error when compiling the K matrix. When it reaches that line, it spits out an error code related to “incompatible function arguments.” I have been looking for something that talks about this error, yet I haven’t encountered anything that solves the issue as relatable ones were seemingly resolved by items which are featured that do exist within the current code (compiling with fem.form, assembly after running the petsc command). Maybe the dolfinx lab I pulled recently causes the code to become outdated (it says it was created six months ago), but I have no idea. Wish I could submit an screenshot as to not clutter too badly but I cannot do so as a new user, so please bear with me.
The lines in question are:

k = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
m = rho * ufl.inner(u, v) * ufl.dx

k_compiled = dolfinx.fem.form(k)
K = dolfinx.fem.petsc.assemble_matrix(k_compiled, bcs=bcs, diag=1)
K.assemble()

specifically, K = dolfinx.fem.petsc.assemlbe_matrix(k_compiled,bcs=bcs,diag=1).

When running the code it throws the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 72
     69 m = rho * ufl.inner(u, v) * ufl.dx
     71 k_compiled = dolfinx.fem.form(k)
---> 72 K = dolfinx.fem.petsc.assemble_matrix(k_compiled, bcs=bcs, diag=1)
     73 K.assemble()
     75 m_compiled = dolfinx.fem.form(m)

File /usr/lib/python3.12/functools.py:909, in singledispatch.<locals>.wrapper(*args, **kw)
    905 if not args:
    906     raise TypeError(f'{funcname} requires at least '
    907                     '1 positional argument')
--> 909 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:400, in assemble_matrix(a, bcs, diag, constants, coeffs, kind)
    350 """Assemble a bilinear form into a matrix.
    351 
    352 The following cases are supported:
   (...)    397     Matrix representing the bilinear form.
    398 """
    399 A = create_matrix(a, kind)
--> 400 assemble_matrix(A, a, bcs, diag, constants, coeffs)  # type: ignore[arg-type]
    401 return A

File /usr/lib/python3.12/functools.py:909, in singledispatch.<locals>.wrapper(*args, **kw)
    905 if not args:
    906     raise TypeError(f'{funcname} requires at least '
    907                     '1 positional argument')
--> 909 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:506, in _(A, a, bcs, diag, constants, coeffs)
    504 coeffs = pack_coefficients(a) if coeffs is None else coeffs  # type: ignore[assignment]
    505 _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else []
--> 506 _cpp.fem.petsc.assemble_matrix(A, a._cpp_object, constants, coeffs, _bcs)
    507 if a.function_spaces[0] is a.function_spaces[1]:
    508     A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore[attr-defined]

TypeError: assemble_matrix(): incompatible function arguments. The following argument types are supported:
    1. assemble_matrix(A: mat, a: dolfinx.cpp.fem.Form_complex128, constants: ndarray[dtype=complex128, shape=(*), order='C', writable=False], coeffs: collections.abc.Mapping[tuple[dolfinx.cpp.fem._IntegralType, int], ndarray[dtype=complex128, shape=(*, *), order='C', writable=False]], bcs: collections.abc.Sequence[dolfinx.cpp.fem.DirichletBC_complex128], unrolled: bool = False) -> None
    2. assemble_matrix(A: mat, a: dolfinx.cpp.fem.Form_complex128, constants: ndarray[dtype=complex128, shape=(*), order='C', writable=False], coeffs: collections.abc.Mapping[tuple[dolfinx.cpp.fem._IntegralType, int], ndarray[dtype=complex128, shape=(*, *), order='C', writable=False]], rows0: ndarray[dtype=int8, shape=(*), order='C', writable=False], rows1: ndarray[dtype=int8, shape=(*), order='C', writable=False], unrolled: bool = False) -> None

Invoked with types: Mat, dolfinx.cpp.fem.Form_complex128, ndarray, dict, list

I’m pretty new to docker/dolfinx, so I understand if this is classifiable as a “skill issue,” or if I am a bit out of my league. For context I am trying to get the bending modes of a plate of arbitrary shape given a fixed edge at y=0, so the code here feels like it would give me some direction.

EDIT: I messed up and posted this in the wrong thread lol. I had meant to put this in a different thread that also had modal analysis code. Nevertheless, running either code returns the same issue and I cannot really see what else breaks past the aforementioned problematic lines.

11

To run on the DOLFINx main branch, I have no trouble running the code once one changes it from:

to

M = assemble_matrix(fem.form(m_form), bcs=[bc], diag=1/62831)
K = assemble_matrix(fem.form(k_form), bcs=[bc], diag=62831)

The code by @MattFerraro returns

Number of converged eigenpairs 8
Mode 0: 13.946300071906201 Hz
Mode 1: 20.08439252108583 Hz
Mode 2: 85.25303934739036 Hz
Mode 3: 118.83155190126726 Hz
Mode 4: 141.7986866065624 Hz
Mode 5: 230.52964645985446 Hz

However, when i look at your output, I see that you have:

which indicates that you are running with complex numbers, which means that you should wrap the constants and diag value as complex numbers, i.e.

# Material constants
from dolfinx import default_scalar_type
E, nu = (2e11), (0.3)  
rho = (7850) 
mu = Constant(mesh, default_scalar_type(E/2./(1+nu)))
lambda_ = Constant(mesh, default_scalar_type(E*nu/(1+nu)/(1-2*nu)))
```
and
```python
M = assemble_matrix(fem.form(m_form), bcs=[bc], diag=default_scalar_type(1/62831))
K = assemble_matrix(fem.form(k_form), bcs=[bc], diag=default_scalar_type(62831))
```
12

It worked, thanks for the help!