More precisely i want to know if there is any pitfall for the implementation of something like seen in this elmer fem video.
[video]- YouTube
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()
geo.AddCircle(([color=#1750eb]0.0[/color], [color=#1750eb]0.0[/color]), [color=#1750eb]0.8[/color], [color=#660099]bc[/color]=[color=#067d17]“outer”[/color])
geo.AddRectangle((-[color=#1750eb]0.05[/color], -[color=#1750eb]0.025[/color]), ([color=#1750eb]0.025[/color], [color=#1750eb]0.275[/color]),
[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 += grad(u)grad(v)dx - omega[color=#1750eb]2[/color]uvdx
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)
geo1.Add(sphere)
m1 = geo1.GenerateMesh()
print(m2)
create an empty mesh
mesh = netgen.meshing.Mesh()
a face-descriptor stores properties associated with a set of surface elements
bc … boundary condition marker,
domin/domout … domain-number in front/back of surface elements (0 = void),
surfnr … number of the surface described by the face-descriptor
fd_outside = mesh.Add(FaceDescriptor(bc=1,domin=1,surfnr=1))
fd_inside = mesh.Add(FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))
copy all boundary points from first mesh to new 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:
pmap1[v] = mesh.Add (m1[v])
copy surface elements from first mesh to new mesh
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:
pmap2[v] = mesh.Add(m2[v])
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 -repost for check if i had made a formating mistake
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()
geo.AddCircle((0.0, 0.0), 0.8, bc=“outer”)
geo.AddRectangle((-0.05, -0.025), (0.025, 0.275),
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 += grad(u)grad(v)dx - omega**2uvdx
a += -omega1juv * 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)
geo1.Add(sphere)
m1 = geo1.GenerateMesh()
print(m2)
create an empty mesh
mesh = netgen.meshing.Mesh()
a face-descriptor stores properties associated with a set of surface elements
bc … boundary condition marker,
domin/domout … domain-number in front/back of surface elements (0 = void),
surfnr … number of the surface described by the face-descriptor
fd_outside = mesh.Add(FaceDescriptor(bc=1,domin=1,surfnr=1))
fd_inside = mesh.Add(FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))
copy all boundary points from first mesh to new 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:
pmap1[v] = mesh.Add (m1[v])
copy surface elements from first mesh to new mesh
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:
pmap2[v] = mesh.Add(m2[v])
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