Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61289596
fracture_pod.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, 17:48
Size
7 KB
Mime Type
text/x-python
Expires
Tue, May 7, 17:48 (2 d)
Engine
blob
Format
Raw Data
Handle
17491063
Attached To
R6746 RationalROMPy
fracture_pod.py
View Options
import
numpy
as
np
from
membrane_fracture_engine
import
MembraneFractureEngine
as
MFE
from
rrompy.reduction_methods.greedy
import
RationalInterpolantGreedy
as
RIG
from
rrompy.reduction_methods.greedy
import
ReducedBasisGreedy
as
RBG
from
rrompy.parameter.parameter_sampling
import
(
QuadratureSampler
as
QS
,
QuadratureSamplerTotal
as
QST
)
verb
=
15
size
=
16
show_sample
=
False
show_norm
=
True
MN
=
1
S
=
(
MN
+
2
)
*
(
MN
+
1
)
//
2
greedyTol
=
1e-1
collinearityTol
=
0.
nTestPoints
=
900
algo
=
"rational"
#algo = "RB"
if
algo
==
"rational"
:
radial
=
""
#radial = "_gaussian"
#radial = "_thinplate"
#radial = "_multiquadric"
rW0
=
5.
radialWeight
=
[
rW0
]
*
2
polybasis
=
"CHEBYSHEV"
polybasis
=
"LEGENDRE"
#polybasis = "MONOMIAL"
errorEstimatorKind
=
'AFFINE'
errorEstimatorKind
=
'DISCREPANCY'
errorEstimatorKind
=
'INTERPOLATORY'
# errorEstimatorKind = 'EIM_DIAGONAL'
errorEstimatorKind
=
'EIM_INTERPOLATORY'
if
size
==
1
:
# below
mu0
=
[
40
**
.
5
,
.
4
]
mutar
=
[
45
**
.
5
,
.
4
]
murange
=
[[
30
**
.
5
,
.
3
],
[
50
**
.
5
,
.
5
]]
elif
size
==
2
:
# top
mu0
=
[
40
**
.
5
,
.
6
]
mutar
=
[
45
**
.
5
,
.
6
]
murange
=
[[
30
**
.
5
,
.
5
],
[
50
**
.
5
,
.
7
]]
elif
size
==
3
:
# interesting
mu0
=
[
40
**
.
5
,
.
5
]
mutar
=
[
45
**
.
5
,
.
5
]
murange
=
[[
30
**
.
5
,
.
3
],
[
50
**
.
5
,
.
7
]]
elif
size
==
4
:
# wide_low
mu0
=
[
40
**
.
5
,
.
2
]
mutar
=
[
45
**
.
5
,
.
2
]
murange
=
[[
10
**
.
5
,
.
1
],
[
70
**
.
5
,
.
3
]]
elif
size
==
5
:
# wide_hi
mu0
=
[
40
**
.
5
,
.
8
]
mutar
=
[
45
**
.
5
,
.
8
]
murange
=
[[
10
**
.
5
,
.
7
],
[
70
**
.
5
,
.
9
]]
elif
size
==
6
:
# top_zoom
mu0
=
[
50
**
.
5
,
.
8
]
mutar
=
[
55
**
.
5
,
.
8
]
murange
=
[[
40
**
.
5
,
.
7
],
[
60
**
.
5
,
.
9
]]
elif
size
==
7
:
# huge
mu0
=
[
50
**
.
5
,
.
5
]
mutar
=
[
55
**
.
5
,
.
5
]
murange
=
[[
10
**
.
5
,
.
2
],
[
90
**
.
5
,
.
8
]]
elif
size
==
11
:
#L below
mu0
=
[
110
**
.
5
,
.
4
]
mutar
=
[
115
**
.
5
,
.
4
]
murange
=
[[
90
**
.
5
,
.
3
],
[
130
**
.
5
,
.
5
]]
elif
size
==
12
:
#L top
mu0
=
[
110
**
.
5
,
.
6
]
mutar
=
[
115
**
.
5
,
.
6
]
murange
=
[[
90
**
.
5
,
.
5
],
[
130
**
.
5
,
.
7
]]
elif
size
==
13
:
#L interesting
mu0
=
[
110
**
.
5
,
.
5
]
mutar
=
[
115
**
.
5
,
.
5
]
murange
=
[[
90
**
.
5
,
.
3
],
[
130
**
.
5
,
.
7
]]
elif
size
==
14
:
#L belowbelow
mu0
=
[
110
**
.
5
,
.
2
]
mutar
=
[
115
**
.
5
,
.
2
]
murange
=
[[
90
**
.
5
,
.
1
],
[
130
**
.
5
,
.
3
]]
elif
size
==
15
:
#L toptop
mu0
=
[
110
**
.
5
,
.
8
]
mutar
=
[
115
**
.
5
,
.
8
]
murange
=
[[
90
**
.
5
,
.
7
],
[
130
**
.
5
,
.
9
]]
elif
size
==
16
:
#L interestinginteresting
mu0
=
[
110
**
.
5
,
.
5
]
mutar
=
[
115
**
.
5
,
.
6
]
murange
=
[[
90
**
.
5
,
.
1
],
[
130
**
.
5
,
.
9
]]
elif
size
==
17
:
#L interestingtop
mu0
=
[
110
**
.
5
,
.
7
]
mutar
=
[
115
**
.
5
,
.
6
]
murange
=
[[
90
**
.
5
,
.
5
],
[
130
**
.
5
,
.
9
]]
elif
size
==
18
:
#L interestingbelow
mu0
=
[
110
**
.
5
,
.
3
]
mutar
=
[
115
**
.
5
,
.
4
]
murange
=
[[
90
**
.
5
,
.
1
],
[
130
**
.
5
,
.
5
]]
elif
size
==
100
:
# tiny
mu0
=
[
32.5
**
.
5
,
.
5
]
mutar
=
[
34
**
.
5
,
.
5
]
murange
=
[[
30
**
.
5
,
.
3
],
[
35
**
.
5
,
.
7
]]
aEff
=
1.05
bEff
=
1.
-
aEff
murangeEff
=
[[(
aEff
*
murange
[
0
][
0
]
**
2.
+
bEff
*
murange
[
1
][
0
]
**
2.
)
**
.
5
,
aEff
*
murange
[
0
][
1
]
+
bEff
*
murange
[
1
][
1
]],
[(
aEff
*
murange
[
1
][
0
]
**
2.
+
bEff
*
murange
[
0
][
0
]
**
2.
)
**
.
5
,
aEff
*
murange
[
1
][
1
]
+
bEff
*
murange
[
0
][
1
]]]
H
=
1.
L
=
.
75
delta
=
.
05
n
=
20
solver
=
MFE
(
mu0
=
mu0
,
H
=
H
,
L
=
L
,
delta
=
delta
,
n
=
n
,
verbosity
=
verb
)
rescalingExp
=
[
2.
,
1.
]
if
algo
==
"rational"
:
params
=
{
'S'
:
S
,
'POD'
:
True
,
'greedyTol'
:
greedyTol
,
'sampler'
:
QS
(
murange
,
'UNIFORM'
,
scalingExp
=
rescalingExp
),
'nTestPoints'
:
nTestPoints
,
'polybasis'
:
polybasis
+
radial
,
'radialDirectionalWeights'
:
radialWeight
,
'errorEstimatorKind'
:
errorEstimatorKind
,
'trainSetGenerator'
:
QST
(
murange
,
'CHEBYSHEV'
,
scalingExp
=
rescalingExp
)}
method
=
RIG
else
:
#if algo == "RB":
params
=
{
'S'
:
S
,
'POD'
:
True
,
'greedyTol'
:
greedyTol
,
'sampler'
:
QS
(
murange
,
'UNIFORM'
,
scalingExp
=
rescalingExp
),
'nTestPoints'
:
nTestPoints
,
'trainSetGenerator'
:
QST
(
murange
,
'CHEBYSHEV'
,
scalingExp
=
rescalingExp
)}
method
=
RBG
approx
=
method
(
solver
,
mu0
=
mu0
,
approxParameters
=
params
,
verbosity
=
verb
)
approx
.
greedy
(
True
)
if
show_sample
:
import
ufl
import
fenics
as
fen
from
rrompy.solver.fenics
import
affine_warping
L
=
mutar
[
1
]
y
=
fen
.
SpatialCoordinate
(
solver
.
V
.
mesh
())[
1
]
warp1
,
warpI1
=
affine_warping
(
solver
.
V
.
mesh
(),
np
.
array
([[
1
,
0
],
[
0
,
2.
*
L
]]))
warp2
,
warpI2
=
affine_warping
(
solver
.
V
.
mesh
(),
np
.
array
([[
1
,
0
],
[
0
,
2.
-
2.
*
L
]]))
warp
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
0.
),
warp1
,
warp2
)
warpI
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
0.
),
warpI1
,
warpI2
)
approx
.
plotApprox
(
mutar
,
[
warp
,
warpI
],
name
=
'u_app'
,
what
=
"REAL"
)
approx
.
plotHF
(
mutar
,
[
warp
,
warpI
],
name
=
'u_HF'
,
what
=
"REAL"
)
approx
.
plotErr
(
mutar
,
[
warp
,
warpI
],
name
=
'err'
,
what
=
"REAL"
)
# approx.plotRes(mutar, [warp, warpI], name = 'res', what = "REAL")
appErr
=
approx
.
normErr
(
mutar
)
solNorm
=
approx
.
normHF
(
mutar
)
resNorm
=
approx
.
normRes
(
mutar
)
RHSNorm
=
approx
.
normRHS
(
mutar
)
print
((
'SolNorm:
\t
{}
\n
Err:
\t
{}
\n
ErrRel:
\t
{}'
)
.
format
(
solNorm
,
appErr
,
np
.
divide
(
appErr
,
solNorm
)))
print
((
'RHSNorm:
\t
{}
\n
Res:
\t
{}
\n
ResRel:
\t
{}'
)
.
format
(
RHSNorm
,
resNorm
,
np
.
divide
(
resNorm
,
RHSNorm
)))
verb
=
approx
.
trainedModel
.
verbosity
approx
.
trainedModel
.
verbosity
=
5
from
plot_zero_set
import
plotZeroSet2
muZeroVals
,
Qvals
=
plotZeroSet2
(
murange
,
murangeEff
,
approx
,
mu0
,
200
,
[
2.
,
1.
])
if
show_norm
:
from
plot_inf_set
import
plotInfSet2
muInfVals
,
normEx
,
normApp
,
normRes
,
normErr
,
beta
=
plotInfSet2
(
murange
,
murangeEff
,
approx
,
mu0
,
50
,
[
2.
,
1.
],
relative
=
True
,
nobeta
=
True
)
try
:
QV
=
approx
.
trainedModel
.
getQVal
(
muInfVals
)
import
warnings
from
matplotlib
import
pyplot
as
plt
mu2x
,
mu2y
=
approx
.
mus
(
0
)
**
2.
,
approx
.
mus
(
1
)
**
1.
murangeExp
=
[[
murange
[
0
][
0
]
**
2.
,
murange
[
0
][
1
]],
[
murange
[
1
][
0
]
**
2.
,
murange
[
1
][
1
]]]
mu1s
=
np
.
unique
([
m
[
0
]
for
m
in
muInfVals
])
mu2s
=
np
.
unique
([
m
[
1
]
for
m
in
muInfVals
])
mu1
=
np
.
power
(
mu1s
,
2.
)
mu2
=
np
.
power
(
mu2s
,
1.
)
Mu1
,
Mu2
=
np
.
meshgrid
(
np
.
real
(
mu1
),
np
.
real
(
mu2
))
Reste
=
(
approx
.
errorEstimator
(
muInfVals
)
*
QV
)
.
reshape
(
normEx
.
shape
)
Rest
=
np
.
log10
(
Reste
)
Restmin
,
Restmax
=
np
.
min
(
Rest
),
np
.
max
(
Rest
)
cmap
=
plt
.
cm
.
jet
warnings
.
simplefilter
(
"ignore"
,
category
=
(
UserWarning
,
np
.
ComplexWarning
))
plt
.
figure
(
figsize
=
(
15
,
7
))
plt
.
jet
()
p
=
plt
.
contourf
(
Mu1
,
Mu2
,
Rest
,
levels
=
np
.
linspace
(
Restmin
,
Restmax
,
50
))
plt
.
plot
(
np
.
real
(
mu2x
),
np
.
real
(
mu2y
),
'kx'
)
for
j
,
xy
in
enumerate
(
zip
(
np
.
real
(
mu2x
),
np
.
real
(
mu2y
))):
plt
.
annotate
(
"{}"
.
format
(
j
),
xy
)
plt
.
plot
([
murangeExp
[
0
][
0
]]
*
2
,
[
murangeExp
[
0
][
1
],
murangeExp
[
1
][
1
]],
'm:'
)
plt
.
plot
([
murangeExp
[
0
][
0
],
murangeExp
[
1
][
0
]],
[
murangeExp
[
1
][
1
]]
*
2
,
'm:'
)
plt
.
plot
([
murangeExp
[
1
][
0
]]
*
2
,
[
murangeExp
[
1
][
1
],
murangeExp
[
0
][
1
]],
'm:'
)
plt
.
plot
([
murangeExp
[
1
][
0
],
murangeExp
[
0
][
0
]],
[
murangeExp
[
0
][
1
]]
*
2
,
'm:'
)
plt
.
title
(
"log10|res_est(mu)|"
)
plt
.
colorbar
(
p
)
plt
.
grid
()
plt
.
show
()
except
:
pass
approx
.
trainedModel
.
verbosity
=
verb
try
:
print
(
np
.
sort
(
approx
.
getPoles
([
None
,
.
5
])
**
2.
))
except
:
pass
Event Timeline
Log In to Comment