diff --git a/examples/2d/pod/cookie_single_pod.py b/examples/2d/pod/cookie_single_pod.py
index eac23b1..d6d190b 100644
--- a/examples/2d/pod/cookie_single_pod.py
+++ b/examples/2d/pod/cookie_single_pod.py
@@ -1,137 +1,153 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
CookieEngineSingle as CES
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
+Delta = -10
MN = 15
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
+PODTol = 1e-6
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
+if samples == "distributed":
+ radial = 0
+# radial = "gaussian"
+# radial = "thinplate"
+# radial = "multiquadric"
+ rW0 = 1.
+ radialWeight = [rW0] * 2
+
+assert Delta <= 0
+
if size == 1: # below
mu0 = [20 ** .5, 1. ** .5]
mutar = [20.5 ** .5, 1.05 ** .5]
murange = [[18.5 ** .5, .85 ** .5], [21.5 ** .5, 1.15 ** .5]]
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]]]
kappa = 20. ** .5
theta = - np.pi / 6.
-n = 40
+n = 30
Rad = 1.
L = np.pi
nX = 2
nY = 1
solver = CES(kappa = kappa, theta = theta, n = n, R = Rad, L = L, nX = nX,
nY = nY, mu0 = mu0, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
- params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
+ params['radialBasis'] = radial
+ params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
- params = {'R':R, 'S':S, 'POD':True}
+ params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
+ 'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
what = "REAL")
approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
- murange, murangeEff, approx, mu0, 40, [2., 2.], clip = clip)
+ murange, murangeEff, approx, mu0, 25,
+ [2., 2.], clip = clip, relative = False)
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 9c02db0..1aff9fa 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,191 +1,194 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 3
+size = 4
show_sample = False
show_norm = True
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
Delta = 0
-MN = 10
+MN = 20
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
-PODTol = 1e-10
+PODTol = 1e-8
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
-sampling = "random"
+#sampling = "random"
if samples == "distributed":
radial = 0
# radial = "gaussian"
# radial = "thinplate"
# radial = "multiquadric"
- rW0 = 10.
- radialWeight = [r * rW0 for r in [40., .3]]
+ rW0 = 5.
+ radialWeight = [rW0] * 2
+
+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 == 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)
rescaling = [lambda x: np.power(x, 2.), lambda x: x]
rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
if algo == "rational":
- params = {'N':MN, 'M':MN - Delta, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
params['radialBasis'] = radial
params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
- params = {'R':(MN + 2 - Delta) * (MN + 1 - Delta) // 2, 'S':S,
+ params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
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',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
- 200, [2., 1.], clip = clip)
+ 50, [2., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
- murange, murangeEff, approx, mu0, 25, [2., 1.], clip = clip)
-
+ murange, murangeEff, approx, mu0, 50,
+ [2., 1.], clip = clip, relative = False,
+ normalizeDen = True)
diff --git a/examples/2d/pod/fracture_pod_nodomain.py b/examples/2d/pod/fracture_pod_nodomain.py
index 7764988..b9211f4 100644
--- a/examples/2d/pod/fracture_pod_nodomain.py
+++ b/examples/2d/pod/fracture_pod_nodomain.py
@@ -1,187 +1,190 @@
import numpy as np
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 3
+size = 1
show_sample = True
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
-M = 5
-N = 6
-R = max(M, N) + 1
+Delta = 0
+MN = 6
+R = MN + 1
S = R
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "distributed":
radial = 0
# radial = "gaussian"
# radial = "thinplate"
# radial = "multiquadric"
- radialWeight = np.array([2.])
+ radialWeight = [2.]
+
+assert Delta <= 0
if size == 1: # below
mu0Aug = [40 ** .5, .4]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 2: # top
mu0Aug = [40 ** .5, .6]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 3: # interesting
mu0Aug = [40 ** .5, .5]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 4: # wide_low
mu0Aug = [40 ** .5, .2]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 5: # wide_hi
mu0Aug = [40 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 6: # top_zoom
mu0Aug = [50 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 55 ** .5
murange = [[40 ** .5], [60 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta, n = n,
verbosity = verb)
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
if algo == "rational":
- params = {'N':N, 'M':M, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
- params['E'] = max(M, N)
+ params['E'] = MN
params['radialBasis'] = radial
params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
- params['S'] = [max(j, M + 1, N + 1) for j in params['S']]
+ params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = solver.lFrac
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',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet1
muZeroVals, Qvals = plotZeroSet1(murange, murangeEff, approx, mu0,
1000, 2.)
if show_norm:
solver._solveBatchSize = 10
from plot_inf_set import plotInfSet1
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet1(
- murange, murangeEff, approx, mu0, 200, 2.)
-
+ murange, murangeEff, approx, mu0,
+ 200, 2., relative = False,
+ normalizeDen = True)
diff --git a/examples/2d/pod/plot_inf_set.py b/examples/2d/pod/plot_inf_set.py
index 1e3b68c..8c01fc2 100644
--- a/examples/2d/pod/plot_inf_set.py
+++ b/examples/2d/pod/plot_inf_set.py
@@ -1,223 +1,323 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
-def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0, exp = 2.):
+def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
+ exp = 2., normalizeDen = False):
if hasattr(approx, "mus"):
mu2x = approx.mus(0) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
mu1 = np.real(np.power(mus, exp))
+ eta = R / beta / E
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
betamin, betamax = np.min(beta), np.max(beta)
+ etamin, etamax = np.min(eta), np.max(eta)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Z)
plt.semilogy(mu1, T, '--')
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [ZTmin, ZTmax], 'b:')
plt.plot(mu2x, [ZTmin] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|u(mu)|, |u_app(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, R)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:')
plt.plot(mu2x, [Rmax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
- plt.title("|res(mu)|")
+ if normalizeDen:
+ plt.title("|Q(mu)res(mu)|")
+ else:
+ plt.title("|res(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, E)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
- plt.title("|err(mu)|")
+ if normalizeDen:
+ plt.title("|Q(mu)err(mu)|")
+ else:
+ plt.title("|err(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, beta)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("beta(mu)")
plt.grid()
plt.show()
-def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.):
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ plt.semilogy(mu1, R / beta)
+ plt.semilogy(mu1, E, '--')
+ for l_ in approx.trainedModel.getPoles():
+ plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
+ plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
+ plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
+ plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
+ plt.xlim(mu1[0], mu1[-1])
+ if normalizeDen:
+ plt.title("|Q(mu)res(mu)/beta(mu)|")
+ else:
+ plt.title("|res(mu)/beta(mu)|")
+ plt.grid()
+ plt.show()
+
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ plt.semilogy(mu1, eta)
+ for l_ in approx.trainedModel.getPoles():
+ plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
+ plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
+ plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
+ plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
+ plt.xlim(mu1[0], mu1[-1])
+ plt.title("eta(mu)")
+ plt.grid()
+ plt.show()
+
+def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.,
+ relative = True, normalizeDen = False):
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
- F = approx.normRHS(mus)
Z = approx.normHF(mus)
T = approx.normApprox(mus)
- R = approx.normRes(mus) / F
- E = approx.normErr(mus) / Z
+ R = approx.normRes(mus)
+ E = approx.normErr(mus)
+ if relative:
+ F = approx.normRHS(mus)
+ R /= F
+ E /= Z
+ if normalizeDen:
+ Qvals = np.abs(approx.trainedModel.getQVal(mus))
+ R *= Qvals
+ E *= Qvals
beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
- plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0, exp)
+ plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
+ exp, normalizeDen)
return mus, Z, T, R, E, beta
def plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange, approx, mu0,
- exps = [2., 2.], clip = -1):
+ exps = [2., 2.], clip = -1, normalizeDen = False):
if hasattr(approx, "mus"):
mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
else:
mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
[murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
mu1s = np.unique([m[0] for m in mus])
mu2s = np.unique([m[1] for m in mus])
mu1 = np.power(mu1s, exps[0])
mu2 = np.power(mu2s, exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
Z = np.log10(Ze)
T = np.log10(Te)
R = np.log10(Re)
E = np.log10(Ee)
beta = np.log10(betae)
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
betamin, betamax = np.min(beta), np.max(beta)
if clip > 0:
ZTmax -= clip * (ZTmax - ZTmin)
cmap = plt.cm.bone
else:
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = UserWarning)
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Z, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Z, [ZTmin])
plt.plot(mu2x, mu2y, 'kx')
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.colorbar(p)
plt.title("log10|u(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, T, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, T, [ZTmin])
plt.plot(mu2x, mu2y, 'kx')
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|u_app(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R, levels = np.linspace(Rmin, Rmax, 50))
plt.plot(mu2x, mu2y, 'kx')
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(mu)|")
+ if normalizeDen:
+ plt.title("log10|Q(mu)res(mu)|")
+ else:
+ plt.title("log10|res(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, E, levels = np.linspace(Emin, Emax, 50))
plt.plot(mu2x, mu2y, 'kx')
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|err(mu)|")
+ if normalizeDen:
+ plt.title("log10|Q(mu)err(mu)|")
+ else:
+ plt.title("log10|err(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, beta,
levels = np.linspace(betamin, betamax, 50))
plt.plot(mu2x, mu2y, 'kx')
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(beta(mu))")
plt.colorbar(p)
plt.grid()
plt.show()
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ p = plt.contourf(Mu1, Mu2, R - beta, levels = np.linspace(Emin, Emax, 50))
+ plt.plot(mu2x, mu2y, 'kx')
+ 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:')
+ if normalizeDen:
+ plt.title("log10|Q(mu)res(mu)/beta(mu)|")
+ else:
+ plt.title("log10|res(mu)/beta(mu)|")
+ plt.colorbar(p)
+ plt.grid()
+ plt.show()
+
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ p = plt.contourf(Mu1, Mu2, R - beta - E, 50)
+ plt.plot(mu2x, mu2y, 'kx')
+ 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|eta(mu)|")
+ plt.colorbar(p)
+ plt.grid()
+ plt.show()
+
def plotInfSet2(murange, murangeEff, approx, mu0, nSamples = 200,
- exps = [2., 2.], clip = -1):
+ exps = [2., 2.], clip = -1, relative = True,
+ normalizeDen = False):
mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
nSamples)
mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
- Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
Ze = approx.normHF(mus).reshape((nSamples, nSamples))
Te = approx.normApprox(mus).reshape((nSamples, nSamples))
- Re = approx.normRes(mus).reshape((nSamples, nSamples)) / Fe
- Ee = approx.normErr(mus).reshape((nSamples, nSamples)) / Ze
+ Re = approx.normRes(mus).reshape((nSamples, nSamples))
+ Ee = approx.normErr(mus).reshape((nSamples, nSamples))
+ if relative:
+ Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
+ Re /= Fe
+ Ee /= Ze
+ if normalizeDen:
+ Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
+ (nSamples, nSamples)))
+ Re *= Qvals
+ Ee *= Qvals
betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus).reshape(
(nSamples, nSamples))
plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange,
- approx, mu0, exps, clip)
+ approx, mu0, exps, clip, normalizeDen)
return mus, Ze, Te, Re, Ee, betae
diff --git a/examples/2d/pod/synthetic_pod.py b/examples/2d/pod/synthetic_pod.py
index ce5ef62..a466a0d 100644
--- a/examples/2d/pod/synthetic_pod.py
+++ b/examples/2d/pod/synthetic_pod.py
@@ -1,157 +1,158 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
SyntheticBivariateEngine as SBE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 10
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
-PODTol = 1e-10
+PODTol = 1e-6
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
if samples == "distributed":
radial = 0
-# radial = "gaussian"
+ radial = "gaussian"
# radial = "thinplate"
# radial = "multiquadric"
rW0 = 10.
- radialWeight = [r * rW0 for r in [5., 5.]]
+ radialWeight = [rW0] * 2
if size == 1: # small
mu0 = [10. ** .5, 15. ** .5]
mutar = [12. ** .5, 14. ** .5]
murange = [[5. ** .5, 10. ** .5], [15 ** .5, 20 ** .5]]
if size == 2: # large
mu0 = [15. ** .5, 17.5 ** .5]
mutar = [18. ** .5, 22. ** .5]
murange = [[5. ** .5, 10. ** .5], [25 ** .5, 25 ** .5]]
if size == 3: # medium
mu0 = [17.5 ** .5, 15 ** .5]
mutar = [20. ** .5, 18. ** .5]
murange = [[10. ** .5, 10. ** .5], [25 ** .5, 20 ** .5]]
assert Delta <= 0
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]]]
kappa = 20. ** .5
theta = - np.pi / 6.
n = 20
L = np.pi
solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
mu0 = mu0, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
params['radialBasis'] = radial
params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
what = "REAL")
approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
- murange, murangeEff, approx, mu0, 50, [2., 2.], clip = clip)
+ murange, murangeEff, approx, mu0, 50,
+ [2., 2.], clip = clip, relative = False)
diff --git a/tests/hfengines/fracture.py b/tests/hfengines/fracture.py
index 55c3aaa..a7624a2 100644
--- a/tests/hfengines/fracture.py
+++ b/tests/hfengines/fracture.py
@@ -1,58 +1,59 @@
# 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 .
#
import pytest
import numpy as np
import ufl
import fenics as fen
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.solver.fenics import affine_warping
-@pytest.mark.xfail(raises = RuntimeError('Invalid DISPLAY variable'),
- reason = "no graphical interface")
+@pytest.mark.xfail(reason = "no graphical interface")
def test_fracture():
mu0 = [45. ** .5, .6]
solver2 = MFE(mu0 = mu0, H = 1., L = .75, delta = .05, n = 20,
verbosity = 0)
uh2 = solver2.solve(mu0)[0]
solver1 = MFEND(mu0 = mu0, H = 1., L = .75, delta = .05, n = 20,
verbosity = 0)
uh1 = solver1.solve(mu0[0])[0]
L = mu0[1]
y = fen.SpatialCoordinate(solver1.V.mesh())[1]
warp1, warpI1 = affine_warping(solver1.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver1.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)
solver1.plotmesh([warp, warpI], show = False, figsize = (7, 7))
solver1.plot(uh1, [warp, warpI], what = 'REAL', show = False)
from matplotlib import pyplot as plt
plt.close('all')
- assert np.isclose(solver1.norm(solver1.residual(uh1, mu0[0]))[0],
- solver2.norm(solver2.residual(uh2, mu0))[0], atol = 1e-5)
+ assert np.isclose(
+ solver1.norm(solver1.residual(uh1, mu0[0]), dual = True)[0],
+ solver2.norm(solver2.residual(uh2, mu0), dual = True)[0],
+ atol = 1e-5)
assert np.isclose(solver1.norm(uh1 - uh2), 0., atol = 1e-6)
diff --git a/tests/hfengines/helmholtz_elasticity.py b/tests/hfengines/helmholtz_elasticity.py
index 85a2b14..6d3279c 100644
--- a/tests/hfengines/helmholtz_elasticity.py
+++ b/tests/hfengines/helmholtz_elasticity.py
@@ -1,66 +1,69 @@
# 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 .
#
import numpy as np
from rrompy.hfengines.vector_linear_problem import (
LinearElasticityHelmholtzProblemEngine,
LinearElasticityHelmholtzProblemEngineDamped,
LinearElasticityHelmholtzArchwayFrequency)
from rod_3d import rod3Dsolver
def test_helmholtz_elastic_arch():
solver = LinearElasticityHelmholtzArchwayFrequency(kappa = 10, n = 30,
rho_ = 1e4, T = 1e5, lambda_ = 4e6, mu_ = 7e5,
R = 2e1, r = 1.5e1, verbosity = 0)
mu = 10
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 3188.9960782143194, rtol = 1e-5)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 3.025504915e-05,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 3.6771479e-12, rtol = 1e-1)
def test_helmholtz_elastic_rod():
solverBase = rod3Dsolver()
solver = LinearElasticityHelmholtzProblemEngine()
solver.V = solverBase.V
solver.lambda_ = solverBase.lambda_
solver.mu_ = solverBase.mu_
solver.forcingTerm = solverBase.forcingTerm
solver.DirichletBoundary = solverBase.DirichletBoundary
solver.NeumannBoundary = solverBase.NeumannBoundary
mu = 10
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 0.17836028624665373, rtol = 1e-5)
- assert np.isclose(solver.norm(solver.residual(uh, 10)[0]), 7.030048088e-08,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 8.070977e-07, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
def test_helmholtz_elastic_rod_damped():
solverBase = rod3Dsolver()
solver = LinearElasticityHelmholtzProblemEngineDamped()
solver.V = solverBase.V
solver.lambda_ = solverBase.lambda_
solver.mu_ = solverBase.mu_
solver.forcingTerm = solverBase.forcingTerm
solver.DirichletBoundary = solverBase.DirichletBoundary
solver.NeumannBoundary = solverBase.NeumannBoundary
solver.eta = 10
mu = 10
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 0.17646530119044376, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, 10)[0]), 6.802444e-08,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, 10)[0], dual = True),
+ 6.7057338e-07, rtol = 1e-1)
diff --git a/tests/hfengines/helmholtz_external.py b/tests/hfengines/helmholtz_external.py
index 009ab89..3e74d16 100644
--- a/tests/hfengines/helmholtz_external.py
+++ b/tests/hfengines/helmholtz_external.py
@@ -1,67 +1,69 @@
# 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 .
#
import pytest
import numpy as np
from rrompy.hfengines.linear_problem import (
HelmholtzCavityScatteringProblemEngine, HelmholtzBoxScatteringProblemEngine)
def test_helmholtz_square_scattering():
solver = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
n = 20, verbosity = 0)
mu = 5
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 19.9362, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 4.25056407e-13,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 4.25056407e-13, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
def test_helmholtz_scattering_copy(capsys):
solver1 = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
n = 20, verbosity = 0)
mu = 5
uh1 = solver1.solve(mu)[0]
solver2 = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
n = 20, verbosity = 100)
assert solver1.As[0] is not None and solver1.bs[0] is not None
assert solver2.As[0] is None and solver2.bs[0] is None
solver2.setAs(solver1.As)
solver2.setbs(solver1.bs)
uh2 = solver2.solve(mu)[0]
assert np.allclose(uh1, uh2, rtol = 1e-8)
out, err = capsys.readouterr()
assert ("Assembling operator term" not in out
and "Assembling forcing term" not in out)
assert len(err) == 0
-@pytest.mark.xfail(raises = RuntimeError('Invalid DISPLAY variable'),
- reason = "no graphical interface")
+@pytest.mark.xfail(reason = "no graphical interface")
def test_helmholtz_box_scattering():
solver = HelmholtzBoxScatteringProblemEngine(R = 2, kappa = 10.,
theta = np.pi * 30 / 180, n = 20, verbosity = 0)
mu = 15
uh = solver.solve(mu)[0]
solver.plotmesh(show = False, figsize = (7, 7))
assert np.isclose(solver.norm(uh), 62.113, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 9.62989935e-13,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 9.62989935e-13, rtol = 1e-1)
from matplotlib import pyplot as plt
plt.close('all')
diff --git a/tests/hfengines/helmholtz_internal.py b/tests/hfengines/helmholtz_internal.py
index 00724bb..920dfcb 100644
--- a/tests/hfengines/helmholtz_internal.py
+++ b/tests/hfengines/helmholtz_internal.py
@@ -1,98 +1,99 @@
# 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 .
#
import pytest
import os, shutil
import numpy as np
from rrompy.hfengines.linear_problem import (
HelmholtzSquareBubbleDomainProblemEngine, HelmholtzSquareBubbleProblemEngine,
HelmholtzSquareTransmissionProblemEngine)
def test_helmholtz_square_io():
- solver = HelmholtzSquareBubbleProblemEngine(kappa = 4, theta = 1., n = 50,
+ solver = HelmholtzSquareBubbleProblemEngine(kappa = 4, theta = 1., n = 20,
verbosity = 0)
mu = 5
uh = solver.solve(mu)[0]
- assert np.isclose(solver.norm(uh), 913.396, rtol = 1e-3)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 1.19934e-11,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(uh), 145.0115, rtol = 1e-3)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 1.19934e-11, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".pvd" and x[:9] == "outSquare")]
filesOutData = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".vtu" and x[:9] == "outSquare")]
for fileOut in filesOut:
os.remove("./.pytest_cache/" + fileOut)
for fileOut in filesOutData:
os.remove("./.pytest_cache/" + fileOut)
solver.outParaview(uh, what = ["MESH", "ABS"],
filename = ".pytest_cache/outSquare",
forceNewFile = False)
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".pvd" and x[:9] == "outSquare")]
filesOutData = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".vtu" and x[:9] == "outSquare")]
assert len(filesOut) == 1
assert len(filesOutData) == 1
os.remove("./.pytest_cache/" + filesOut[0])
os.remove("./.pytest_cache/" + filesOutData[0])
def test_helmholtz_transmission_io():
solver = HelmholtzSquareTransmissionProblemEngine(nT = 1, nB = 2,
- theta = np.pi * 40 / 180, kappa = 4., n = 50, verbosity = 0)
+ theta = np.pi * 40 / 180, kappa = 4., n = 20, verbosity = 0)
mu = 5.
uh = solver.solve(mu)[0]
- assert np.isclose(solver.norm(uh), 43.9268, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 3.7288565e-12,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(uh), 138.6609, rtol = 1e-2)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 3.7288565e-12, rtol = 1e-1)
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
solver.outParaviewTimeDomain(uh, omega = mu,
filename = ".pytest_cache/outTrans",
forceNewFile = False, folder = True)
filesOut = [x for x in os.listdir("./.pytest_cache/outTrans") if
(x[-4:] == ".pvd" and x[:8] == "outTrans")]
filesOutData = [x for x in os.listdir("./.pytest_cache/outTrans") if
(x[-4:] == ".vtu" and x[:8] == "outTrans")]
assert len(filesOut) == 1
assert len(filesOutData) == 20
shutil.rmtree("./.pytest_cache/outTrans")
-@pytest.mark.xfail(raises = RuntimeError('Invalid DISPLAY variable'),
- reason = "no graphical interface")
+@pytest.mark.xfail(reason = "no graphical interface")
def test_helmholtz_domain_io():
solver = HelmholtzSquareBubbleDomainProblemEngine(kappa = 4, theta = 1.,
n = 10, mu0 = 1.5, verbosity = 0)
mu = 1.5
uh = solver.solve(mu)[0]
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
solver.plot(uh, save = "./.pytest_cache/outDomain", show = False)
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".eps" and x[:9] == "outDomain")]
assert len(filesOut) == 1
os.remove("./.pytest_cache/" + filesOut[0])
- assert np.isclose(solver.norm(uh), 8.9947, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 6.14454989e-13,
- rtol = 1e-1)
-
+ assert np.isclose(solver.norm(uh), 10.07843, rtol = 1e-2)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 6.14454989e-13, rtol = 1e-1)
diff --git a/tests/hfengines/laplace.py b/tests/hfengines/laplace.py
index 8bd21f7..a76dbc1 100644
--- a/tests/hfengines/laplace.py
+++ b/tests/hfengines/laplace.py
@@ -1,39 +1,42 @@
# 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 .
#
import numpy as np
from rrompy.hfengines.linear_problem import LaplaceDiskGaussian
from rrompy.hfengines.linear_problem.bidimensional import LaplaceDiskGaussian2
def test_laplace_disk():
solver = LaplaceDiskGaussian(n = 20, verbosity = 0)
mu = 1.5
solver.setSolver("BICG", {"tol" : 1e-15})
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 1.053403077447029, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 5.27345e-13,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 3.4591353e-08, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
def test_laplace_disk_2():
solver = LaplaceDiskGaussian2(n = 20, verbosity = 0)
mu = [[0., 1.5]]
mu = [0., 1.5]
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 1.0534030774205372, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 5.27345e-13,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 2.40043363e-08, rtol = 1e-1)
diff --git a/tests/hfengines/linear_elasticity.py b/tests/hfengines/linear_elasticity.py
index 92e10da..309716a 100644
--- a/tests/hfengines/linear_elasticity.py
+++ b/tests/hfengines/linear_elasticity.py
@@ -1,38 +1,41 @@
# 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 .
#
import numpy as np
from rrompy.hfengines.vector_linear_problem import (
LinearElasticityBeamPoissonRatio)
from rod_3d import rod3Dsolver
def test_elastic_beam():
solver = LinearElasticityBeamPoissonRatio(n = 10, rho_ = 1e3, g = 3,
E = 1e6, nu0 = .45, length = 5, verbosity = 0)
mu = .45
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 5.866888537228743e-08, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 8.4545952e-13,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 8.4545952e-13, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
def test_elastic_rod():
solver = rod3Dsolver()
uh = solver.solve()[0]
assert np.isclose(solver.norm(uh), 0.15563476339534466, rtol = 1e-5)
- assert np.isclose(solver.norm(solver.residual(uh)[0]), 5.708389944e-08,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh)[0], dual = True),
+ 1.2210129e-07, rtol = 1e-1)
diff --git a/tests/hfengines/matrix.py b/tests/hfengines/matrix.py
index db2f54a..d0959e5 100644
--- a/tests/hfengines/matrix.py
+++ b/tests/hfengines/matrix.py
@@ -1,62 +1,65 @@
# 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 .
#
import numpy as np
import scipy.sparse as sp
from rrompy.hfengines.base import MatrixEngineBase as MEB
def test_deterministic():
N = 100
verb = 0
solver = MEB(verbosity = verb)
solver.npar = 1
solver.nAs = 2
mu = 10. + .5j
solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- sp.eye(N)]
solver.nbs = 1
solver.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
uh = solver.solve(mu)[0]
- assert np.isclose(np.linalg.norm(solver.residual(uh, mu)[0]), 1.088e-15,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 1.088e-15, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ solver.norm(solver.residual(uh, mu, duality = False)[0],
+ dual = True, duality = False), rtol = 1e-1)
def test_random():
N = 100
verb = 0
solver = MEB(verbosity = verb)
solver.npar = 1
solver.nAs = 2
mu = 1. + .5j
np.random.seed(420)
solver.setSolver("SOLVE")
fftB = np.fft.fft(np.eye(N)) * N**-.5
solver.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
- np.eye(N)]
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
uh = solver.solve(mu)[0]
- assert np.isclose(np.linalg.norm(solver.residual(uh, mu)[0]), 7.18658e-14,
- rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ 7.18658e-14, rtol = 1e-1)
diff --git a/tests/utilities/radial_fitting.py b/tests/utilities/radial_fitting.py
index e1cd8c7..766e0a9 100644
--- a/tests/utilities/radial_fitting.py
+++ b/tests/utilities/radial_fitting.py
@@ -1,128 +1,161 @@
# 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 .
#
import numpy as np
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.poly_fitting.radial_basis import (radialGaussian,
thinPlateSpline,
multiQuadric,
polybases, polyfitname,
polydomcoeff,
radialFunction,
- polyval, polyvander)
+ polyval, polyvander,
+ homogeneizedpolyvander)
+from rrompy.utilities.poly_fitting.polynomial import homogeneizedToFull
from rrompy.parameter import checkParameterList
def test_monomial_gaussian():
polyrbname = "MONOMIAL_GAUSSIAN"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "polyfit_gaussian"
assert np.isclose(domcoeff, 1., rtol = 1e-5)
directionalWeights = np.array([5.])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
cRB = radialFunction(directionalWeights = directionalWeights,
supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
globalCoeffs = cRBCoeffs[4 :])
ySupp = 1 + 2. * xSupp.data - .5 * xSupp.data ** 2.
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
yyman = 1. + 2. * xx - .5 * xx ** 2.
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (5. * (xx - xc)) ** 2.
rbj = radialGaussian(r2j)
assert np.allclose(rbj, np.exp(-.5 * r2j))
yyman += cRB.localCoeffs[j] * rbj
ySupp += cRB.localCoeffs[j] * radialGaussian((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
def test_legendre_thinplate():
polyrbname = "LEGENDRE_THINPLATE"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "legfit_thinplate"
assert np.isclose(domcoeff, 63. / 8, rtol = 1e-5)
directionalWeights = np.array([.5])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
cRB = radialFunction(directionalWeights = directionalWeights,
supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
globalCoeffs = cRBCoeffs[4 :])
ySupp = 1 + 2. * xSupp.data - .5 * (.5 * (3. * xSupp.data ** 2. - 1.))
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
yyman = 1. + 2. * xx - .5 * (.5 * (3. * xx ** 2. - 1.))
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (directionalWeights[0] * (xx - xc)) ** 2.
rbj = thinPlateSpline(r2j)
assert np.allclose(rbj, .5 * r2j * np.log(np.finfo(float).eps + r2j))
yyman += cRB.localCoeffs[j] * rbj
ySupp += cRB.localCoeffs[j] * thinPlateSpline((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
def test_chebyshev_multiquadric():
polyrbname = "CHEBYSHEV_MULTIQUADRIC"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "chebfit_multiquadric"
assert np.isclose(domcoeff, 16, rtol = 1e-5)
directionalWeights = np.array([1.])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
cRB = radialFunction(directionalWeights = directionalWeights,
supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
globalCoeffs = cRBCoeffs[4 :])
ySupp = 1 + 2. * xSupp.data - .5 * (2. * xSupp.data ** 2. - 1.)
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
yyman = 1. + 2. * xx - .5 * (2. * xx ** 2. - 1.)
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (directionalWeights[0] * (xx - xc)) ** 2.
rbj = multiQuadric(r2j)
assert np.allclose(rbj, np.power(r2j + 1, -.5))
yyman += cRB.localCoeffs[j] * rbj
ySupp += cRB.localCoeffs[j] * multiQuadric((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
+def test_total_degree_2d():
+ values = lambda x, y: (x - 3.) ** 2. * y - (x + 1.) * y ** 2.
+
+ polyrbname = "CHEBYSHEV_GAUSSIAN"
+
+ xs, ys = np.meshgrid(np.linspace(0., 4., 5), np.linspace(0., 4., 4))
+ xySupp = np.concatenate((xs.reshape(-1, 1), ys.reshape(-1, 1)), axis = 1)
+ zs = values(xs, ys)
+ zSupp = zs.flatten()
+
+ deg = 3
+ directionalWeights = [2., 1.]
+ VanT, _, reidxs = homogeneizedpolyvander(xySupp, deg, polyrbname,
+ directionalWeights = directionalWeights)
+ VanT = VanT[reidxs]
+ VanT = VanT[:, reidxs]
+ cFit = np.linalg.solve(VanT, np.pad(zSupp, (0, len(VanT) - len(zSupp)),
+ 'constant'))
+ globCoeff = homogeneizedToFull([deg + 1] * 2, 2, cFit[len(zSupp) :])
+ cRB = radialFunction(directionalWeights = directionalWeights,
+ supportPoints = xySupp,
+ localCoeffs = cFit[: len(zSupp)],
+ globalCoeffs = globCoeff)
+ xx, yy = np.meshgrid(np.linspace(0., 4., 100), np.linspace(0., 4., 100))
+ xxyy = np.concatenate((xx.reshape(-1, 1), yy.reshape(-1, 1)), axis = 1)
+ zz = polyval(xxyy, cRB, polyrbname).reshape(xx.shape)
+ zzex = values(xx, yy)
+
+ error = np.abs(zz - zzex)
+ print(np.max(error))
+ assert np.max(error) < 1e-10