Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60741663
HelmholtzSolver.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
Thu, May 2, 07:50
Size
6 KB
Mime Type
text/x-python
Expires
Sat, May 4, 07:50 (2 d)
Engine
blob
Format
Raw Data
Handle
17404605
Attached To
R6746 RationalROMPy
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
FenicsHelmholtzEngine
as
HF
from
context
import
FenicsHSEngine
as
HS
testNo
=
1
if
testNo
==
1
:
PI
=
np
.
pi
def
boundary
(
x
,
on_boundary
):
return
on_boundary
nu
=
12
**
.
5
theta
=
PI
/
3
x
,
y
=
sp
.
symbols
(
'x[0] x[1]'
,
real
=
True
)
wex
=
16
/
PI
**
4
*
x
*
y
*
(
x
-
PI
)
*
(
y
-
PI
)
phiex
=
nu
*
(
x
*
np
.
cos
(
theta
)
+
y
*
np
.
sin
(
theta
))
uex
=
wex
*
sp
.
exp
(
-
1.j
*
phiex
)
fex
=
-
uex
.
diff
(
x
,
2
)
-
uex
.
diff
(
y
,
2
)
-
nu
**
2
*
uex
nx
=
ny
=
40
mesh
=
fen
.
RectangleMesh
(
fen
.
Point
(
0
,
0
),
fen
.
Point
(
PI
,
PI
),
nx
,
ny
)
forcingTerm
=
[
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
re
(
fex
))),
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
im
(
fex
)))]
solver
=
HF
.
FenicsHelmholtzEngine
(
mesh
=
mesh
,
wavenumber
=
nu
+
1.j
,
forcingTerm
=
forcingTerm
,
FEDegree
=
3
,
\
DirichletBoundary
=
boundary
,
DirichletDatum
=
0
)
uh
=
solver
.
solve
()
plotter
=
HS
.
FenicsHSEngine
(
solver
.
V
)
print
(
plotter
.
norm
(
uh
,
np
.
real
(
nu
)))
plotter
.
plot
(
uh
)
###########
elif
testNo
==
2
:
def
boundary
(
x
,
on_boundary
):
return
on_boundary
PI
=
np
.
pi
n1
=
4
**.
5
n2
=
1
**.
5
kappa
=
3
theta
=
PI
*
70
/
180
d1
,
d2
=
np
.
cos
(
theta
),
np
.
sin
(
theta
)
K1
=
kappa
*
n1
*
d1
if
kappa
*
n2
>=
K1
:
K2
=
((
kappa
*
n2
)
**
2
-
K1
**
2
)
**.
5
else
:
K2
=
1.j
*
(
K1
**
2
-
(
kappa
*
n2
)
**
2
)
**.
5
R
=
(
kappa
*
n1
*
d2
-
K2
)
/
(
kappa
*
n1
*
d2
+
K2
)
T
=
R
+
1
x
,
y
=
sp
.
symbols
(
'x[0] x[1]'
,
real
=
True
)
uex1
=
T
*
sp
.
exp
(
1.j
*
(
K1
*
x
+
K2
*
y
))
uex2
=
sp
.
exp
(
1.j
*
kappa
*
n1
*
(
d1
*
x
+
d2
*
y
))
+
R
*
sp
.
exp
(
1.j
*
kappa
*
n1
*
(
d1
*
x
-
d2
*
y
))
# Exact solution
uexRe
=
fen
.
Expression
(
'x[1]>=0 ? {0} : {1}'
.
format
(
\
sp
.
printing
.
ccode
(
sp
.
re
(
uex1
)),
sp
.
printing
.
ccode
(
sp
.
re
(
uex2
))),
degree
=
4
)
uexIm
=
fen
.
Expression
(
'x[1]>=0 ? {0} : {1}'
.
format
(
\
sp
.
printing
.
ccode
(
sp
.
im
(
uex1
)),
sp
.
printing
.
ccode
(
sp
.
im
(
uex2
))),
degree
=
4
)
# refraction index
n2Re
=
fen
.
Expression
(
'x[1]<0 ? n1r : n2r'
,
n1r
=
n1
.
real
,
n2r
=
n2
.
real
,
degree
=
4
)
n2Im
=
fen
.
Expression
(
'x[1]<0 ? n1i : n2i'
,
n1i
=
n1
.
imag
,
n2i
=
n2
.
imag
,
degree
=
4
)
# Create mesh and define function space
nx
=
ny
=
50
mesh
=
fen
.
RectangleMesh
(
fen
.
Point
(
-
PI
/
2
,
-
PI
/
2
),
fen
.
Point
(
PI
/
2
,
PI
/
2
),
nx
,
ny
)
solver
=
HF
.
FenicsHelmholtzEngine
(
mesh
=
mesh
,
wavenumber
=
kappa
,
refractionIndex
=
(
n2Re
,
n2Im
),
\
forcingTerm
=
0
,
FEDegree
=
3
,
DirichletBoundary
=
boundary
,
\
DirichletDatum
=
(
uexRe
,
uexIm
))
uh
=
solver
.
solve
()
plotter
=
HS
.
FenicsHSEngine
(
solver
.
V
)
print
(
plotter
.
norm
(
uh
,
kappa
))
plotter
.
plot
(
uh
)
###########
elif
testNo
==
3
:
import
mshr
from
matplotlib
import
pyplot
as
plt
PI
=
np
.
pi
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
=
-
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
=
50
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
)
plt
.
jet
()
plt
.
figure
()
fen
.
plot
(
mesh
)
DirichletTerm
=
[
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
re
(
u0ex
))),
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
im
(
u0ex
)))]
solver
=
HF
.
FenicsHelmholtzScatteringEngine
(
mesh
=
mesh
,
wavenumber
=
kappa
,
forcingTerm
=
0
,
FEDegree
=
3
,
\
DirichletBoundary
=
Dboundary
,
RobinBoundary
=
'rest'
,
\
DirichletDatum
=
DirichletTerm
)
baseRe
,
baseIm
=
DirichletTerm
baseRe
=
fen
.
project
(
fen
.
Expression
(
baseRe
,
degree
=
4
),
solver
.
V
)
baseIm
=
fen
.
project
(
fen
.
Expression
(
baseIm
,
degree
=
4
),
solver
.
V
)
uinc
=
np
.
array
(
baseRe
.
vector
())
+
1.j
*
np
.
array
(
baseIm
.
vector
())
uh
=
solver
.
solve
()
plotter
=
HS
.
FenicsHSEngine
(
solver
.
V
)
print
(
plotter
.
norm
(
uh
,
kappa
))
print
(
plotter
.
norm
(
uh
-
uinc
,
kappa
))
plotter
.
plot
(
uh
)
plotter
.
plot
(
uh
-
uinc
)
###########
elif
testNo
==
4
:
import
mshr
from
matplotlib
import
pyplot
as
plt
PI
=
np
.
pi
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
=
-
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
=
40
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
)
plt
.
jet
()
plt
.
figure
()
fen
.
plot
(
mesh
)
DirichletTerm
=
[
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
re
(
u0ex
))),
sp
.
printing
.
ccode
(
sp
.
simplify
(
sp
.
im
(
u0ex
)))]
solver
=
HF
.
FenicsHelmholtzScatteringAugmentedEngine
(
mesh
=
mesh
,
wavenumber
=
kappa
,
forcingTerm
=
0
,
FEDegree
=
3
,
\
DirichletBoundary
=
Dboundary
,
RobinBoundary
=
'rest'
,
\
DirichletDatum
=
DirichletTerm
,
constraintType
=
'MASS'
)
baseRe
,
baseIm
=
DirichletTerm
baseRe
=
fen
.
project
(
fen
.
Expression
(
baseRe
,
degree
=
4
),
solver
.
V
)
baseIm
=
fen
.
project
(
fen
.
Expression
(
baseIm
,
degree
=
4
),
solver
.
V
)
uinc
=
np
.
array
(
baseRe
.
vector
())
+
1.j
*
np
.
array
(
baseIm
.
vector
())
uinc
=
np
.
concatenate
((
uinc
,
kappa
*
uinc
))
uh
=
solver
.
solve
()
plotter
=
HS
.
FenicsHSAugmentedEngine
(
solver
.
V
,
2
)
print
(
plotter
.
norm
(
uh
,
kappa
))
print
(
plotter
.
norm
(
uh
-
uinc
,
kappa
))
plotter
.
plot
(
uh
)
plotter
.
plot
(
uh
-
uinc
)
Event Timeline
Log In to Comment