3D beam eigenvalues

NGSolve
1

Hello everyone,
I wrote a code in order to calculate eigenvalue for a 3D beam using NGSolve. I got inspired by Sebastian Hirnschall, Christoph Leonhartsberger and Rafael Dorigo’s work “Computing igenfrequencies using Finite Element Methods”.
I have a problem which is having negative values in stuiffness matrix when computing “scipy.linalg.eigh”.
Did any of you have had this problem please ? Could you please help me solve the problem I appreciate in advance.
Here is the code:

from netgen.csg import *
from ngsolve import *
import math
import scipy.linalg
from numpy import *
from netgen.NgOCC import *
from netgen.occ import *

def getF(l,shift=0):
        return (sqrt(abs(l - shift)) / (2 * math.pi))

maxh1 = 10
E, nu, rho = 210000, 0.3, 7810e-9
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
maxNumIterations = 50
eps = 1e-6
geometryPath = "Geometrie.stp"
geometry=OCCGeometry(geometryPath)
mesh = Mesh(geometry.GenerateMesh(maxh=maxh1)).Curve(3)

fes = VectorH1(mesh, order=3, dirichlet="Face3")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(2 * mu * InnerProduct(0.5 * (grad(u) + grad(u).trans), 0.5 * (grad(v) + grad(v).trans)) * dx + lam * div(u) * div(v) * dx)
pre = Preconditioner(a, "multigrid", inverse="sparsecholesky")
m = BilinearForm(rho * u * v * dx)
a.Assemble()
m.Assemble()


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])
res = [1] * (maxNumIterations)
numIterations = 0
while True:
    numIterations += 1
    vecs[0:num] = a.mat * uvecs[0:num] - (m.mat * uvecs[0:num]).Scale(lams)
    vecs[num:2 * num] = pre * vecs[0:num]

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

    vecs[0:num] = uvecs[0:num]
    vecs.Orthogonalize()
    print(vecs[0])
    asmall = InnerProduct(vecs, a.mat * vecs)
    msmall = InnerProduct(vecs, m.mat * vecs)
    print(asmall)   
    ev, evec = scipy.linalg.eigh(a=asmall, b=msmall)
    lams = Vector(ev[0:num])
    print(numIterations, ":", [getF(l, 0) for l in lams])
    print("res:", res[numIterations-1], "\n")
    if (numIterations==1):
        tmp = MultiVector(u.vec, 2 * num)
        tmp[0:2*num] = vecs
        vecs = MultiVector(u.vec, 3 * num)
        vecs[0:2*num] = tmp[0:2 * num]

    uvecs[0:num] = vecs * Matrix(evec[:, 0:num])
    vecs[2 * num:3 * num] = vecs[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)

Thank you
Sabri

2

can you also add the geometry to replicate the problem?

3

Hello Christopher, here is the code with the good format and the beam geometry.

from netgen.csg import *
from ngsolve import *
import math
import scipy.linalg
from numpy import *
from netgen.occ import *

def getF(l,shift=0):
        return (sqrt(abs(l - shift)) / (2 * math.pi))

maxh1 = 10
E, nu, rho = 210000, 0.3, 7810e-9
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
maxNumIterations = 50
eps = 1e-6

box = Box((0,0,0), (3,0.6,1))
box.faces.Min(X).name = "fix"
mesh = Mesh(OCCGeometry(box).GenerateMesh(maxh=0.1)).Curve(3)

fes = VectorH1(mesh, order=3, dirichlet="fix")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(2 * mu * InnerProduct(0.5 * (grad(u) + grad(u).trans), 0.5 * (grad(v) + grad(v).trans)) * dx + lam * div(u) * div(v) * dx)
pre = Preconditioner(a, "multigrid", inverse="sparsecholesky")
m = BilinearForm(rho * u * v * dx)
a.Assemble()
m.Assemble()

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])
res = [1] * (maxNumIterations)
numIterations = 0
while True:
    numIterations += 1
    vecs[0:num] = a.mat * uvecs[0:num] - (m.mat * uvecs[0:num]).Scale(lams)
    vecs[num:2 * num] = pre * vecs[0:num]

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

    res[numIterations - 1] = s

    vecs[0:num] = uvecs[0:num]
    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, ":", [getF(l, 0) for l in lams])
    print("res:", res[numIterations-1], "\n")

    if (numIterations==1):
        tmp = MultiVector(u.vec, 2 * num)
        tmp[0:2*num] = vecs
        vecs = MultiVector(u.vec, 3 * num)
        vecs[0:2*num] = tmp[0:2 * num]

    uvecs[0:num] = vecs * Matrix(evec[:, 0:num])
    vecs[2 * num:3 * num] = vecs[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)

Thank you in advance.
Sabri

4

I think you run into numerical instabilities of scipy eigh.
What works for me (I run then later down in line 76 in some non matching matrix sizes) is orthogonalizing the vecs wrt the mass matrix:

vecs.Orthogonalize(m.mat)

but I think scaling alone should also help.

5

Hello Christopher,
Thank you for your help.
I got rid of the line causing the new error and it worked for me.
Thank you again
Best regards
Sabri

6

Dear Christopher,

I have a problem vizualising the eigenmodes using NETGEN. All I get is a grey drawing of the mesh eventhough the frequencies are correct and the script works fine.

image
image1869×821 124 KB

Do you, please, have an idea where the problem could come from ? I appreciate your help in advance.
Sabri

7

Can you send your mwe?

8 
from ngsolve import *
import netgen.gui
from netgen.NgOCC import *
import scipy.linalg
import math

def getF(l,shift=0):
        return (sqrt(abs(l - shift)) / (2 * math.pi))

maxh1 = 10
E, nu, rho = 210000, 0.3, 7810e-9
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
maxNumIterations = 50
eps = 1e-3

box = Box((0,0,0), (3,0.6,1))
box.faces.Min(X).name = "fix"
mesh = Mesh(OCCGeometry(box).GenerateMesh(maxh=maxh1)).Curve(3)

fes = VectorH1(mesh, order=3, dirichlet="fix")

u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(2 * mu * InnerProduct(0.5 * (grad(u) + grad(u).trans), 0.5 * (grad(v) + grad(v).trans)) * dx + lam * div(u) * div(v) * dx)
pre = Preconditioner(a, "multigrid", inverse="sparsecholesky")
m = BilinearForm(rho * u * v * dx)
a.Assemble()
m.Assemble()

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])
res = [1] * (maxNumIterations)
numIterations = 0
while True:
    numIterations += 1
    vecs[0:num] = a.mat * uvecs[0:num] - (m.mat * uvecs[0:num]).Scale(lams)
    vecs[num:2 * num] = pre * vecs[0:num]

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

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

    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, ":", [getF(l, 0) for l in lams])
    print("res:", res[numIterations-1], "\n")
    
    uvecs[0:num] = 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)
input('press esc to quit \n')```
9

go to visual menu
select clipping plane sol: Vector Funtion
then go clipping menu
select enable clipping
and then with the multidim component selector you can select the eigenfunctions to the corresponding eigenvalue

image
image613×352 16.5 KB

10

Thank you Christopher, I activated Draw Surface Vectors in visualization menu and it worked. I really appreciate your help.
Best regards.
Sabri