Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81249947
linear_elasticity_beam_poisson_ratio_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
Wed, Sep 4, 23:33
Size
5 KB
Mime Type
text/x-python
Expires
Fri, Sep 6, 23:33 (2 d)
Engine
blob
Format
Raw Data
Handle
20536910
Attached To
R6746 RationalROMPy
linear_elasticity_beam_poisson_ratio_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
fenics
as
fen
from
rrompy.hfengines.vector_linear_problem
import
\
LinearElasticityProblemEngine
from
rrompy.solver.fenics
import
fenZERO
,
fenZEROS
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.solver.fenics
import
fenics2Sparse
,
fenics2Vector
class
LinearElasticityBeamPoissonRatioEngine
(
LinearElasticityProblemEngine
):
"""
Solver for linear elasticity problem of a beam subject to its own weight,
with parametric Poisson's ratio.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = rho_ * g in \Omega
u = 0 on \Gamma_D
\partial_nu = 0 on \Gamma_N
"""
def
__init__
(
self
,
n
:
int
,
rho_
:
float
,
g
:
float
,
E
:
float
,
nu0
:
float
,
length
:
float
,
*
args
,
**
kwargs
):
super
()
.
__init__
(
nu0
,
*
args
,
**
kwargs
)
self
.
_affinePoly
=
False
self
.
nAs
,
self
.
nbs
=
3
,
2
self
.
E_
=
E
mesh
=
fen
.
RectangleMesh
(
fen
.
Point
(
0.
,
0.
),
fen
.
Point
(
length
,
1
),
n
,
max
(
int
(
n
/
length
),
1
))
self
.
V
=
fen
.
VectorFunctionSpace
(
mesh
,
"P"
,
1
)
self
.
forcingTerm
=
[
fen
.
Constant
((
0.
,
-
rho_
*
g
/
E
)),
fenZEROS
(
2
)]
self
.
DirichletBoundary
=
lambda
x
,
on_b
:
on_b
and
fen
.
near
(
x
[
0
],
0.
)
self
.
NeumannBoundary
=
"REST"
def
buildA
(
self
):
"""Build terms of operator of linear system."""
if
self
.
As
[
0
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A0."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
2
),
self
.
DirichletBoundary
)
a0Re
=
fenZERO
*
fen
.
inner
(
self
.
u
,
self
.
v
)
*
fen
.
dx
self
.
As
[
0
]
=
fenics2Sparse
(
a0Re
,
{},
DirichletBC0
,
1
)
self
.
thAs
[
0
]
=
self
.
getMonomialSingleWeight
([
0
])
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
if
self
.
As
[
1
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A1."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
2
),
self
.
DirichletBoundary
)
epsilon
=
lambda
u
:
.
5
*
(
fen
.
grad
(
u
)
+
fen
.
nabla_grad
(
u
))
a1Re
=
self
.
E_
*
fen
.
inner
(
epsilon
(
self
.
u
),
epsilon
(
self
.
v
))
*
fen
.
dx
self
.
As
[
1
]
=
fenics2Sparse
(
a1Re
,
{},
DirichletBC0
,
0
)
self
.
thAs
[
1
]
=
[(
'x'
,
'()'
,
0
,
'*'
,
-
2.
,
'+'
,
1.
),
(
-
2.0
,),
None
]
epsilon
=
lambda
u
:
.
5
*
(
fen
.
grad
(
u
)
+
fen
.
nabla_grad
(
u
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
if
self
.
As
[
2
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A2."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
2
),
self
.
DirichletBoundary
)
a2Re
=
self
.
E_
*
fen
.
div
(
self
.
u
)
*
fen
.
div
(
self
.
v
)
*
fen
.
dx
self
.
As
[
2
]
=
fenics2Sparse
(
a2Re
,
{},
DirichletBC0
,
0
)
self
.
thAs
[
2
]
=
self
.
getMonomialSingleWeight
([
1
])
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
def
buildb
(
self
):
"""Build terms of operator of linear system."""
if
self
.
thbs
[
0
]
is
None
:
self
.
thbs
=
[
None
]
*
self
.
nbs
if
self
.
bs
[
0
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling forcing term b0."
,
20
)
L0Re
=
fen
.
inner
(
fenZEROS
(
2
),
self
.
v
)
*
fen
.
dx
DBCR
=
fen
.
DirichletBC
(
self
.
V
,
self
.
DirichletDatum
[
0
],
self
.
DirichletBoundary
)
DBCI
=
fen
.
DirichletBC
(
self
.
V
,
self
.
DirichletDatum
[
1
],
self
.
DirichletBoundary
)
self
.
bs
[
0
]
=
(
fenics2Vector
(
L0Re
,
{},
DBCR
,
1
)
+
1.j
*
fenics2Vector
(
L0Re
,
{},
DBCI
,
1
))
self
.
thbs
[
0
]
=
self
.
getMonomialSingleWeight
([
0
])
if
self
.
bs
[
1
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling forcing term b1."
,
20
)
fRe
,
fIm
=
self
.
forcingTerm
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
fRe
],
[
"forcingTermReal"
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
fIm
],
[
"forcingTermImag"
]))
L1Re
=
fen
.
inner
(
fRe
,
self
.
v
)
*
fen
.
dx
L1Im
=
fen
.
inner
(
fIm
,
self
.
v
)
*
fen
.
dx
DBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZEROS
(
2
),
self
.
DirichletBoundary
)
self
.
bs
[
1
]
=
(
fenics2Vector
(
L1Re
,
parsRe
,
DBC0
,
1
)
+
1.j
*
fenics2Vector
(
L1Im
,
parsIm
,
DBC0
,
1
))
self
.
thbs
[
1
]
=
[(
'x'
,
'()'
,
0
,
'**'
,
2.
,
'*'
,
-
2.
,
'-'
,
(
'x'
,
'()'
,
0
),
'+'
,
1.
),
(
'x'
,
'()'
,
0
,
'*'
,
-
4.
,
'-'
,
1.
),
(
-
2.0
,),
None
]
vbMng
(
self
,
"DEL"
,
"Done assembling forcing term."
,
20
)
Event Timeline
Log In to Comment