Modal analysis using dolfin-x

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
4 Likes