Page MenuHomec4science
Files F90032034

HelmholtzSolver.py
No OneTemporary
Actions

File Metadata

Created
Mon, Oct 28, 16:11

HelmholtzSolver.py
View Options

#!/usr/bin/env python3
import numpy as np
from context import utilities
from context import FreeFemHelmholtzEngine as HF
from context import FreeFemHelmholtzScatteringEngine as HFS
from context import FreeFemHelmholtzScatteringAugmentedEngine as HFSA
from context import FreeFemHSEngine as HS
from context import FreeFemHSAugmentedEngine as HSA
testNo = 2
if testNo == 1:
from FreeFem_snippets import SquareHomogeneousBubble
bdr, V, f = SquareHomogeneousBubble(kappa = 12 ** .5,
theta = np.pi / 3,
n = 40)
solver = HF(V, (12 + 1.j)**.5, forcingTerm = f, DirichletBoundary = bdr)
uh = solver.solve()
plotter = HS(solver.V)
plotter.plotmesh(save = True)
print(plotter.norm(uh))
plotter.plot(uh, save = True)
###########
elif testNo == 2:
from FreeFem_snippets import SquareTransmissionDirichlet
bdr, V, n, u0 = SquareTransmissionDirichlet(nT = 1, nB = 2,
theta = np.pi * 20 / 180,
kappa = 4., n = 40)
solver = HF(V, 4., refractionIndex = n, DirichletBoundary = bdr,
DirichletDatum = u0)
uh = solver.solve()
plotter = HS(solver.V)
print(plotter.norm(uh, 4.))
plotter.plot(uh)
###########
elif testNo in [3, 4]:
A = 10
kappa = 12**.5
theta = - np.pi * 90 / 180
V = ("int n = 10, R = 5;\n"
"real kappa = {}, theta = {}, A = {};\n"
"border S1(t = 0, 1){{x = 2*t-1.; y = -.5; label = 1;}}\n"
"border S2(t = 0, 1){{x = 1.; y = -.5+t; label = 1;}}\n"
"border S3(t = 0, 1){{x = 1-.2*t; y = .5; label = 1;}}\n"
"border S4(t = 0, 1){{x = .8; y = .5-.8*t; label = 1;}}\n"
"border S5(t = 0, 1){{x = .8-1.6*t; y = -.3; label = 1;}}\n"
"border S6(t = 0, 1){{x = -.8; y = -.3+.8*t; label = 1;}}\n"
"border S7(t = 0, 1){{x = -.8-.2*t; y = .5; label = 1;}}\n"
"border S8(t = 0, 1){{x = -1.; y = .5-t; label = 1;}}\n"
"border Inf(t = 0, 1){{x = R * cos(2 * pi * t); "
"y = R * sin(2 * pi * t); label = 2;}}\n"
"mesh Th = buildmesh(S1(-4*n) + S2(-3*n) + S3(-n) + S4(-2*n) + "
"S5(-4*n) + S6(-2*n) + S7(-n) + S8(-3*n) + Inf(10*n));\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);").format(kappa, theta, A)
u0 = "- {0} * exp(1.i * {1} * (cos({2}) * x + sin({2}) * y))".format(A,
kappa, theta)
if testNo == 3:
solver = HFS(V, kappa, DirichletBoundary = "1", DirichletDatum = u0,
RobinBoundary = "2")
plotter = HS(solver.V)
uinc = solver.liftDirichletData()
else:
solver = HFSA(V, kappa, DirichletBoundary = "1", DirichletDatum = u0,
RobinBoundary = "2")
plotter = HSA(solver.V, 2)
uinc = utilities.listify(solver.liftDirichletData(), 2)
uinc[1] = kappa * uinc[0]
uh = solver.solve()
plotter.plotmesh()
print(plotter.norm(uh, kappa))
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh)
plotter.plot(uh - uinc, name = 'u_tot')

Event Timeline

c4science · Help