Page MenuHomec4science
Files F62203287

HelmholtzSolver_FF.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 11, 14:43

HelmholtzSolver_FF.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
import numpy as np
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
PI = np.pi
testNo = 1
if testNo == 1:
V = ("int n = 40;\n"
"real nu = 12^.5, theta = pi / 3;\n"
"mesh Th = square(n, n, [pi * x, pi * y]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);")
f = ("16 / pi^4 * exp(-1i * nu*(cos(theta) * x + sin(theta) * y)) * ( 2i *"
" nu *cos(theta) * (2 * x * y^2 - 2 * pi * x * y - pi * y^2 + pi^2 * "
"y) + 2i * nu * sin(theta) * (2 *x^2 * y - pi * x^2 - 2 * pi * x * y "
"+ pi^2 * x) - (2 * y^2 - 2 * pi * y + 2 * x^2 - 2 * pi * x))")
solver = HF(V, 12**.5 + 1.j, forcingTerm = f, DirichletBoundary = "1,2,3,4")
uh = solver.solve()
plotter = HS(solver.V)
print(plotter.norm(uh, 12**.5))
plotter.plot(uh, figspecs = "fill = 1, value = 1", save = True)
###########
elif testNo == 2:
kappa = 3
theta = np.pi * 70 / 180
k1 = kappa * 2 * np.cos(theta)
if kappa > k1:
k2 = (kappa**2. - k1**2.)**.5
else:
k2 = 1.j * (k1**2. - kappa**2.)**.5
R = (kappa * 2 * np.sin(theta) - k2) / (kappa * 2 * np.sin(theta) + k2)
T = R + 1
V = ("int n = 50;\n"
"real kappa = {}, theta = {};\n"
"mesh Th = square(n, n, [pi * (x - .5), pi * (y - .5)]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);").format(kappa, theta)
n = "(y >= 0 ? 1 : 2)"
u0 = ("(y >= 0 ? {0}*exp(1.i*({1}*x+{2}*y)) : exp(2.i*{3}*({4}*x+{5}*y)) + "
"{6}*exp(2.i*{3}*({4}*x-{5}*y)))".format(T, k1, k2, kappa,
np.cos(theta), np.sin(theta), R))
solver = HF(V, kappa, refractionIndex = n, DirichletBoundary = "1,2,3,4", DirichletDatum = u0)
uh = solver.solve()
plotter = HS(solver.V)
print(plotter.norm(uh, kappa))
plotter.plot(uh, figspecs = "fill = 1, value = 1")
###########
elif testNo == 3:
A = 10
kappa = 12**.5
theta = - 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)
solver = HFS(V, kappa, DirichletBoundary = "1", DirichletDatum = u0, RobinBoundary = "2")
uinc = solver.liftDirichletData()
uh = solver.solve()
plotter = HS(solver.V)
plotter.plotmesh()
print(plotter.norm(uh, kappa))
plotter.plot(uh, figspecs = 'fill = 1, value = 1')
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh - uinc, name = 'u_tot', figspecs = 'fill = 1, value = 1')
###########
elif testNo == 4:
A = 10
kappa = 12**.5
theta = - PI * 90 / 180
V = ("int n = 5, 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)
solver = HFSA(V, kappa, DirichletBoundary = "1", DirichletDatum = u0, RobinBoundary = "2")
uinc = solver.liftDirichletData()
uinc[1] = kappa * uinc[0]
uh = solver.solve()
plotter = HSA(solver.V, 2)
plotter.plotmesh()
print(plotter.norm(uh, kappa))
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh, figspecs = 'fill = 1, value = 1')
plotter.plot(uh - uinc, name = 'u_tot', figspecs = 'fill = 1, value = 1')

Event Timeline

c4science · Help