Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82719012
helmholtz_square_transmission_problem_engine.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 13, 01:42
Size
3 KB
Mime Type
text/x-python
Expires
Sun, Sep 15, 01:42 (2 d)
Engine
blob
Format
Raw Data
Handle
20739578
Attached To
R6746 RationalROMPy
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
import
ufl
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
,
degree_threshold
:
int
=
np
.
inf
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
super
()
.
__init__
(
degree_threshold
=
degree_threshold
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
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
)
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
=
fen
.
SpatialCoordinate
(
mesh
)[:]
TR
,
TI
=
np
.
real
(
T
),
np
.
imag
(
T
)
if
np
.
isclose
(
np
.
imag
(
Ky
),
0.
):
u0RT
=
(
TR
*
fen
.
cos
(
Kx
*
x
+
np
.
real
(
Ky
)
*
y
)
-
TI
*
fen
.
sin
(
Kx
*
x
+
np
.
real
(
Ky
)
*
y
))
u0IT
=
(
TR
*
fen
.
sin
(
Kx
*
x
+
np
.
real
(
Ky
)
*
y
)
+
TI
*
fen
.
cos
(
Kx
*
x
+
np
.
real
(
Ky
)
*
y
))
else
:
u0RT
=
fen
.
exp
(
-
np
.
imag
(
Ky
)
*
y
)
*
(
TR
*
fen
.
cos
(
Kx
*
x
)
-
TI
*
fen
.
sin
(
Kx
*
x
))
u0IT
=
fen
.
exp
(
-
np
.
imag
(
Ky
)
*
y
)
*
(
TR
*
fen
.
sin
(
Kx
*
x
)
+
TI
*
fen
.
cos
(
Kx
*
x
))
u0RB
=
(
fen
.
cos
(
kappa
*
nB
*
(
dx
*
x
+
dy
*
y
))
+
(
TR
-
1
)
*
fen
.
cos
(
kappa
*
nB
*
(
dx
*
x
-
dy
*
y
))
-
TI
*
fen
.
sin
(
kappa
*
nB
*
(
dx
*
x
-
dy
*
y
)))
u0IB
=
(
fen
.
sin
(
kappa
*
nB
*
(
dx
*
x
+
dy
*
y
))
+
(
TR
-
1
)
*
fen
.
sin
(
kappa
*
nB
*
(
dx
*
x
-
dy
*
y
))
+
TI
*
fen
.
cos
(
kappa
*
nB
*
(
dx
*
x
-
dy
*
y
)))
u0R
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
0.
),
u0RT
,
u0RB
)
u0I
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
0.
),
u0IT
,
u0IB
)
self
.
refractionIndex
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
0.
),
fen
.
Constant
(
nT
),
fen
.
Constant
(
nB
))
self
.
DirichletDatum
=
[
u0R
,
u0I
]
Event Timeline
Log In to Comment