Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79441175
diapason_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, 10:21
Size
6 KB
Mime Type
text/x-python
Expires
Wed, Aug 28, 10:21 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20186862
Attached To
R6746 RationalROMPy
diapason_engine.py
View Options
import
numpy
as
np
import
fenics
as
fen
import
ufl
from
rrompy.hfengines.vector_linear_problem
import
\
LinearElasticityHelmholtzProblemEngine
as
LEHPE
from
rrompy.hfengines.vector_linear_problem
import
\
LinearElasticityHelmholtzProblemEngineDamped
as
LEHPED
class
DiapasonEngine
(
LEHPE
):
def
__init__
(
self
,
kappa
:
float
,
c
:
float
,
rho
:
float
,
E
:
float
,
nu
:
float
,
T
:
float
,
theta
:
float
,
phi
:
float
,
meshNo
:
int
=
1
,
deg
:
int
=
1
,
degree_threshold
:
int
=
np
.
inf
,
homogeneized
:
bool
=
False
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
super
()
.
__init__
(
mu0
=
kappa
,
degree_threshold
=
degree_threshold
,
homogeneized
=
homogeneized
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
mesh
=
fen
.
Mesh
(
"../data/mesh/diapason_{}.xml"
.
format
(
meshNo
))
subdomains
=
fen
.
MeshFunction
(
"size_t"
,
mesh
,
(
"../data/mesh/diapason_{}_physical_"
"region.xml"
)
.
format
(
meshNo
))
meshBall
=
fen
.
SubMesh
(
mesh
,
subdomains
,
2
)
meshFork
=
fen
.
SubMesh
(
mesh
,
subdomains
,
1
)
Hball
=
np
.
max
(
meshBall
.
coordinates
()[:,
1
])
#.00257
Ltot
=
np
.
max
(
mesh
.
coordinates
()[:,
1
])
#.1022
Lhandle
=
np
.
max
(
meshFork
.
coordinates
()[:,
1
])
#.026
Lrod
=
Ltot
-
Lhandle
#.0762
L2var
=
(
Lrod
/
4.
)
**
2.
Ehandle_ratio
=
3.
rhohandle_ratio
=
1.5
kWave
=
(
np
.
cos
(
theta
)
*
np
.
cos
(
phi
),
np
.
sin
(
phi
),
np
.
sin
(
theta
)
*
np
.
cos
(
phi
))
x
,
y
,
z
=
fen
.
SpatialCoordinate
(
mesh
)[:]
yCorr
=
y
-
Ltot
compPlane
=
kWave
[
0
]
*
x
+
kWave
[
1
]
*
y
+
kWave
[
2
]
*
z
xPlane
,
yPlane
,
zPlane
=
(
xx
-
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
xOrtho
,
yOrtho
,
zOrtho
=
(
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
forcingBase
=
(
T
/
(
2.
*
np
.
pi
*
L2var
)
**.
5
*
fen
.
exp
(
-
(
xPlane
**
2.
+
yPlane
**
2.
+
zPlane
**
2.
)
/
(
2.
*
L2var
)))
forcingWeight
=
np
.
real
(
kappa
)
/
c
*
(
xOrtho
+
yOrtho
+
zOrtho
)
neumannDatum
=
[
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
cos
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
)),
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
sin
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
))]
lambda_
=
E
*
nu
/
(
1.
+
nu
)
/
(
1.
-
2.
*
nu
)
mu_
=
E
/
(
1.
+
nu
)
lambda_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
lambda_
),
fen
.
Constant
(
Ehandle_ratio
*
lambda_
))
mu_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
mu_
),
fen
.
Constant
(
Ehandle_ratio
*
mu_
))
rho_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
rho
),
fen
.
Constant
(
rhohandle_ratio
*
rho
))
self
.
lambda_
=
lambda_eff
self
.
mu_
=
mu_eff
self
.
rho_
=
rho_eff
self
.
V
=
fen
.
VectorFunctionSpace
(
mesh
,
"P"
,
deg
)
self
.
DirichletBoundary
=
lambda
x
,
on_b
:
on_b
and
x
[
1
]
<
Hball
self
.
NeumannBoundary
=
"REST"
self
.
forcingTerm
=
neumannDatum
class
DiapasonEngineDamped
(
LEHPED
):
def
__init__
(
self
,
kappa
:
float
,
c
:
float
,
rho
:
float
,
E
:
float
,
nu
:
float
,
T
:
float
,
theta
:
float
,
phi
:
float
,
dampingEta
:
float
,
meshNo
:
int
=
1
,
deg
:
int
=
1
,
degree_threshold
:
int
=
np
.
inf
,
homogeneized
:
bool
=
False
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
super
()
.
__init__
(
mu0
=
kappa
,
degree_threshold
=
degree_threshold
,
homogeneized
=
homogeneized
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
self
.
eta
=
dampingEta
mesh
=
fen
.
Mesh
(
"../data/mesh/diapason_{}.xml"
.
format
(
meshNo
))
subdomains
=
fen
.
MeshFunction
(
"size_t"
,
mesh
,
(
"../data/mesh/diapason_{}_physical_"
"region.xml"
)
.
format
(
meshNo
))
meshBall
=
fen
.
SubMesh
(
mesh
,
subdomains
,
2
)
meshFork
=
fen
.
SubMesh
(
mesh
,
subdomains
,
1
)
Hball
=
np
.
max
(
meshBall
.
coordinates
()[:,
1
])
#.00257
Ltot
=
np
.
max
(
mesh
.
coordinates
()[:,
1
])
#.1022
Lhandle
=
np
.
max
(
meshFork
.
coordinates
()[:,
1
])
#.026
Lrod
=
Ltot
-
Lhandle
#.0762
L2var
=
(
Lrod
/
4.
)
**
2.
Ehandle_ratio
=
3.
rhohandle_ratio
=
1.5
kWave
=
(
np
.
cos
(
theta
)
*
np
.
cos
(
phi
),
np
.
sin
(
phi
),
np
.
sin
(
theta
)
*
np
.
cos
(
phi
))
x
,
y
,
z
=
fen
.
SpatialCoordinate
(
mesh
)[:]
yCorr
=
y
-
Ltot
compPlane
=
kWave
[
0
]
*
x
+
kWave
[
1
]
*
y
+
kWave
[
2
]
*
z
xPlane
,
yPlane
,
zPlane
=
(
xx
-
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
xOrtho
,
yOrtho
,
zOrtho
=
(
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
forcingBase
=
(
T
/
(
2.
*
np
.
pi
*
L2var
)
**.
5
*
fen
.
exp
(
-
(
xPlane
**
2.
+
yPlane
**
2.
+
zPlane
**
2.
)
/
(
2.
*
L2var
)))
forcingWeight
=
np
.
real
(
kappa
)
/
c
*
(
xOrtho
+
yOrtho
+
zOrtho
)
neumannDatum
=
[
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
cos
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
)),
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
sin
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
))]
lambda_
=
E
*
nu
/
(
1.
+
nu
)
/
(
1.
-
2.
*
nu
)
mu_
=
E
/
(
1.
+
nu
)
lambda_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
lambda_
),
fen
.
Constant
(
Ehandle_ratio
*
lambda_
))
mu_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
mu_
),
fen
.
Constant
(
Ehandle_ratio
*
mu_
))
rho_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
rho
),
fen
.
Constant
(
rhohandle_ratio
*
rho
))
self
.
lambda_
=
lambda_eff
self
.
mu_
=
mu_eff
self
.
rho_
=
rho_eff
self
.
V
=
fen
.
VectorFunctionSpace
(
mesh
,
"P"
,
deg
)
self
.
DirichletBoundary
=
lambda
x
,
on_b
:
on_b
and
x
[
1
]
<
Hball
self
.
NeumannBoundary
=
"REST"
self
.
forcingTerm
=
neumannDatum
Event Timeline
Log In to Comment