# Eigenmodes / eigenvalues 3D calculation

Hi there,

I would like to calculate eigenvalues and eigenmodes of a 3D tetrahedrale mesh.

In their report the explain how they made a 3D eigenvalue solver.
It’s seem clear until I tried to use their codes.

I’ve always this error :
Process finished with exit code -1073741819 (0xC0000005)

It’s maybe caused by my python version or other stuff like that, so I’ve tried to remake the code myself, following their report…

Same error.

Is there a simpler solution to calculate eigenvalue of a 3D mesh ?
or a solution to solve this problem ?

Thanks

did you also find Sections 2.2 and 2.4.1 in the NGSolve tutorials on eigenvalue-solvers ?

https://docu.ngsolve.org/nightly/i-tutorials/unit-2.2-eigenvalues/pinvit.html

Joachim

1 Like

You might also find useful the SLEPc wrapper inside ngsPETSc, here is a tutorial:

1 Like

I’ve try to mix this tutorial with the work that you did with your student.
Unfortunately, I don’t get it.
I wrote this code, and I no longer understanding how to make it work from the line 53.

Can you help me ?

Thanks a lot.
Thomas

from ngsolve import *
import netgen.gui
from netgen.occ import *
import scipy.linalg
import math

# Variable
global E, nu, rho
E = 210000
nu = 0.3
rho = 7810

# Degree of Element
global EleDeg
EleDeg = 2

global mu
global lam
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))

# Geometry creation and processing functions
box = Box((-100, -10, -10), (100, 10, 10))
geo = OCCGeometry(box)

mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
fes = VectorH1(mesh, order=EleDeg, dirichlet=".*")
u, v = fes.TnT()

A = BilinearForm(fes)
A += lam * div(u) * div(v) * dx

M = BilinearForm(fes)
M += rho * InnerProduct(u, v) * dx

pre = Preconditioner(A, "multigrid", inverse="sparsecholesky")

A.Assemble()
M.Assemble()

maxNumIterations = 50
eps = 1e-6
num = 10

u = GridFunction(fes, multidim=num)

uvecs = MultiVector(u.vec, num)
vecs = MultiVector(u.vec, 2 * num)

for v in vecs[0:num]:
v.SetRandom()
uvecs[:] = pre * vecs[0:num]
lams = Vector(num * [1])
numIterations = 0
res = [1] * (maxNumIterations)

for numIterations in range(maxNumIterations):
while True:
numIterations += 1
vecs[0:num] = A.mat * uvecs - (M.mat * uvecs).Scale(lams)
vecs[num:2 * num] = pre * vecs[0:num]

r = InnerProduct(vecs[num], vecs[0])
for j in range(1, num):
tmp = InnerProduct(vecs[num + j], vecs[j])
if r < tmp:
r = tmp
res[numIterations - 1] = r

vecs[0:num] = uvecs
vecs.Orthogonalize()

asmall = InnerProduct(vecs, A.mat * vecs)
msmall = InnerProduct(vecs, M.mat * vecs)

ev, evec = scipy.linalg.eigh(a=asmall, b=msmall)
lams = Vector(ev[0:num])
print(numIterations, ":", [lam / math.pi**2 for lam in lams])

uvecs[:] = vecs * Matrix(evec[:, 0:num])

if abs(res[numIterations - 1]) < eps or numIterations >= maxNumIterations:
break

for j in range(num):
u.vecs[j][:] = 0.0
u.vecs[j].data += uvecs[j]

Draw(u, mesh, "mode")
SetVisualization(deformation=True)