Page MenuHomec4science
Files F73218595

HelmholtzSolver.py
No OneTemporary
Actions

File Metadata

Created
Fri, Jul 19, 17:24

HelmholtzSolver.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
import fenics as fen
import numpy as np
import sympy as sp
from context import utilities
from context import FenicsHelmholtzEngine as HF
from context import FenicsHelmholtzScatteringEngine as HFS
from context import FenicsHelmholtzScatteringAugmentedEngine as HFSA
from context import FenicsHSEngine as HS
from context import FenicsHSAugmentedEngine as HSA
testNo = 3
if testNo == 1:
from FEniCS_snippets import SquareHomogeneousBubble
boundary, mesh, forcingTerm = SquareHomogeneousBubble(kappa = 12 ** .5,
theta = np.pi / 3,
n = 40)
solver = HF(mesh = mesh, wavenumber = (12 + 1.j)**.5,
forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = boundary, DirichletDatum = 0)
uh = solver.solve()
plotter = HS(solver.V)
plotter.plotmesh(save = True)
print(plotter.norm(uh))
plotter.plot(uh, save = True)
###########
elif testNo == 2:
from FEniCS_snippets import SquareTransmissionDirichlet
boundary, mesh, n, u0 = SquareTransmissionDirichlet(nT = 1, nB = 2,
theta = np.pi * 20 / 180,
kappa = 4., n = 40)
solver = HF(mesh = mesh, wavenumber = 4., refractionIndex = n,
FEDegree = 3, DirichletBoundary = boundary,
DirichletDatum = u0)
uh = solver.solve()
plotter = HS(solver.V)
print(plotter.norm(uh))
plotter.plot(uh)
###########
elif testNo in [3, 4]:
import mshr
R = 5
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
A = 10
kappa = 12**.5
theta = - np.pi * 90 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - A * sp.exp(1.j * phiex)
npoints = 30
scatterer = mshr.Polygon([fen.Point(-1, -.5),
fen.Point(1, -.5),
fen.Point(1, .5),
fen.Point(.8, .5),
fen.Point(.8, -.3),
fen.Point(-.8, -.3),
fen.Point(-.8, .5),
fen.Point(-1, .5),
])
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer,
npoints)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))),
sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
if testNo == 3:
solver = HFS(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm)
plotter = HS(solver.V)
uinc = solver.liftDirichletData()
uh = solver.solve()
else:
solver = HFSA(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
constraintType = 'MASS')
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