Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/hfengines/scipy/helmholtz_square_transmission_problem_engine.py
0c9a0bbd9f4edevelop

Last Change
Raw File

helmholtz_square_transmission_problem_engine.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import fenics as fen
from .helmholtz_problem_engine import HelmholtzProblemEngine
__all__ = ['HelmholtzSquareTransmissionProblemEngine']
class HelmholtzSquareTransmissionProblemEngine(HelmholtzProblemEngine):
"""
Solver for square transmission Helmholtz problems with parametric
wavenumber.
- \Delta u - omega^2 * n^2 * u = 0 in \Omega
u = 0 on \Gamma_D
with exact solution a transmitted plane wave.
"""
def __init__(self, nT:float, nB:float, kappa:float, theta:float, n:int):
super().__init__()
self.omega = kappa
mesh = fen.RectangleMesh(fen.Point(-np.pi/2, -np.pi/2),
fen.Point(np.pi/2, np.pi/2), n, n)
self.V = fen.FunctionSpace(mesh, "P", 3)
import sympy as sp
dx, dy = np.cos(theta), np.sin(theta)
Kx = kappa * nB * dx
Ky = kappa * (nT**2. - (nB * dx)**2. + 0.j)**.5
T = 2 * kappa * nB * dy / (Ky + kappa * nB * dy)
x, y = sp.symbols('x[0] x[1]', real=True)
uT = T * sp.exp(1.j * (Kx*x + Ky*y))
uB = ( sp.exp(1.j * kappa * nB * (dx*x + dy*y))
+ (T - 1) * sp.exp(1.j * kappa * nB * (dx*x - dy*y)))
uRe = fen.Expression('x[1]>=0 ? {} : {}'.format(
sp.printing.ccode(sp.re(uT)),
sp.printing.ccode(sp.re(uB))),
degree = 3)
uIm = fen.Expression('x[1]>=0 ? {} : {}'.format(
sp.printing.ccode(sp.im(uT)),
sp.printing.ccode(sp.im(uB))),
degree = 3)
self.refractionIndex = fen.Expression('x[1] >= 0 ? nT : nB',
nT = nT, nB = nB, degree = 0)
self.DirichletDatum = [uRe, uIm]
c4science · Help