Eigenmodes / eigenvalues 3D calculation

NGSolve
1

Hi there,

I would like to calculate eigenvalues and eigenmodes of a 3D tetrahedrale mesh.
I’ve read a very interssant work of 3 students :

image
image694×994 59.6 KB

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

2

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
3

You might also find useful the SLEPc wrapper inside ngsPETSc, here is a tutorial:
https://ngspetsc.readthedocs.io/en/uzerbinati-docs/SLEPcEPS/poisson.py.html

1 Like
4

Thanks for the link,
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 += 2 * mu * InnerProduct(0.5 * (grad(u) + grad(u).trans), 0.5 * (grad(v) + grad(v).trans)) * dx
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)