# Is it possible to do Vibroacoustics ?

More precisely i want to know if there is any pitfall for the implementation of something like seen in this elmer fem video.

It would be great if i can use ngsolve for this task

Hi Kai,

I guess you want to solve it in two steps ?
Step 1 eigenvalue problem for the tuning forg.
Step 2 solve the acoustic problem

For step 1 you may watch this video from a studentâ€™s seminar:

For step 2 you may want to use PML for the radiation condition, see
https://docu.ngsolve.org/latest/i-tutorials/unit-1.7-helmholtz/pml.html

Best, Joachim

Ok, i have played around with some examples where the wave equation is used. I finally end up with a code that have the boundaries of a rectangle as source for the wave equation.

[code][code][color=#0033b3]from [/color]ngsolve [color=#0033b3]import [/color]*
[color=#0033b3]from [/color]netgen.geom2d [color=#0033b3]import [/color]SplineGeometry

[color=#8c8c8c]# Geometry[/color]
geo = SplineGeometry()
[color=#660099]leftdomain[/color]=[color=#1750eb]0[/color], [color=#660099]rightdomain[/color]=[color=#1750eb]1[/color], [color=#660099]bc[/color]=[color=#067d17]â€śscatâ€ť[/color])
mesh = Mesh(geo.GenerateMesh([color=#660099]maxh[/color]=[color=#1750eb]0.1[/color]))

fes = H1(mesh, [color=#660099]order[/color]=[color=#1750eb]5[/color], [color=#660099]complex[/color]=[color=#0033b3]True[/color])
u, v = fes.TnT()

[color=#8c8c8c]# Wavenumber & source[/color]
omega = [color=#1750eb]15[/color]
pulse = [color=#1750eb]5e4[/color]exp(-([color=#1750eb]40[/color]**[color=#1750eb]2[/color])((x-[color=#1750eb]0.5[/color])(x-[color=#1750eb]0.5[/color]) + (y-[color=#1750eb]0.5[/color])(y-[color=#1750eb]0.5[/color])))

[color=#8c8c8c]# Forms[/color]
a = BilinearForm(fes)
a += -omega
[color=#1750eb]1j[/color]uv * ds([color=#067d17]â€śouterâ€ť[/color])
a.Assemble()

f = LinearForm(fes)
f += [color=#1750eb]1 [/color]* v * ds([color=#067d17]â€śscatâ€ť[/color])
f.Assemble();

gfu = GridFunction(fes, [color=#660099]name[/color]=[color=#067d17]â€śuâ€ť[/color])
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)[/code][/code]

I now want to use my computed eigen modes for the source geometry, that i compute as

```[code]a, b, fes = build_elasticity_system(mesh, steel, dirichlet) u = ngsolve.GridFunction(fes, multidim=num) lams = ngsolve.ArnoldiSolver(a.mat, b.mat, fes.FreeDofs(), u.vecs, shift)```[/code]

How can i achieve this ?

I further have some problems for the 3D version of all this, i do not really understand how to asign boundary conditions to specific parts of a herachical mesh. My code so far to generate the acoustic structure mesh is

[code][code]import ngsolve

from ngsolve import *

from netgen.stl import *
from netgen.csg import *
from netgen.meshing import *

geo2 = STLGeometry(â€śâ€¦/test_data/tuningfork.stlâ€ť)
m2 = geo2.GenerateMesh()

centerx = 0
centery = 0
centerz = 0

count = 0

for e in m2.Elements2D():
for v in e.vertices:
centerx = centerx + m2[v][0]
centery = centery + m2[v][1]
centerz = centerz + m2[v][2]
count = count + 1

centerx = centerx / count
centery = centery / count
centerz = centerz / count

R = 0
for e in m2.Elements2D():
for v in e.vertices:
x = -centerx + m2[v][0]
y = -centery + m2[v][1]
z = -centerz + m2[v][2]
tmp = sqrt(xx + yy + z*z)
if tmp > R:
R = tmp

print("R : " + str(R))

geo1 = CSGeometry()

# Set the mesh size on the sphere surface to 0.1

sphere = Sphere(Pnt(centerx, centery, centerz), 2.5*R)
m1 = geo1.GenerateMesh()

print(m2)

# create an empty mesh

mesh = netgen.meshing.Mesh()

# pmap1 maps point-numbers from old to new mesh

pmap1 = {}
for e in m1.Elements2D():
for v in e.vertices:
if v not in pmap1:

# we have to map point-numbers:

for e in m1.Elements2D():
mesh.Add(Element2D(fd_outside, [pmap1[v] for v in e.vertices]))

# same for the second mesh:

pmap2 = {}
for e in m2.Elements2D():
for v in e.vertices:
if v not in pmap2:

n = len(pmap2)

for e in m2.Elements2D():
mesh.Add(Element2D(fd_inside, [pmap2[v] for v in e.vertices]))

mesh.GenerateVolumeMesh(maxh=10.0)
mesh.Save(â€śacousticmesh.volâ€ť)

mesh = ngsolve.Mesh(mesh)[/code][/code] So i would be very happy if some one could give me some hints to put this all together

hm â€¦ my code is not visible in firefox â€¦

Ok, i have played around with some examples where the wave equation is used. I finally end up with a code that have the boundaries of a rectangle as source for the wave equation.

from ngsolve import *
from netgen.geom2d import SplineGeometry

# Geometry

geo = SplineGeometry()
leftdomain=0, rightdomain=1, bc=â€śscatâ€ť)
mesh = Mesh(geo.GenerateMesh(maxh=0.1))

fes = H1(mesh, order=5, complex=True)
u, v = fes.TnT()

# Wavenumber & source

omega = 15
pulse = 5e4exp(-(40**2)((x-0.5)(x-0.5) + (y-0.5)(y-0.5)))

# Forms

a = BilinearForm(fes)
a += -omega
1juv * ds(â€śouterâ€ť)
a.Assemble()

f = LinearForm(fes)
f += 1 * v * ds(â€śscatâ€ť)
f.Assemble();

gfu = GridFunction(fes, name=â€śuâ€ť)
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)

I now want to use my computed eigen modes for the source geometry, that i compute as

a, b, fes = build_elasticity_system(mesh, steel, dirichlet)
u = ngsolve.GridFunction(fes, multidim=num)
lams = ngsolve.ArnoldiSolver(a.mat, b.mat, fes.FreeDofs(), u.vecs, shift)

How can i achieve this ?

I further have some problems for the 3D version of all this, i do not really understand how to asign boundary conditions to specific parts of a herachical mesh. My code so far to generate the acoustic structure mesh is

import ngsolve

from ngsolve import *

from netgen.stl import *
from netgen.csg import *
from netgen.meshing import *

geo2 = STLGeometry(â€śâ€¦/test_data/tuningfork.stlâ€ť)
m2 = geo2.GenerateMesh()

centerx = 0
centery = 0
centerz = 0

count = 0

for e in m2.Elements2D():
for v in e.vertices:
centerx = centerx + m2[v][0]
centery = centery + m2[v][1]
centerz = centerz + m2[v][2]
count = count + 1

centerx = centerx / count
centery = centery / count
centerz = centerz / count

R = 0
for e in m2.Elements2D():
for v in e.vertices:
x = -centerx + m2[v][0]
y = -centery + m2[v][1]
z = -centerz + m2[v][2]
tmp = sqrt(xx + yy + z*z)
if tmp > R:
R = tmp

print("R : " + str(R))

geo1 = CSGeometry()

# Set the mesh size on the sphere surface to 0.1

sphere = Sphere(Pnt(centerx, centery, centerz), 2.5*R)
m1 = geo1.GenerateMesh()

print(m2)

# create an empty mesh

mesh = netgen.meshing.Mesh()

# pmap1 maps point-numbers from old to new mesh

pmap1 = {}
for e in m1.Elements2D():
for v in e.vertices:
if v not in pmap1:

# we have to map point-numbers:

for e in m1.Elements2D():
mesh.Add(Element2D(fd_outside, [pmap1[v] for v in e.vertices]))

# same for the second mesh:

pmap2 = {}
for e in m2.Elements2D():
for v in e.vertices:
if v not in pmap2:

n = len(pmap2)

for e in m2.Elements2D():
mesh.Add(Element2D(fd_inside, [pmap2[v] for v in e.vertices]))

mesh.GenerateVolumeMesh(maxh=10.0)
mesh.Save(â€śacousticmesh.volâ€ť)

mesh = ngsolve.Mesh(mesh)

So i would be very happy if some one could give me some hints to put this all together

looks like the [code] tag does not work