Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/examples/Python_FreeFem++/HelmholtzSolver.py
13582719ce94develop

Last Change
Raw File

HelmholtzSolver.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
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
PI = np.pi
testNo = 3
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))")
kappa = (12 + 1.j)**.5
solver = HF(V, kappa, forcingTerm = f, DirichletBoundary = "1,2,3,4")
uh = solver.solve()
plotter = HS(solver.V)
print(plotter.norm(uh, 12**.5))
plotter.plot(uh)
###########
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)
###########
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)
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh - uinc, name = 'u_tot')
###########
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 = utilities.listify(solver.liftDirichletData(), 2)
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)
plotter.plot(uh - uinc, name = 'u_tot')
c4science · Help