Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FEniCS/HelmholtzSolver.py
6e6c171cf128master

Last Change
Raw File

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')
c4science · Help