Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61165105
linear_elasticity_helmholtz_problem_engine_augmented.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, May 4, 22:57
Size
12 KB
Mime Type
text/x-python
Expires
Mon, May 6, 22:57 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17474669
Attached To
R6746 RationalROMPy
linear_elasticity_helmholtz_problem_engine_augmented.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
from
scipy.sparse
import
eye
,
bmat
,
block_diag
from
collections.abc
import
Iterable
from
.linear_elasticity_helmholtz_problem_engine
import
(
LinearElasticityHelmholtzProblemEngine
,
LinearElasticityHelmholtzProblemEngineDamped
)
from
rrompy.solver.fenics
import
(
augmentedElasticNormMatrix
,
augmentedElasticDualNormMatrix
)
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.exception_manager
import
RROMPyException
from
rrompy.parameter
import
parameterMap
as
pMap
__all__
=
[
'LinearElasticityHelmholtzProblemEngineAugmented'
,
'LinearElasticityHelmholtzProblemEngineDampedAugmented'
]
class
LinearElasticityHelmholtzProblemEngineAugmented
(
LinearElasticityHelmholtzProblemEngine
):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * mu * v = f in \Omega
mu * u = v in \overline{\Omega}
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
dsToBeSet: Whether ds needs to be set.
"""
def
__init__
(
self
,
*
args
,
**
kwargs
):
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
parameterMap
=
pMap
(
1.
,
self
.
npar
)
@property
def
spacedim
(
self
):
if
(
hasattr
(
self
,
"bs"
)
and
isinstance
(
self
.
bs
,
Iterable
)
and
self
.
bs
[
0
]
is
not
None
):
return
len
(
self
.
bs
[
0
])
return
2
*
super
()
.
spacedim
def
buildEnergyNormForm
(
self
):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng
(
self
,
"INIT"
,
"Assembling energy matrix."
,
20
)
self
.
energyNormMatrix
=
augmentedElasticNormMatrix
(
self
.
V
,
self
.
lambda_
[
0
],
self
.
mu_
[
0
])
vbMng
(
self
,
"DEL"
,
"Done assembling energy matrix."
,
20
)
def
buildEnergyNormDualForm
(
self
):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng
(
self
,
"INIT"
,
"Assembling energy dual matrix."
,
20
)
self
.
energyNormDualMatrix
=
augmentedElasticDualNormMatrix
(
self
.
V
,
self
.
lambda_
[
0
],
self
.
mu_
[
0
],
compressRank
=
self
.
_energyDualNormCompress
)
vbMng
(
self
,
"DEL"
,
"Done assembling energy dual matrix."
,
20
)
def
buildA
(
self
):
"""Build terms of operator of linear system."""
ANone
=
any
([
A
is
None
for
A
in
self
.
As
])
if
not
ANone
:
return
self
.
nAs
=
2
super
()
.
buildA
()
I
=
eye
(
self
.
spacedim
//
2
)
self
.
As
[
0
]
=
block_diag
((
self
.
As
[
0
],
I
),
format
=
"csr"
)
self
.
As
[
1
]
=
bmat
([[
None
,
self
.
As
[
1
]],
[
-
I
,
None
]],
format
=
"csr"
)
def
buildb
(
self
):
"""Build terms of operator of linear system."""
bNone
=
any
([
b
is
None
for
b
in
self
.
bs
])
if
not
bNone
:
return
self
.
nbs
=
1
dim
=
self
.
spacedim
//
2
super
()
.
buildb
()
self
.
bs
[
0
]
=
np
.
pad
(
self
.
bs
[
0
],
(
0
,
dim
),
"constant"
)
def
plot
(
self
,
u
,
warping
=
None
,
is_state
=
False
,
name
=
"u"
,
save
=
None
,
what
=
'all'
,
forceNewFile
=
True
,
saveFormat
=
"eps"
,
saveDPI
=
100
,
show
=
True
,
colorMap
=
"jet"
,
fenplotArgs
=
{},
**
figspecs
):
uh
=
u
[:
self
.
spacedim
//
2
]
if
is_state
or
self
.
isCEye
else
u
return
super
()
.
plot
(
uh
,
warping
,
is_state
,
name
,
save
,
what
,
forceNewFile
,
saveFormat
,
saveDPI
,
show
,
colorMap
,
fenplotArgs
,
**
figspecs
)
def
outParaview
(
self
,
u
,
warping
=
None
,
is_state
=
False
,
name
=
"u"
,
filename
=
"out"
,
time
=
0.
,
what
=
'all'
,
forceNewFile
=
True
,
folder
=
False
,
filePW
=
None
):
if
not
is_state
and
not
self
.
isCEye
:
raise
RROMPyException
((
"Cannot output to Paraview non-state "
"object."
))
return
super
()
.
outParaview
(
u
[:
self
.
spacedim
//
2
],
warping
,
is_state
,
name
,
filename
,
time
,
what
,
forceNewFile
,
folder
,
filePW
)
def
outParaviewTimeDomain
(
self
,
u
,
omega
,
warping
=
None
,
is_state
=
False
,
timeFinal
=
None
,
periodResolution
=
20
,
name
=
"u"
,
filename
=
"out"
,
forceNewFile
=
True
,
folder
=
False
):
if
not
is_state
and
not
self
.
isCEye
:
raise
RROMPyException
((
"Cannot output to Paraview non-state "
"object."
))
return
super
()
.
outParaviewTimeDomain
(
u
[:
self
.
spacedim
//
2
],
omega
,
warping
,
is_state
,
timeFinal
,
periodResolution
,
name
,
filename
,
forceNewFile
,
folder
)
class
LinearElasticityHelmholtzProblemEngineDampedAugmented
(
LinearElasticityHelmholtzProblemEngineDamped
):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * (mu - i * eta) * v = f in \Omega
mu * u = v in \overline{\Omega}
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
eta: Value of eta.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
dsToBeSet: Whether ds needs to be set.
"""
def
__init__
(
self
,
*
args
,
**
kwargs
):
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
nAs
=
2
self
.
_weight0
=
1.
@property
def
spacedim
(
self
):
if
(
hasattr
(
self
,
"bs"
)
and
isinstance
(
self
.
bs
,
Iterable
)
and
self
.
bs
[
0
]
is
not
None
):
return
len
(
self
.
bs
[
0
])
return
2
*
super
()
.
spacedim
def
buildEnergyNormForm
(
self
):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng
(
self
,
"INIT"
,
"Assembling energy matrix."
,
20
)
self
.
energyNormMatrix
=
augmentedElasticNormMatrix
(
self
.
V
,
self
.
lambda_
[
0
],
self
.
mu_
[
0
])
vbMng
(
self
,
"DEL"
,
"Done assembling energy matrix."
,
20
)
def
buildEnergyNormDualForm
(
self
):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng
(
self
,
"INIT"
,
"Assembling energy dual matrix."
,
20
)
self
.
energyNormDualMatrix
=
augmentedElasticDualNormMatrix
(
self
.
V
,
self
.
lambda_
[
0
],
self
.
mu_
[
0
],
compressRank
=
self
.
_energyDualNormCompress
)
vbMng
(
self
,
"DEL"
,
"Done assembling energy dual matrix."
,
20
)
def
buildA
(
self
):
"""Build terms of operator of linear system."""
ANone
=
any
([
A
is
None
for
A
in
self
.
As
])
if
not
ANone
:
return
self
.
nAs
=
3
super
()
.
buildA
()
self
.
_nAs
=
2
I
=
eye
(
self
.
spacedim
//
2
)
self
.
As
[
0
]
=
bmat
([[
self
.
As
[
0
],
self
.
_weight0
*
self
.
As
[
1
]],
[
None
,
I
]],
format
=
"csr"
)
self
.
As
[
1
]
=
bmat
([[(
1.
-
self
.
_weight0
)
*
self
.
As
[
1
],
self
.
As
[
2
]],
[
-
I
,
None
]],
format
=
"csr"
)
self
.
thAs
.
pop
()
self
.
As
.
pop
()
def
buildb
(
self
):
"""Build terms of operator of linear system."""
bNone
=
any
([
b
is
None
for
b
in
self
.
bs
])
if
not
bNone
:
return
self
.
nbs
=
1
dim
=
self
.
spacedim
//
2
super
()
.
buildb
()
self
.
bs
[
0
]
=
np
.
pad
(
self
.
bs
[
0
],
(
0
,
dim
),
"constant"
)
def
plot
(
self
,
u
,
warping
=
None
,
is_state
=
False
,
name
=
"u"
,
save
=
None
,
what
=
'all'
,
forceNewFile
=
True
,
saveFormat
=
"eps"
,
saveDPI
=
100
,
show
=
True
,
colorMap
=
"jet"
,
fenplotArgs
=
{},
**
figspecs
):
uh
=
u
[:
self
.
spacedim
//
2
]
if
is_state
or
self
.
isCEye
else
u
return
super
()
.
plot
(
uh
,
warping
,
is_state
,
name
,
save
,
what
,
forceNewFile
,
saveFormat
,
saveDPI
,
show
,
colorMap
,
fenplotArgs
,
**
figspecs
)
def
outParaview
(
self
,
u
,
warping
=
None
,
is_state
=
False
,
name
=
"u"
,
filename
=
"out"
,
time
=
0.
,
what
=
'all'
,
forceNewFile
=
True
,
folder
=
False
,
filePW
=
None
):
if
not
is_state
and
not
self
.
isCEye
:
raise
RROMPyException
((
"Cannot output to Paraview non-state "
"object."
))
return
super
()
.
outParaview
(
u
[:
self
.
spacedim
//
2
],
warping
,
is_state
,
name
,
filename
,
time
,
what
,
forceNewFile
,
folder
,
filePW
)
def
outParaviewTimeDomain
(
self
,
u
,
omega
,
warping
=
None
,
is_state
=
False
,
timeFinal
=
None
,
periodResolution
=
20
,
name
=
"u"
,
filename
=
"out"
,
forceNewFile
=
True
,
folder
=
False
):
if
not
is_state
and
not
self
.
isCEye
:
raise
RROMPyException
((
"Cannot output to Paraview non-state "
"object."
))
return
super
()
.
outParaviewTimeDomain
(
u
[:
self
.
spacedim
//
2
],
omega
,
warping
,
is_state
,
timeFinal
,
periodResolution
,
name
,
filename
,
forceNewFile
,
folder
)
Event Timeline
Log In to Comment