Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88604823
anisotropic_square_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
Sat, Oct 19, 17:15
Size
3 KB
Mime Type
text/x-python
Expires
Mon, Oct 21, 17:15 (2 d)
Engine
blob
Format
Raw Data
Handle
21796907
Attached To
R6746 RationalROMPy
anisotropic_square_engine.py
View Options
import
numpy
as
np
import
fenics
as
fen
import
mshr
import
ufl
from
rrompy.hfengines.linear_problem
import
HelmholtzProblemEngine
from
rrompy.solver.fenics
import
fenONE
,
fenZERO
,
fenics2Sparse
class
AnisotropicSquareEngine
(
HelmholtzProblemEngine
):
def
__init__
(
self
,
k2
:
float
,
L2
:
float
,
n
:
int
):
super
()
.
__init__
(
mu0
=
[
k2
,
L2
])
self
.
_affinePoly
=
True
self
.
nAs
=
3
self
.
npar
=
2
self
.
rescalingExp
=
[
1.
,
1.
]
self
.
V
=
fen
.
FunctionSpace
(
fen
.
UnitSquareMesh
(
n
,
n
),
"P"
,
1
)
eps
=
1e-6
self
.
DirichletBoundary
=
lambda
x
,
on_boundary
:
(
on_boundary
and
x
[
1
]
<
eps
)
self
.
NeumannBoundary
=
"REST"
x
,
y
=
fen
.
SpatialCoordinate
(
self
.
V
.
mesh
())[:]
self
.
NeumannDatum
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
1.
-
eps
),
fen
.
cos
(
np
.
pi
*
x
),
fenZERO
)
self
.
forcingTerm
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
.
5
),
fenONE
,
fenZERO
)
*
(
5
*
ufl
.
conditional
(
ufl
.
lt
(
x
,
.
1
),
fenONE
,
fenZERO
)
-
5
*
ufl
.
conditional
(
ufl
.
And
(
ufl
.
gt
(
x
,
.
2
),
ufl
.
lt
(
x
,
.
3
)),
fenONE
,
fenZERO
)
+
10
*
ufl
.
conditional
(
ufl
.
And
(
ufl
.
gt
(
x
,
.
45
),
ufl
.
lt
(
x
,
.
55
)),
fenONE
,
fenZERO
)
-
5
*
ufl
.
conditional
(
ufl
.
And
(
ufl
.
gt
(
x
,
.
7
),
ufl
.
lt
(
x
,
.
8
)),
fenONE
,
fenZERO
)
+
5
*
ufl
.
conditional
(
ufl
.
gt
(
x
,
.
9
),
fenONE
,
fenZERO
))
def
buildA
(
self
):
"""Build terms of operator of linear system."""
if
self
.
thAs
[
0
]
is
None
:
self
.
thAs
=
self
.
getMonomialWeights
(
self
.
nAs
)
if
self
.
As
[
0
]
is
None
:
self
.
autoSetDS
()
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
a0
=
fen
.
dot
(
self
.
u
.
dx
(
0
),
self
.
v
.
dx
(
0
))
*
fen
.
dx
self
.
As
[
0
]
=
fenics2Sparse
(
a0
,
{},
DirichletBC0
,
1
)
if
self
.
As
[
1
]
is
None
:
self
.
autoSetDS
()
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
a1
=
-
fen
.
dot
(
self
.
u
,
self
.
v
)
*
fen
.
dx
self
.
As
[
1
]
=
fenics2Sparse
(
a1
,
{},
DirichletBC0
,
0
)
if
self
.
As
[
2
]
is
None
:
self
.
autoSetDS
()
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
a2
=
fen
.
dot
(
self
.
u
.
dx
(
1
),
self
.
v
.
dx
(
1
))
*
fen
.
dx
self
.
As
[
2
]
=
fenics2Sparse
(
a2
,
{},
DirichletBC0
,
0
)
def
AnisotropicSquareEnginePoles
(
L2
:
float
,
k2min
:
float
,
k2max
:
float
):
poles
=
[]
for
alpha
in
np
.
arange
(
np
.
ceil
((
k2max
)
**
.
5
/
np
.
pi
)):
p
=
(
np
.
pi
*
alpha
)
**
2.
pkmin
=
np
.
ceil
(
max
(
0.
,
(
k2min
-
p
)
*
4
/
L2
)
**
.
5
/
np
.
pi
)
pkmin
+=
1
-
(
pkmin
%
2
)
pkmax
=
np
.
floor
(
max
(
0.
,
(
k2max
-
p
)
*
4
/
L2
)
**
.
5
/
np
.
pi
)
for
beta
in
np
.
arange
(
pkmin
,
pkmax
+
1
,
2
):
poles
+=
[
p
+
L2
*
(
np
.
pi
*
beta
/
2.
)
**
2.
]
return
np
.
unique
(
poles
)
Event Timeline
Log In to Comment