Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62976631
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
Thu, May 16, 21:48
Size
9 KB
Mime Type
text/x-python
Expires
Sat, May 18, 21:48 (2 d)
Engine
blob
Format
Raw Data
Handle
17715269
Attached To
R6746 RationalROMPy
fracture_pod.py
View Options
import
numpy
as
np
from
rrompy.hfengines.linear_problem.bidimensional
import
\
MembraneFractureEngine
as
MFE
from
rrompy.reduction_methods.standard
import
RationalInterpolant
as
RI
from
rrompy.reduction_methods.standard
import
ReducedBasis
as
RB
from
rrompy.reduction_methods.pivoted
import
RationalInterpolantPivoted
as
RIP
#from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from
rrompy.parameter.parameter_sampling
import
(
QuadratureSampler
as
QS
,
QuadratureSamplerTotal
as
QST
,
ManualSampler
as
MS
,
RandomSampler
as
RS
)
verb
=
10
size
=
2
show_sample
=
True
show_norm
=
True
Delta
=
0
MN
=
5
R
=
(
MN
+
2
)
*
(
MN
+
1
)
//
2
STensorized
=
(
MN
+
1
)
**
2
PODTol
=
1e-10
SPivot
=
MN
+
1
SMarginal
=
4
MMarginal
=
SMarginal
-
1
matchingWeight
=
1.
cutOffTolerance
=
1.
*
np
.
inf
cutOffType
=
"MAGNITUDE"
samples
=
"centered"
samples
=
"standard"
#samples = "pivoted"
algo
=
"rational"
algo
=
"RB"
sampling
=
"quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM
=
"quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype
=
"TOTAL"
polydegreetype
=
"FULL"
if
samples
in
[
"standard"
,
"pivoted"
]:
radial
=
""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0
=
5.
radialWeight
=
[
rW0
]
if
samples
==
"pivoted"
:
radialM
=
""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0
=
2.
radialWeightM
=
[
rMW0
]
assert
Delta
<=
0
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.
#25
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
=
{
'N'
:
MN
,
'M'
:
MN
+
Delta
,
'S'
:
R
,
'POD'
:
True
,
'polydegreetype'
:
polydegreetype
}
if
samples
==
"standard"
:
params
[
'polybasis'
]
=
"CHEBYSHEV"
+
radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params
[
'radialDirectionalWeights'
]
=
radialWeight
*
2
method
=
RI
elif
samples
==
"centered"
:
params
[
'polybasis'
]
=
"MONOMIAL"
method
=
RI
elif
samples
==
"pivoted"
:
params
[
'S'
]
=
SPivot
params
[
'SMarginal'
]
=
SMarginal
params
[
'MMarginal'
]
=
MMarginal
params
[
'polybasisPivot'
]
=
"CHEBYSHEV"
+
radial
params
[
'polybasisMarginal'
]
=
"MONOMIAL"
+
radialM
params
[
'radialDirectionalWeightsPivot'
]
=
radialWeight
params
[
'radialDirectionalWeightsMarginal'
]
=
radialWeightM
params
[
'matchingWeight'
]
=
matchingWeight
params
[
'cutOffTolerance'
]
=
cutOffTolerance
params
[
"cutOffType"
]
=
cutOffType
method
=
RIP
else
:
#if algo == "RB":
params
=
{
'R'
:(
MN
+
2
+
Delta
)
*
(
MN
+
1
+
Delta
)
//
2
,
'S'
:
R
,
'POD'
:
True
,
'PODTolerance'
:
PODTol
}
if
samples
==
"standard"
:
method
=
RB
elif
samples
==
"centered"
:
method
=
RB
elif
samples
==
"pivoted"
:
raise
# params['S'] = SPivot
# params['R'] = SPivot
# params['SMarginal'] = SMarginal
# params['MMarginal'] = MMarginal
# params['polybasisMarginal'] = "MONOMIAL" + radialM
# params['radialDirectionalWeightsMarginal'] = radialWeightM
# params['matchingWeight'] = matchingWeight
# params['cutOffTolerance'] = cutOffTolerance
# params["cutOffType"] = cutOffType
# method = RBP
if
samples
==
"standard"
:
if
sampling
==
"quadrature"
:
params
[
'sampler'
]
=
QS
(
murange
,
"CHEBYSHEV"
,
scalingExp
=
rescalingExp
)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params
[
'S'
]
=
STensorized
elif
sampling
==
"quadrature_total"
:
params
[
'sampler'
]
=
QST
(
murange
,
"CHEBYSHEV"
,
scalingExp
=
rescalingExp
)
else
:
# if sampling == "random":
params
[
'sampler'
]
=
RS
(
murange
,
"HALTON"
,
scalingExp
=
rescalingExp
)
elif
samples
==
"centered"
:
params
[
'sampler'
]
=
MS
(
murange
,
points
=
[
mu0
],
scalingExp
=
rescalingExp
)
elif
samples
==
"pivoted"
:
if
sampling
==
"quadrature"
:
params
[
'samplerPivot'
]
=
QS
([
murange
[
0
][
0
],
murange
[
1
][
0
]],
"CHEBYSHEV"
)
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif
sampling
==
"quadrature_total"
:
params
[
'samplerPivot'
]
=
QST
([
murange
[
0
][
0
],
murange
[
1
][
0
]],
"CHEBYSHEV"
)
else
:
# if sampling == "random":
params
[
'samplerPivot'
]
=
RS
([
murange
[
0
][
0
],
murange
[
1
][
0
]],
"HALTON"
)
if
samplingM
==
"quadrature"
:
params
[
'samplerMarginal'
]
=
QS
([
murange
[
0
][
1
],
murange
[
1
][
1
]],
"UNIFORM"
)
elif
samplingM
==
"quadrature_total"
:
params
[
'samplerMarginal'
]
=
QST
([
murange
[
0
][
1
],
murange
[
1
][
1
]],
"CHEBYSHEV"
)
else
:
# if samplingM == "random":
params
[
'samplerMarginal'
]
=
RS
([
murange
[
0
][
1
],
murange
[
1
][
1
]],
"HALTON"
)
if
samples
!=
"pivoted"
:
approx
=
method
(
solver
,
mu0
=
mu0
,
approxParameters
=
params
,
verbosity
=
verb
)
else
:
approx
=
method
(
solver
,
mu0
=
mu0
,
directionPivot
=
[
0
],
approxParameters
=
params
,
verbosity
=
verb
)
approx
.
setupApprox
()
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
)))
approx
.
verbosity
=
5
try
:
from
plot_zero_set
import
plotZeroSet2
muZeroVals
,
Qvals
=
plotZeroSet2
(
murange
,
murangeEff
,
approx
,
mu0
,
200
,
[
2.
,
1.
])
except
:
pass
if
show_norm
:
from
plot_inf_set
import
plotInfSet2
muInfVals
,
normEx
,
normApp
,
normRes
,
normErr
,
beta
=
plotInfSet2
(
murange
,
murangeEff
,
approx
,
mu0
,
25
,
[
2.
,
1.
],
relative
=
True
,
nobeta
=
True
)
try
:
print
(
np
.
sort
(
approx
.
getPoles
([
None
,
.
5
])
**
2.
))
except
:
pass
Event Timeline
Log In to Comment