Eigen solver syntax in dolfin-x

1

Hi, I am solving a modal analysis problem. I am not able to write the syntax of eigen solver in dolfin-x The MWE for the same has been attached.

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)
2

I made some changes to your code for the current dolfinx

import dolfinx
import numpy as np
from mpi4py import MPI

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

from dolfinx.mesh import MeshTags, locate_entities, create_box
from dolfinx.mesh import CellType, locate_entities_boundary
from dolfinx.fem import (Constant, DirichletBC, Function, LinearProblem, FunctionSpace, VectorFunctionSpace, 
                         locate_dofs_topological)
import ufl

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

from petsc4py import PETSc
from slepc4py import SLEPc

mesh = create_box(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, (PETSc.ScalarType(0), PETSc.ScalarType(0), PETSc.ScalarType(0)))
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()

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)

and I got the following results

Number of converged eigenpairs 12
24170102233.356747
24159609663.2646
24128225812.786083
24064285982.84818
24062146763.88163
24057006404.662098
24047436042.282703
24002756434.838116
23953242686.51013
23912724277.967308
23802203732.108154
23776394801.64928
3

This is the solution that I was also getting, but this is not correct. Frequency should come in increasing order. These are correct natural frequencies that is coming using dolfin with same problem statement.

Solid FE: 61.76 [Hz]
Solid FE: 122.20 [Hz]
Solid FE: 380.86 [Hz]
Solid FE: 720.322[Hz]
Solid FE: 746.56 [Hz]
Solid FE: 1040.998 [Hz]

With dolfin-x, it should come same

4

Hi Shubham,

The default setting of the slepc4py eigensolver solved for the largest eigenvalue and as such that may be why you are not getting the correct result. It is difficult for me to reproduce a similar result without an example code that produces the dolfin result. Nonetheless, based on an example from this webpage, I’ve got a similar result except for small differences in the first couples of eigenvalues. Below is the dolfinx code:

import sys
import dolfinx
import numpy as np
import slepc4py
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, DirichletBC, Function, VectorFunctionSpace,
                         locate_dofs_topological)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1.0, 1.0, 0.01])],
    [20, 20, 2],
    cell_type=CellType.tetrahedron)

E, nu = (70.0E9), (0.23)
rho = (2500)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = 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, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form = inner(sigma(u), epsilon(v))*dx
m_form = rho*dot(u, v)*dx


K = dolfinx.fem.assemble_matrix(k_form)
K.assemble()

M = dolfinx.fem.assemble_matrix(m_form)
M.assemble()

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(100)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=15)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

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

        # Calculation of eigenfrequency from real part of eigenvalue
        freq_3D = np.sqrt(abs(e_val.real))/2/np.pi
        print("Eigenfrequency: {0:8.5f} [Hz]".format(freq_3D))

Hope this helps.

1 Like
5

Hello sir,
I am attaching the dolfin code for the same problem as you have asked for it. I have to replicate the same using dolfin-x.

from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
%matplotlib inline
mesh = BoxMesh(Point(0.,0.,0.),Point(5,0.6,0.3), 25, 3, 2)
E, nu, rho = (2e11), (0.3), (7850)
mu,lmbda = E/2./(1+nu), E*nu/(1+nu)/(1-2*nu)
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)
def eps(v):
    return sym(grad(v))
V = VectorFunctionSpace(mesh, 'CG', degree=2)
u = TrialFunction(V)
v = TestFunction(V)
left = CompiledSubDomain("near(x[0], 0) && on_boundary")
bc = DirichletBC(V, Constant((0.,0.,0.)), left)
k_form = inner(sigma(u),eps(v))*dx
m_form = rho*dot(u,v)*dx
K = PETScMatrix()
assemble(k_form, tensor=K)
M = PETScMatrix()
assemble(m_form, tensor=M)
bc.zero(M), bc.apply(K)
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.

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


L, B, H = 5., 0.3, 0.6 #Dimensions of the beam
# Exact solution computation
from scipy.optimize import root
from math import cos, cosh
falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]
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

    # Beam eigenfrequency
    if i % 2 == 0: # exact solution should correspond to weak axis bending
        I_bend = H*B**3/12.
    else:          #exact solution should correspond to strong axis bending
        I_bend = B*H**3/12.
    freq_beam = alpha(i/2)**2*sqrt(E*I_bend/(rho*B*H*L**4))/2/pi

    print("Solid FE: {0:8.5f} [Hz]   Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
    
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode, 0)
    eig_l.append(freq_3D)
    eig_v.append(rx)

The results are:

Computing 40 first eigenvalues...
Solid FE:  9.83051 [Hz]   Beam theory:  9.78457 [Hz]
Solid FE: 19.45146 [Hz]   Beam theory: 19.56914 [Hz]
Solid FE: 60.61639 [Hz]   Beam theory: 61.31884 [Hz]
Solid FE: 114.64296 [Hz]   Beam theory: 122.63769 [Hz]
Solid FE: 118.82391 [Hz]   Beam theory: 171.69454 [Hz]
Solid FE: 165.68699 [Hz]   Beam theory: 343.38908 [Hz]
6

This is the closest that I can get to the dolfin solver. I’m not sure of the internal settings of dolfin to obtain a totally similar result. My current implementation contains a couple of eigenvalues that is equal to 1. In the implementation below, I just ignored those eigenvalues.

import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, DirichletBC, Function, VectorFunctionSpace,
                         locate_dofs_topological)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)


L, B, H = 5., 0.3, 0.6

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([L, H, B])],
    [25, 3, 2],
    cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = 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, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form = inner(sigma(u), epsilon(v))*dx
m_form = rho*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()

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()


# Exact solution computation
def falpha(x):
    return np.cos(x) * np.cosh(x) + 1


def alpha(n):
    return root(falpha, (2*n + 1)*np.pi/2)['x'][0]


ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))

            ef += 1

If you run this code in dolfinx you’ll get a similar result as you posted before. I’m not an extensive user of SLEPc nor am I an eigenvalue solver expert. So, I don’t know why we get those 1’s eigenvalue. Anyway, hope this helps. :smiley:

1 Like
7

Thank you so much for the help :smiley:

8

Sir, can you also tell the syntax for writing eigen vectors corresponding to these eigen values/frequencies. The dolfin code presented above had
r, c, rx, cx = eigensolver.get_eigenpair(i)
where rx is the eigen vector corresponding to all nodes of the mesh.

9

The vr PETSc vector should contains the real part of the eigenvector.

10

I tried with vr but not getting correct results. This the MWE for that

import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, DirichletBC, Function, VectorFunctionSpace,
                         locate_dofs_topological)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)


L, B, H = 5., 0.3, 0.6

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([L, H, B])],
    [25, 3, 2],
    cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = 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, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form = inner(sigma(u), epsilon(v))*dx
m_form = rho*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()

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()


# Exact solution computation
def falpha(x):
    return np.cos(x) * np.cosh(x) + 1


def alpha(n):
    return root(falpha, (2*n + 1)*np.pi/2)['x'][0]


eig_v=[]
ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
            eig_v.append(vr)
            ef += 1
print(np.array(eig_v).T)

With this code getting this result

array([[ 7.87128161e-05,  7.87128161e-05,  7.87128161e-05, ...,
         7.87128161e-05,  7.87128161e-05,  7.87128161e-05],
       [-4.79422279e-05, -4.79422279e-05, -4.79422279e-05, ...,
        -4.79422279e-05, -4.79422279e-05, -4.79422279e-05],
       [ 2.09577501e-06,  2.09577501e-06,  2.09577501e-06, ...,
         2.09577501e-06,  2.09577501e-06,  2.09577501e-06],
       ...,
       [-3.39323655e-03, -3.39323655e-03, -3.39323655e-03, ...,
        -3.39323655e-03, -3.39323655e-03, -3.39323655e-03],
       [-1.81604351e-02, -1.81604351e-02, -1.81604351e-02, ...,
        -1.81604351e-02, -1.81604351e-02, -1.81604351e-02],
       [-2.21205070e-03, -2.21205070e-03, -2.21205070e-03, ...,
        -2.21205070e-03, -2.21205070e-03, -2.21205070e-03]])

but the correct one from dolfin code presented above is
Screenshot 2022-01-06 083818

11

Hi Shubham,

Have you tried to plot the eigenmode and check whether it is correct?

Adeeb

12

Hello sir,
sorry for the late reply.
I tried writing eigenmode in xdmf. But getting some error. Can you please help.

import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, dirichletbc, Function, VectorFunctionSpace,
                         locate_dofs_topological,assemble_matrix, form)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)

L, B, H = 5., 0.3, 0.6

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([L, H, B])],
    [25, 3, 2],
    cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = 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, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form= form(inner(sigma(u), epsilon(v)) * dx)
K = assemble_matrix(k_form, bcs=[bc])
K.assemble()
m_form = form(rho*dot(u, v)*dx)
M = assemble_matrix(m_form, bcs=[bc])
M.assemble()
shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()


# Exact solution computation
def falpha(x):
    return np.cos(x) * np.cosh(x) + 1


def alpha(n):
    return root(falpha, (2*n + 1)*np.pi/2)['x'][0]


import dolfinx.io


eig_v=[]
ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
            eig_v.append(vr.copy())
            
            
            with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "eigenvector.xdmf", "w") as xdmf:
                xdmf.write_mesh(mesh)
                vr.name = "eigenvector"
                xdmf.write_function(vr.copy()) 
            ef += 1
13

You cannot write a PETSc vector directly to XDMF. You need to wrap it as a dolfinx.fem.Function:

ur = dolfinx.fem.Function(V)
ur.vector.array[:] = vr
xdmf.write_function(ur.copy())
2 Likes
14

Thanks sir,
I have modified the code, but getting only one vector in xdmf. I need to get all the eigenmodes in one xdmf file.

import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, dirichletbc, Function, VectorFunctionSpace,
                         locate_dofs_topological,assemble_matrix, form)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)

L, B, H = 5., 0.3, 0.6

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([L, H, B])],
    [25, 3, 2],
    cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = 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, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form= form(inner(sigma(u), epsilon(v)) * dx)
K = assemble_matrix(k_form, bcs=[bc])
K.assemble()
m_form = form(rho*dot(u, v)*dx)
M = assemble_matrix(m_form, bcs=[bc])
M.assemble()
shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()


# Exact solution computation
def falpha(x):
    return np.cos(x) * np.cosh(x) + 1


def alpha(n):
    return root(falpha, (2*n + 1)*np.pi/2)['x'][0]


import dolfinx.io


eig_v=[]
ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
            eig_v.append(vr.copy())
            
            
            with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "eigenvector.xdmf", "w") as xdmf: 
                xdmf.write_mesh(mesh)
                ur = dolfinx.fem.Function(V)
                ur.vector.array[:] = vr
                xdmf.write_function(ur.copy()) 
            ef += 1
15

You need to move the declartion of the XDMFFile outside of the for loop, or use the append mode. Currently you are rewriting it at every entry to the loop.

1 Like
16

I changed this. but every mode eigenvector is giving white colour contour only.
the code is

import dolfinx.io

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "eigenvector.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
eig_v=[]
ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:

    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
            

                
            ur = dolfinx.fem.Function(V)
            ur.vector.array[:] = vr
            xdmf.write_function(ur.copy()) 
            eig_v.append(vr.copy())
            ef += 1

the output is coming like this for all the modes

abcd
abcd1908×977 101 KB

17

See: Using Paraview for visualization — FEniCSx tutorial

18

File->Open.Xdmf3ReaderT->Filters->Alphabetic shows the list of options, which are not selectable. I am unable to click on ExtractBlock. The screenshot is attached.

new
new1913×1050 214 KB

19

Press «apply» before trying to extract the block (green buttom center left on your screenshot)

1 Like
20

I tried with that also, but there is no view after applying ExtractBlock.

another
another1911×993 86.1 KB

It may be possible that syntax for writing the file is not correct. you may please go through this.