Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79433441
helmholtz_base_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
Mon, Aug 26, 09:20
Size
12 KB
Mime Type
text/x-python
Expires
Wed, Aug 28, 09:20 (2 d)
Engine
blob
Format
Raw Data
Handle
20186023
Attached To
R6746 RationalROMPy
helmholtz_base_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
scipy.sparse
as
scsp
import
fenics
as
fen
from
rrompy.hfengines.fenics.generic_problem_engine
import
GenericProblemEngine
from
rrompy.utilities.base.types
import
Np1D
,
Np2D
from
rrompy.utilities.base.fenics
import
fenZERO
,
fenONE
,
bdrTrue
,
bdrFalse
__all__
=
[
'HelmholtzBaseProblemEngine'
]
class
HelmholtzBaseProblemEngine
(
GenericProblemEngine
):
"""
ABSTRACT
Solver for generic Helmholtz problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
omega: Value of omega.
diffusivity: Value of a.
refractionIndex: Value of n.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
b0: Numpy array representation of b0.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
"""
Aterms
=
2
omega
=
1.
_diffusivity
=
[
fenONE
,
fenZERO
]
_refractionIndex
=
[
fenONE
,
fenZERO
]
_forcingTerm
=
[
fenZERO
,
fenZERO
]
_DirichletDatum
=
[
fenZERO
,
fenZERO
]
_NeumannDatum
=
[
fenZERO
,
fenZERO
]
_RobinDatumG
=
[
fenZERO
,
fenZERO
]
def
__init__
(
self
):
self
.
V
=
fen
.
FunctionSpace
(
fen
.
UnitSquareMesh
(
10
,
10
),
"P"
,
1
)
self
.
DirichletBoundary
=
"ALL"
@property
def
V
(
self
):
"""Value of V."""
return
self
.
_V
@V.setter
def
V
(
self
,
V
):
if
hasattr
(
self
,
"A0"
):
del
self
.
A0
if
hasattr
(
self
,
"A1"
):
del
self
.
A1
if
hasattr
(
self
,
"b0"
):
del
self
.
b0
if
not
type
(
V
)
.
__name__
==
'FunctionSpace'
:
raise
Exception
(
"V type not recognized."
)
self
.
_V
=
V
self
.
u
=
fen
.
TrialFunction
(
V
)
self
.
v
=
fen
.
TestFunction
(
V
)
self
.
dsToBeSet
=
True
@property
def
diffusivity
(
self
):
"""Value of a."""
return
self
.
_diffusivity
@diffusivity.setter
def
diffusivity
(
self
,
diffusivity
):
if
hasattr
(
self
,
"A0"
):
del
self
.
A0
if
not
isinstance
(
diffusivity
,
(
list
,
tuple
,)):
diffusivity
=
[
diffusivity
,
fenZERO
]
self
.
_diffusivity
=
diffusivity
@property
def
refractionIndex
(
self
):
"""Value of n."""
return
self
.
_refractionIndex
@refractionIndex.setter
def
refractionIndex
(
self
,
refractionIndex
):
if
hasattr
(
self
,
"A1"
):
del
self
.
A1
if
not
isinstance
(
refractionIndex
,
(
list
,
tuple
,)):
refractionIndex
=
[
refractionIndex
,
fenZERO
]
self
.
_refractionIndex
=
refractionIndex
@property
def
forcingTerm
(
self
):
"""Value of f."""
return
self
.
_forcingTerm
@forcingTerm.setter
def
forcingTerm
(
self
,
forcingTerm
):
if
hasattr
(
self
,
"b0"
):
del
self
.
b0
if
not
isinstance
(
forcingTerm
,
(
list
,
tuple
,)):
forcingTerm
=
[
forcingTerm
,
fenZERO
]
self
.
_forcingTerm
=
forcingTerm
@property
def
DirichletDatum
(
self
):
"""
Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
DirichletBCIm.
"""
return
self
.
_DirichletDatum
@DirichletDatum.setter
def
DirichletDatum
(
self
,
DirichletDatum
):
if
hasattr
(
self
,
"b0"
):
del
self
.
b0
if
not
isinstance
(
DirichletDatum
,
(
list
,
tuple
,)):
DirichletDatum
=
[
DirichletDatum
,
fenZERO
]
self
.
_DirichletDatum
=
DirichletDatum
@property
def
NeumannDatum
(
self
):
"""Value of g1."""
return
self
.
_NeumannDatum
@NeumannDatum.setter
def
NeumannDatum
(
self
,
NeumannDatum
):
if
hasattr
(
self
,
"b0"
):
del
self
.
b0
if
not
isinstance
(
NeumannDatum
,
(
list
,
tuple
,)):
NeumannDatum
=
[
NeumannDatum
,
fenZERO
]
self
.
_NeumannDatum
=
NeumannDatum
@property
def
RobinDatumG
(
self
):
"""Value of g2."""
return
self
.
_RobinDatumG
@RobinDatumG.setter
def
RobinDatumG
(
self
,
RobinDatumG
):
if
hasattr
(
self
,
"b0"
):
del
self
.
b0
if
not
isinstance
(
RobinDatumG
,
(
list
,
tuple
,)):
RobinDatumG
=
[
RobinDatumG
,
fenZERO
]
self
.
_RobinDatumG
=
RobinDatumG
@property
def
DirichletBoundary
(
self
):
"""Function handle to DirichletBoundary."""
return
self
.
_DirichletBoundary
@DirichletBoundary.setter
def
DirichletBoundary
(
self
,
DirichletBoundary
):
if
hasattr
(
self
,
"A0"
):
del
self
.
A0
if
hasattr
(
self
,
"A1"
):
del
self
.
A1
if
isinstance
(
DirichletBoundary
,
(
str
,)):
if
DirichletBoundary
.
upper
()
==
"NONE"
:
self
.
_DirichletBoundary
=
bdrFalse
self
.
DREST
=
0
elif
DirichletBoundary
.
upper
()
==
"ALL"
:
self
.
_DirichletBoundary
=
bdrTrue
self
.
NeumannBoundary
=
"NONE"
self
.
RobinBoundary
=
"NONE"
self
.
DREST
=
0
elif
DirichletBoundary
.
upper
()
==
"REST"
:
if
self
.
NREST
+
self
.
RREST
>
0
:
raise
Exception
(
"Only 1 'REST' wildcard can be specified."
)
self
.
_DirichletBoundary
=
lambda
x
,
on_boundary
:
(
on_boundary
and
not
self
.
NeumannBoundary
(
x
,
on_boundary
)
and
not
self
.
RobinBoundary
(
x
,
on_boundary
))
self
.
DREST
=
1
else
:
raise
Exception
(
"DirichletBoundary label not recognized."
)
elif
callable
(
DirichletBoundary
):
self
.
_DirichletBoundary
=
DirichletBoundary
self
.
DREST
=
0
else
:
raise
Exception
(
"DirichletBoundary type not recognized."
)
@property
def
NeumannBoundary
(
self
):
"""Function handle to NeumannBoundary."""
return
self
.
_NeumannBoundary
@NeumannBoundary.setter
def
NeumannBoundary
(
self
,
NeumannBoundary
):
if
hasattr
(
self
,
"A1"
):
del
self
.
A1
self
.
dsToBeSet
=
True
if
isinstance
(
NeumannBoundary
,
(
str
,)):
if
NeumannBoundary
.
upper
()
==
"NONE"
:
self
.
_NeumannBoundary
=
bdrFalse
self
.
NREST
=
0
elif
NeumannBoundary
.
upper
()
==
"ALL"
:
self
.
_NeumannBoundary
=
bdrTrue
self
.
DirichletBoundary
=
"NONE"
self
.
RobinBoundary
=
"NONE"
self
.
NREST
=
0
elif
NeumannBoundary
.
upper
()
==
"REST"
:
if
self
.
DREST
+
self
.
RREST
>
0
:
raise
Exception
(
"Only 1 'REST' wildcard can be specified."
)
self
.
_NeumannBoundary
=
lambda
x
,
on_boundary
:
(
on_boundary
and
not
self
.
DirichletBoundary
(
x
,
on_boundary
)
and
not
self
.
RobinBoundary
(
x
,
on_boundary
))
self
.
NREST
=
1
else
:
raise
Exception
(
"NeumannBoundary label not recognized."
)
elif
callable
(
NeumannBoundary
):
self
.
_NeumannBoundary
=
NeumannBoundary
self
.
NREST
=
0
else
:
raise
Exception
(
"DirichletBoundary type not recognized."
)
@property
def
RobinBoundary
(
self
):
"""Function handle to RobinBoundary."""
return
self
.
_RobinBoundary
@RobinBoundary.setter
def
RobinBoundary
(
self
,
RobinBoundary
):
if
hasattr
(
self
,
"A0"
):
del
self
.
A0
if
hasattr
(
self
,
"A1"
):
del
self
.
A1
self
.
dsToBeSet
=
True
if
isinstance
(
RobinBoundary
,
(
str
,)):
if
RobinBoundary
.
upper
()
==
"NONE"
:
self
.
_RobinBoundary
=
bdrFalse
self
.
RREST
=
0
elif
RobinBoundary
.
upper
()
==
"ALL"
:
self
.
_RobinBoundary
=
bdrTrue
self
.
DirichletBoundary
=
"NONE"
self
.
NeumannBoundary
=
"NONE"
self
.
RREST
=
0
elif
RobinBoundary
.
upper
()
==
"REST"
:
if
self
.
DREST
+
self
.
NREST
>
0
:
raise
Exception
(
"Only 1 'REST' wildcard can be specified."
)
self
.
_RobinBoundary
=
lambda
x
,
on_boundary
:
(
on_boundary
and
not
self
.
DirichletBoundary
(
x
,
on_boundary
)
and
not
self
.
NeumannBoundary
(
x
,
on_boundary
))
self
.
RREST
=
1
else
:
raise
Exception
(
"RobinBoundary label not recognized."
)
return
elif
callable
(
RobinBoundary
):
self
.
_RobinBoundary
=
RobinBoundary
self
.
RREST
=
0
else
:
raise
Exception
(
"RobinBoundary type not recognized."
)
def
autoSetDS
(
self
):
"""Set FEniCS boundary measure based on boundary function handles."""
if
self
.
dsToBeSet
:
NB
=
self
.
NeumannBoundary
RB
=
self
.
RobinBoundary
class
NBoundary
(
fen
.
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
NB
(
x
,
on_boundary
)
class
RBoundary
(
fen
.
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
RB
(
x
,
on_boundary
)
boundary_markers
=
fen
.
MeshFunction
(
"size_t"
,
self
.
V
.
mesh
(),
self
.
V
.
mesh
()
.
topology
()
.
dim
()
-
1
)
NBoundary
()
.
mark
(
boundary_markers
,
0
)
RBoundary
()
.
mark
(
boundary_markers
,
1
)
self
.
ds
=
fen
.
Measure
(
"ds"
,
domain
=
self
.
V
.
mesh
(),
subdomain_data
=
boundary_markers
)
self
.
dsToBeSet
=
False
def
energyNormMatrix
(
self
)
->
Np2D
:
"""Sparse matrix (in CSR format) representative of scalar product."""
normMatFen
=
fen
.
assemble
((
fen
.
dot
(
fen
.
grad
(
self
.
u
),
fen
.
grad
(
self
.
v
))
+
np
.
abs
(
self
.
omega
)
**
2
*
fen
.
dot
(
self
.
u
,
self
.
v
))
*
fen
.
dx
)
normMat
=
fen
.
as_backend_type
(
normMatFen
)
.
mat
()
return
scsp
.
csr_matrix
(
normMat
.
getValuesCSR
()[::
-
1
],
shape
=
normMat
.
size
)
def
liftDirichletData
(
self
)
->
Np1D
:
"""Lift Dirichlet datum."""
solLRe
=
fen
.
interpolate
(
self
.
DirichletDatum
[
0
],
self
.
V
)
solLIm
=
fen
.
interpolate
(
self
.
DirichletDatum
[
1
],
self
.
V
)
return
np
.
array
(
solLRe
.
vector
())
+
1.j
*
np
.
array
(
solLIm
.
vector
())
def
assembleb
(
self
):
"""Assemble RHS blocks of linear system."""
if
not
hasattr
(
self
,
"b0"
):
fRe
,
fIm
=
self
.
forcingTerm
g1Re
,
g1Im
=
self
.
NeumannDatum
g2Re
,
g2Im
=
self
.
RobinDatumG
L0Re
=
(
fen
.
dot
(
fRe
,
self
.
v
)
*
fen
.
dx
+
fen
.
dot
(
g1Re
,
self
.
v
)
*
self
.
ds
(
0
)
+
fen
.
dot
(
g2Re
,
self
.
v
)
*
self
.
ds
(
1
))
L0Im
=
(
fen
.
dot
(
fIm
,
self
.
v
)
*
fen
.
dx
+
fen
.
dot
(
g1Im
,
self
.
v
)
*
self
.
ds
(
0
)
+
fen
.
dot
(
g2Im
,
self
.
v
)
*
self
.
ds
(
1
))
b0Re
=
fen
.
assemble
(
L0Re
)
b0Im
=
fen
.
assemble
(
L0Im
)
fen
.
DirichletBC
(
self
.
V
,
self
.
DirichletDatum
[
0
],
self
.
DirichletBoundary
)
.
apply
(
b0Re
)
fen
.
DirichletBC
(
self
.
V
,
self
.
DirichletDatum
[
1
],
self
.
DirichletBoundary
)
.
apply
(
b0Im
)
self
.
b0
=
np
.
array
(
b0Re
[:]
+
1.j
*
b0Im
[:],
dtype
=
np
.
complex
)
self
.
bs
=
[
self
.
b0
]
Event Timeline
Log In to Comment