Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61177607
linear_elasticity_helmholtz_problem_engine_damped.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
Sun, May 5, 00:55
Size
10 KB
Mime Type
text/x-python
Expires
Tue, May 7, 00:55 (2 d)
Engine
blob
Format
Raw Data
Handle
17476339
Attached To
R6746 RationalROMPy
linear_elasticity_helmholtz_problem_engine_damped.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
csr_matrix
import
fenics
as
fen
from
.linear_elasticity_helmholtz_problem_engine
import
\
LinearElasticityHelmholtzProblemEngine
from
rrompy.utilities.base.types
import
ScOp
,
paramVal
from
rrompy.solver.fenics
import
fenZERO
,
fenZEROS
from
rrompy.utilities.base
import
verbosityDepth
__all__
=
[
'LinearElasticityHelmholtzProblemEngineDamped'
]
class
LinearElasticityHelmholtzProblemEngineDamped
(
LinearElasticityHelmholtzProblemEngine
):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * (mu^2 - i * eta * mu) * u = f in \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.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
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).
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.
dsToBeSet: Whether ds needs to be set.
"""
nAs
=
3
rescalingExp
=
1.
def
__init__
(
self
,
degree_threshold
:
int
=
np
.
inf
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
super
()
.
__init__
(
degree_threshold
=
degree_threshold
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
self
.
eta
=
fenZERO
@property
def
eta
(
self
):
"""Value of eta."""
return
self
.
_eta
@eta.setter
def
eta
(
self
,
eta
):
self
.
resetAs
()
if
not
isinstance
(
eta
,
(
list
,
tuple
,)):
eta
=
[
eta
,
fenZERO
]
self
.
_eta
=
eta
def
A
(
self
,
mu
:
paramVal
=
(),
der
:
int
=
0
)
->
ScOp
:
"""Assemble (derivative of) operator of linear system."""
mu
=
self
.
checkParameter
(
mu
)
Anull
=
self
.
checkAInBounds
(
der
)
if
Anull
is
not
None
:
return
Anull
self
.
autoSetDS
()
if
der
<=
0
and
self
.
As
[
0
]
is
None
:
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"INIT"
,
"Assembling operator term A0."
,
timestamp
=
self
.
timestamp
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
self
.
V
.
mesh
()
.
topology
()
.
dim
()),
self
.
DirichletBoundary
)
lambda_Re
,
lambda_Im
=
self
.
lambda_
mu_Re
,
mu_Im
=
self
.
mu_
hRe
,
hIm
=
self
.
RobinDatumH
termNames
=
[
"lambda_"
,
"mu_"
,
"RobinDatumH"
]
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
(
[
lambda_Re
,
mu_Re
,
hRe
],
[
x
+
"Real"
for
x
in
termNames
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
(
[
lambda_Im
,
mu_Re
,
hIm
],
[
x
+
"Imag"
for
x
in
termNames
]))
epsilon
=
lambda
u
:
0.5
*
(
fen
.
grad
(
u
)
+
fen
.
nabla_grad
(
u
))
sigma
=
lambda
u
,
l_
,
m_
:
(
l_
*
fen
.
div
(
u
)
*
fen
.
Identity
(
u
.
geometric_dimension
())
+
2.
*
m_
*
epsilon
(
u
))
a0Re
=
(
fen
.
inner
(
sigma
(
self
.
u
,
lambda_Re
,
mu_Re
),
epsilon
(
self
.
v
))
*
fen
.
dx
+
hRe
*
fen
.
inner
(
self
.
u
,
self
.
v
)
*
self
.
ds
(
1
))
a0Im
=
(
fen
.
inner
(
sigma
(
self
.
u
,
lambda_Im
,
mu_Im
),
epsilon
(
self
.
v
))
*
fen
.
dx
+
hIm
*
fen
.
inner
(
self
.
u
,
self
.
v
)
*
self
.
ds
(
1
))
A0Re
=
fen
.
assemble
(
a0Re
,
form_compiler_parameters
=
parsRe
)
A0Im
=
fen
.
assemble
(
a0Im
,
form_compiler_parameters
=
parsIm
)
DirichletBC0
.
apply
(
A0Re
)
DirichletBC0
.
zero
(
A0Im
)
A0ReMat
=
fen
.
as_backend_type
(
A0Re
)
.
mat
()
A0ImMat
=
fen
.
as_backend_type
(
A0Im
)
.
mat
()
A0Rer
,
A0Rec
,
A0Rev
=
A0ReMat
.
getValuesCSR
()
A0Imr
,
A0Imc
,
A0Imv
=
A0ImMat
.
getValuesCSR
()
self
.
As
[
0
]
=
(
csr_matrix
((
A0Rev
,
A0Rec
,
A0Rer
),
shape
=
A0ReMat
.
size
)
+
1.j
*
csr_matrix
((
A0Imv
,
A0Imc
,
A0Imr
),
shape
=
A0ImMat
.
size
))
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"DEL"
,
"Done assembling operator term."
,
timestamp
=
self
.
timestamp
)
if
der
<=
1
and
self
.
As
[
1
]
is
None
:
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"INIT"
,
"Assembling operator term A1."
,
timestamp
=
self
.
timestamp
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
self
.
V
.
mesh
()
.
topology
()
.
dim
()),
self
.
DirichletBoundary
)
rho_Re
,
rho_Im
=
self
.
rho_
eta_Re
,
eta_Im
=
self
.
eta
termNames
=
[
"rho_"
,
"eta"
]
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
rho_Re
,
eta_Re
],
[
x
+
"Real"
for
x
in
termNames
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
rho_Im
,
eta_Im
],
[
x
+
"Imag"
for
x
in
termNames
]))
a1Re
=
-
((
eta_Re
*
rho_Im
+
eta_Im
*
rho_Re
)
*
fen
.
inner
(
self
.
u
,
self
.
v
))
*
fen
.
dx
a1Im
=
((
eta_Re
*
rho_Re
-
eta_Im
*
rho_Im
)
*
fen
.
inner
(
self
.
u
,
self
.
v
))
*
fen
.
dx
A1Re
=
fen
.
assemble
(
a1Re
,
form_compiler_parameters
=
parsRe
)
A1Im
=
fen
.
assemble
(
a1Im
,
form_compiler_parameters
=
parsIm
)
DirichletBC0
.
zero
(
A1Re
)
DirichletBC0
.
zero
(
A1Im
)
A1ReMat
=
fen
.
as_backend_type
(
A1Re
)
.
mat
()
A1ImMat
=
fen
.
as_backend_type
(
A1Im
)
.
mat
()
A1Rer
,
A1Rec
,
A1Rev
=
A1ReMat
.
getValuesCSR
()
A1Imr
,
A1Imc
,
A1Imv
=
A1ImMat
.
getValuesCSR
()
self
.
As
[
1
]
=
(
csr_matrix
((
A1Rev
,
A1Rec
,
A1Rer
),
shape
=
A1ReMat
.
size
)
+
1.j
*
csr_matrix
((
A1Imv
,
A1Imc
,
A1Imr
),
shape
=
A1ImMat
.
size
))
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"DEL"
,
"Done assembling operator term."
,
timestamp
=
self
.
timestamp
)
if
der
<=
2
and
self
.
As
[
2
]
is
None
:
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"INIT"
,
"Assembling operator term A2."
,
timestamp
=
self
.
timestamp
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
self
.
V
.
mesh
()
.
topology
()
.
dim
()),
self
.
DirichletBoundary
)
rho_Re
,
rho_Im
=
self
.
rho_
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
rho_Re
],
[
"rho_Real"
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
rho_Im
],
[
"rho_Imag"
]))
a2Re
=
-
rho_Re
*
fen
.
inner
(
self
.
u
,
self
.
v
)
*
fen
.
dx
a2Im
=
-
rho_Im
*
fen
.
inner
(
self
.
u
,
self
.
v
)
*
fen
.
dx
A2Re
=
fen
.
assemble
(
a2Re
,
form_compiler_parameters
=
parsRe
)
A2Im
=
fen
.
assemble
(
a2Im
,
form_compiler_parameters
=
parsIm
)
DirichletBC0
.
zero
(
A2Re
)
DirichletBC0
.
zero
(
A2Im
)
A2ReMat
=
fen
.
as_backend_type
(
A2Re
)
.
mat
()
A2ImMat
=
fen
.
as_backend_type
(
A2Im
)
.
mat
()
A2Rer
,
A2Rec
,
A2Rev
=
A2ReMat
.
getValuesCSR
()
A2Imr
,
A2Imc
,
A2Imv
=
A2ImMat
.
getValuesCSR
()
self
.
As
[
2
]
=
(
csr_matrix
((
A2Rev
,
A2Rec
,
A2Rer
),
shape
=
A2ReMat
.
size
)
+
1.j
*
csr_matrix
((
A2Imv
,
A2Imc
,
A2Imr
),
shape
=
A2ImMat
.
size
))
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"DEL"
,
"Done assembling operator term."
,
timestamp
=
self
.
timestamp
)
if
der
==
0
:
return
self
.
As
[
0
]
+
mu
(
0
)
*
self
.
As
[
1
]
+
mu
(
0
)
**
2.
*
self
.
As
[
2
]
if
der
==
1
:
return
self
.
As
[
1
]
+
2
*
mu
(
0
)
*
self
.
As
[
2
]
return
self
.
As
[
2
]
Event Timeline
Log In to Comment