diff --git a/examples/2d/base/cookie_single.py b/examples/2d/base/cookie_single.py
deleted file mode 100644
index f153e01..0000000
--- a/examples/2d/base/cookie_single.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import numpy as np
-import fenics as fen
-from rrompy.hfengines.linear_problem.bidimensional import \
- CookieEngineSingle as CES
-
-verb = 100
-
-kappa = 15. ** .5
-theta = - np.pi / 6.
-n = 30
-R = 1.
-L = np.pi
-nX = 3
-nY = 2
-
-mu0 = [25. ** .5, 1.]
-mutar = [25. ** .5, 1.5]
-
-solver = CES(kappa = kappa, theta = theta, n = n, R = R, L = L, nX = nX,
- nY = nY, mu0 = mu0, verbosity = verb)
-uh = solver.solve(mutar)[0]
-solver.plotmesh(figsize = (7.5, 4.5))
-fen.plot(fen.project(solver.cookieIn,
- fen.FunctionSpace(solver.V.mesh(), "DG", 0)))
-print(solver.norm(uh))
-solver.plot(uh)
-solver.outParaviewTimeDomain(uh, mutar[0], filename = 'out', folder = True,
- forceNewFile = False)
-
diff --git a/examples/2d/base/fracture.py b/examples/2d/base/fracture.py
index b98b5e5..8c660db 100644
--- a/examples/2d/base/fracture.py
+++ b/examples/2d/base/fracture.py
@@ -1,46 +1,45 @@
import numpy as np
import ufl
import fenics as fen
-from rrompy.hfengines.linear_problem.bidimensional import \
- MembraneFractureEngine as MFE
+from membrane_fracture_engine import MembraneFractureEngine as MFE
from rrompy.solver.fenics import affine_warping
verb = 100
mu0 = [45. ** .5, .7]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
#solver.plotmesh(figsize = (7.5, 4.5))
#solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
#solver.plot(uh, what = 'REAL', figsize = (8, 5))
#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], filename = 'out', folder = True,
# forceNewFile = False)
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * mu0[1]]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * mu0[1]]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
#solver.plotmesh([warp, warpI], figsize = (7.5, 4.5))
#solver.plot(u0, [warp, warpI], what = 'REAL', figsize = (8, 5))
solver.plot(uh, [warp, warpI], what = 'REAL', figsize = (8, 5))
#solver.plot(solver.residual(mu0, uh)[0], [warp, warpI], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], [warp, warpI],
# filename = 'outW', folder = True,
# forceNewFile = False)
diff --git a/examples/2d/base/fracture_nodomain.py b/examples/2d/base/fracture_nodomain.py
index 400c730..ffbc047 100644
--- a/examples/2d/base/fracture_nodomain.py
+++ b/examples/2d/base/fracture_nodomain.py
@@ -1,47 +1,47 @@
import numpy as np
import ufl
import fenics as fen
-from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
+from membrane_fracture_nodomain_engine import MembraneFractureNoDomainEngine \
as MFEND
from rrompy.solver.fenics import affine_warping
verb = 100
mu0Aug = [45. ** .5, .6]
mu0Aug = [45. ** .5, .1]
mu0 = mu0Aug[0]
H = 1.
L = .75
delta = .05
n = 50
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
solver.plotmesh(figsize = (7.5, 4.5))
solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
solver.plot(uh, what = 'REAL', figsize = (8, 5))
solver.plot(solver.residual(mu0, uh)[0], name = 'res',
what = 'REAL', figsize = (8, 5))
solver.outParaviewTimeDomain(uh, mu0, filename = 'outND', folder = True)
##
L = mu0Aug[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)
solver.plotmesh([warp, warpI], figsize = (7.5, 4.5))
solver.plot(u0, [warp, warpI], what = 'REAL', figsize = (8, 5))
solver.plot(uh, [warp, warpI], what = 'REAL', figsize = (8, 5))
solver.plot(solver.residual(mu0, uh)[0], [warp, warpI], name = 'res',
what = 'REAL', figsize = (8, 5))
solver.outParaviewTimeDomain(uh, mu0, [warp, warpI],
filename = 'outNDW', folder = True)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py b/examples/2d/base/helmholtz_square_domain_problem_engine.py
similarity index 95%
rename from rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py
rename to examples/2d/base/helmholtz_square_domain_problem_engine.py
index 3626d52..c12868a 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py
+++ b/examples/2d/base/helmholtz_square_domain_problem_engine.py
@@ -1,122 +1,122 @@
# 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 numpy.polynomial.polynomial import polyfit as fit
import fenics as fen
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
-__all__ = ['HelmholtzSquareDomainProblemEngine']
-
class HelmholtzSquareDomainProblemEngine(HelmholtzProblemEngine):
"""
Solver for square Helmholtz problems with parametric laplacian.
- \dxx u - mu_2**2 \dyy u - mu_1**2 * u = f(mu_2) in \Omega = [0,\pi]^2
u = 0 on \partial\Omega
"""
def __init__(self, kappa:float, theta:float, n:int,
mu0 : paramVal = [12. ** .5, 1.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = False
self.nAs, self.nbs = 3, 11
self.npar = 2
self.rescalingExp = [2., 1.]
self._theta = theta
self._kappa = kappa
pi = np.pi
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
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, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
self.thAs[0] = self.getMonomialSingleWeight([0, 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, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
self.thAs[1] = self.getMonomialSingleWeight([1, 0])
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, fenZERO,
self.DirichletBoundary)
a2Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
self.thAs[2] = self.getMonomialSingleWeight([0, 2])
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
bDEIMCoeffs = None
for j in range(self.nbs):
if self.bs[j] is None:
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
if bDEIMCoeffs is None:
- bDEIM = np.empty((self.nbs, self.spacedim),
- dtype = np.complex)
muDEIM = np.linspace(.5, 4., bDEIM.shape[0])
for jj, muD in enumerate(muDEIM):
c, s = np.cos(self._theta), np.sin(self._theta)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
fRe = fen.cos(self._kappa * (c * x + s / muD * y))
fIm = fen.sin(self._kappa * (c * x + s / muD * y))
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTerm{}Real".format(jj)]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTerm{}Imag".format(jj)]))
LR = fen.dot(fRe, self.v) * fen.dx
LI = fen.dot(fIm, self.v) * fen.dx
DBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1)
+ bjj = (fenics2Vector(LR, parsRe, DBC0, 1)
+ 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
+ if jj == 0:
+ bDEIM = np.empty((self.nbs, len(bjj)),
+ dtype = np.complex)
+ bDEIM[jj] = bjj
bDEIMCoeffs = fit(muDEIM, bDEIM, len(muDEIM) - 1)
self.bs[j] = bDEIMCoeffs[j]
self.thbs[j] = self.getMonomialSingleWeight([0, j])
vbMng(self, "DEL", "Done assembling forcing term.", 20)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py b/examples/2d/base/helmholtz_square_simplified_domain_problem_engine.py
similarity index 98%
rename from rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py
rename to examples/2d/base/helmholtz_square_simplified_domain_problem_engine.py
index 967952f..0022352 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py
+++ b/examples/2d/base/helmholtz_square_simplified_domain_problem_engine.py
@@ -1,88 +1,86 @@
# 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 fenics as fen
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['HelmholtzSquareSimplifiedDomainProblemEngine']
-
class HelmholtzSquareSimplifiedDomainProblemEngine(HelmholtzProblemEngine):
"""
Solver for square Helmholtz problems with parametric laplacian.
- \dxx u - mu_2**2 \dyy u - mu_1**2 * u = f in \Omega_mu = [0,\pi]^2
u = 0 on \partial\Omega
"""
def __init__(self, kappa:float, theta:float, n:int,
mu0 : paramVal = [12. ** .5, 1.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.nAs = 3
self.npar = 2
self.rescalingExp = [2., 2.]
pi = np.pi
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
fen.sin(kappa * (c * x + s * y))]
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
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, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
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, fenZERO,
self.DirichletBoundary)
a2Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/rrompy/hfengines/linear_problem/tridimensional/matrix_dynamical_passive.py b/examples/2d/base/matrix_dynamical_passive_engine.py
similarity index 84%
rename from rrompy/hfengines/linear_problem/tridimensional/matrix_dynamical_passive.py
rename to examples/2d/base/matrix_dynamical_passive_engine.py
index fcec560..33e2095 100644
--- a/rrompy/hfengines/linear_problem/tridimensional/matrix_dynamical_passive.py
+++ b/examples/2d/base/matrix_dynamical_passive_engine.py
@@ -1,61 +1,60 @@
# 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.utilities.base.types import ListAny, paramVal
-from rrompy.hfengines.base.matrix_engine_base import MatrixEngineBase
+from rrompy.hfengines.base.linear_affine_engine import LinearAffineEngine
+from rrompy.hfengines.base.numpy_engine_base import NumpyEngineBase
from rrompy.parameter import checkParameter
-__all__ = ['MatrixDynamicalPassive']
-
-class MatrixDynamicalPassive(MatrixEngineBase):
+class MatrixDynamicalPassiveEngine(LinearAffineEngine, NumpyEngineBase):
def __init__(self, mu0 : paramVal = [0., 10.], n : int = 1000,
omega0 : ListAny = [1.], domega0 : ListAny = [1.],
b : float = 10., verbosity : int = 10,
timestamp : bool = True):
super().__init__(verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.npar = 1 + len(omega0)
self.mu0 = checkParameter(mu0, self.npar)
self.nAs, self.nbs = 1 + self.npar, 1
Asize = 2 + 2 * self.npar + n
dataA0 = np.concatenate((*tuple([tuple([np.imag(o), - np.real(o),
np.real(o), np.imag(o)]) \
for o in omega0]),
[1., -200., 200., 1., 1., -400., 400., 1.],
1. + np.arange(n)))
rowA0 = np.concatenate((np.repeat(np.arange(2 * len(omega0) + 4), 2),
2 + 2 * self.npar + np.arange(n)))
cB0 = np.repeat(np.arange(0, 2 * len(omega0) + 4, 2), 4)
cB1 = np.tile([0, 1], [1, 2 * len(omega0) + 4]).reshape(-1)
colA0 = np.concatenate((cB0 + cB1, 2 + 2 * self.npar + np.arange(n)))
self.As = [sp.csr_matrix((dataA0, (rowA0, colA0)), dtype = np.complex,
- shape = (Asize, Asize))]
+ shape = (Asize, Asize))]
self.As = self.As + [1.j * sp.eye(Asize, dtype = np.complex)]
for j, do0 in enumerate(domega0):
self.As = self.As + [sp.csr_matrix(([- do0, do0],
([2*j, 2*j+1], [2*j+1, 2*j])),
dtype = np.complex,
shape = (Asize, Asize))]
- self.bs = [np.concatenate((b * np.ones(2 + 2 * self.npar),
- np.ones(n)))]
- self.thAs = self.getMonomialWeights(self.nAs)
- self.thbs = self.getMonomialWeights(self.nbs)
+ self.setbs([np.concatenate((b * np.ones(2 + 2 * self.npar),
+ np.ones(n)))])
+ self.setthAs(self.getMonomialWeights(self.nAs))
+ self.setthbs(self.getMonomialWeights(self.nbs))
diff --git a/examples/2d/base/matrix_passive.py b/examples/2d/base/matrix_passive.py
index cc1331a..abbbf8e 100644
--- a/examples/2d/base/matrix_passive.py
+++ b/examples/2d/base/matrix_passive.py
@@ -1,13 +1,12 @@
-from rrompy.hfengines.linear_problem.tridimensional import \
- MatrixDynamicalPassive as MDP
+from matrix_dynamical_passive_engine import MatrixDynamicalPassiveEngine as MDP
verb = 100
n = 100
b = 10
mu0 = [0., 10]
mutar = [4., 15]
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
uh = solver.solve(mutar)[0]
solver.plot(uh)
diff --git a/examples/2d/base/membrane_fracture_engine.py b/examples/2d/base/membrane_fracture_engine.py
new file mode 120000
index 0000000..9e382c7
--- /dev/null
+++ b/examples/2d/base/membrane_fracture_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/membrane_fracture_engine.py
\ No newline at end of file
diff --git a/examples/2d/base/membrane_fracture_nodomain_engine.py b/examples/2d/base/membrane_fracture_nodomain_engine.py
new file mode 120000
index 0000000..45ce4c8
--- /dev/null
+++ b/examples/2d/base/membrane_fracture_nodomain_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/membrane_fracture_nodomain_engine.py
\ No newline at end of file
diff --git a/examples/2d/base/square.py b/examples/2d/base/square.py
index 546b09c..139831e 100644
--- a/examples/2d/base/square.py
+++ b/examples/2d/base/square.py
@@ -1,13 +1,14 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_domain_problem_engine import \
HelmholtzSquareDomainProblemEngine as HSDPE
+
verb = 100
mu0 = [4 ** .5, 1.5 ** .5]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 50,
verbosity = verb)
uh = solver.solve(mu0)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
solver.plot(solver.residual(mu0, uh)[0], name = 'res')
diff --git a/examples/2d/base/square_simplified.py b/examples/2d/base/square_simplified.py
index 4a05610..7e7b836 100644
--- a/examples/2d/base/square_simplified.py
+++ b/examples/2d/base/square_simplified.py
@@ -1,13 +1,13 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_simplified_domain_problem_engine import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
verb = 0
mu0 = [4 ** .5, 1.5 ** .5]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 50,
verbosity = verb)
uh = solver.solve(mu0)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
solver.plot(solver.residual(mu0, uh)[0], name = 'res')
diff --git a/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py b/examples/2d/base/synthetic_bivariate_engine.py
similarity index 99%
rename from rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
rename to examples/2d/base/synthetic_bivariate_engine.py
index 1d10c7e..7758a90 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
+++ b/examples/2d/base/synthetic_bivariate_engine.py
@@ -1,107 +1,105 @@
# 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 fenics as fen
import ufl
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenONE, fenZERO
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['SyntheticBivariateEngine']
-
class SyntheticBivariateEngine(HelmholtzProblemEngine):
def __init__(self, kappa:float, theta:float, n:int, L : int = 2.,
mu0 : paramVal = [12. ** .5, 15. ** .5],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.nAs = 3
self.npar = 2
self.rescalingExp = [2., 2.]
mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(L, L), n, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
x, y = fen.SpatialCoordinate(mesh)[:]
self._above = ufl.conditional(ufl.ge(y, .5 * L), fenONE, fenZERO)
self._below = fenONE - self._above
c, s = np.cos(theta), np.sin(theta)
self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
fen.sin(kappa * (c * x + s * y))]
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
termNames = ["diffusivity", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[aRe, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[aIm, hIm],
[x + "Imag" for x in termNames]))
a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ hRe * fen.dot(self.u, self.v) * self.ds(1))
a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ hIm * fen.dot(self.u, self.v) * self.ds(1))
self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 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, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * self._above * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * self._above * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
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, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a2Re = - n2Re * self._below * fen.dot(self.u, self.v) * fen.dx
a2Im = - n2Im * self._below * fen.dot(self.u, self.v) * fen.dx
self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/examples/2d/base/synthetic_solve.py b/examples/2d/base/synthetic_solve.py
index 6e58ce7..9c574e3 100644
--- a/examples/2d/base/synthetic_solve.py
+++ b/examples/2d/base/synthetic_solve.py
@@ -1,26 +1,25 @@
import numpy as np
import fenics as fen
-from rrompy.hfengines.linear_problem.bidimensional import \
- SyntheticBivariateEngine as SBE
+from synthetic_bivariate_engine import SyntheticBivariateEngine as SBE
verb = 100
kappa = 15. ** .5
theta = - np.pi / 6.
n = 30
L = np.pi
mu0 = [15. ** .5, 20. ** .5]
mutar = [15. ** .5, 20. ** .5]
solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
mu0 = mu0, verbosity = verb)
uh = solver.solve(mutar)[0]
solver.plotmesh(figsize = (7.5, 4.5))
fen.plot(fen.project(solver._above,
fen.FunctionSpace(solver.V.mesh(), "DG", 0)))
print(solver.norm(uh))
solver.plot(uh)
#solver.outParaviewTimeDomain(uh, mutar[0], filename = 'out', folder = True,
# forceNewFile = False)
diff --git a/examples/2d/greedy/fracture_pod.py b/examples/2d/greedy/fracture_pod.py
index 1732615..f752364 100644
--- a/examples/2d/greedy/fracture_pod.py
+++ b/examples/2d/greedy/fracture_pod.py
@@ -1,221 +1,220 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
- MembraneFractureEngine as MFE
+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{}\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)))
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
diff --git a/examples/2d/greedy/membrane_fracture_engine.py b/examples/2d/greedy/membrane_fracture_engine.py
new file mode 120000
index 0000000..9e382c7
--- /dev/null
+++ b/examples/2d/greedy/membrane_fracture_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/membrane_fracture_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/cookie_single_pod.py b/examples/2d/pod/cookie_single_pod.py
deleted file mode 100644
index 98a4f3c..0000000
--- a/examples/2d/pod/cookie_single_pod.py
+++ /dev/null
@@ -1,124 +0,0 @@
-import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
- CookieEngineSingle as CES
-from rrompy.reduction_methods.standard import RationalInterpolant as RI
-from rrompy.reduction_methods.standard import ReducedBasis as RB
-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 = 5
-R = (MN + 2) * (MN + 1) // 2
-STensorized = (MN + 1) ** 2
-PODTol = 1e-6
-
-samples = "centered"
-samples = "standard"
-algo = "rational"
-#algo = "RB"
-sampling = "quadrature"
-sampling = "quadrature_total"
-sampling = "random"
-
-if samples == "standard":
- radial = ""
-# 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 = 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)
-
-rescalingExp = [2.] * 2
-if algo == "rational":
- params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
- if samples == "standard":
- params['polybasis'] = "CHEBYSHEV" + radial
-# params['polybasis'] = "LEGENDRE" + radial
-# params['polybasis'] = "MONOMIAL" + radial
- params['radialDirectionalWeights'] = radialWeight
- elif samples == "centered":
- params['polybasis'] = "MONOMIAL"
- method = RI
-else: #if algo == "RB":
- params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':R,
- 'POD':True, 'PODTolerance':PODTol}
- method = RB
-
-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)
-
-approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
-
-approx.setupApprox()
-if show_sample:
- L = mutar[1]
- approx.plotApprox(mutar, name = 'u_app', what = "REAL")
- approx.plotHF(mutar, name = 'u_HF', what = "REAL")
- approx.plotErr(mutar, name = 'err', what = "REAL")
-# approx.plotRes(mutar, name = 'res', 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:
- from plot_inf_set import plotInfSet2
- muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
- 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 c32c8a7..24a2046 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,265 +1,264 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
- MembraneFractureEngine as MFE
+from membrane_fracture_engine 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{}\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)))
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
diff --git a/examples/2d/pod/fracture_pod_nodomain.py b/examples/2d/pod/fracture_pod_nodomain.py
index 913786a..da958d7 100644
--- a/examples/2d/pod/fracture_pod_nodomain.py
+++ b/examples/2d/pod/fracture_pod_nodomain.py
@@ -1,155 +1,155 @@
import numpy as np
-from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
+from membrane_fracture_nodomain_engine import MembraneFractureNoDomainEngine \
as MFEND
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
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
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 6
R = MN + 1
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
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)
rescalingExp = 2.
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
method = RB
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)
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)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
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', 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{}\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:
from plot_inf_set import plotInfSet1
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet1(
murange, murangeEff, approx, mu0,
200, 2., relative = False,
normalizeDen = True)
print(approx.getPoles() ** 2.)
diff --git a/examples/2d/pod/helmholtz_square_domain_problem_engine.py b/examples/2d/pod/helmholtz_square_domain_problem_engine.py
new file mode 120000
index 0000000..fc97c18
--- /dev/null
+++ b/examples/2d/pod/helmholtz_square_domain_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/2d/base/helmholtz_square_domain_problem_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/helmholtz_square_simplified_domain_problem_engine.py b/examples/2d/pod/helmholtz_square_simplified_domain_problem_engine.py
new file mode 120000
index 0000000..1d2e2e0
--- /dev/null
+++ b/examples/2d/pod/helmholtz_square_simplified_domain_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/2d/base/helmholtz_square_simplified_domain_problem_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/matrix_dynamical_passive_engine.py b/examples/2d/pod/matrix_dynamical_passive_engine.py
new file mode 120000
index 0000000..baf88ed
--- /dev/null
+++ b/examples/2d/pod/matrix_dynamical_passive_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/2d/base/matrix_dynamical_passive_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/matrix_passive_pod.py b/examples/2d/pod/matrix_passive_pod.py
index 294c429..9b70b6d 100644
--- a/examples/2d/pod/matrix_passive_pod.py
+++ b/examples/2d/pod/matrix_passive_pod.py
@@ -1,191 +1,190 @@
import numpy as np
-from rrompy.hfengines.linear_problem.tridimensional import \
- MatrixDynamicalPassive as MDP
+from matrix_dynamical_passive_engine import MatrixDynamicalPassiveEngine as MDP
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 = 3
show_sample = True
show_norm = True
Delta = 0
MN = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
SMarginal = 3
MMarginal = SMarginal - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
algo = "rational"
algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0]
if samples == "pivoted":
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
matchingWeight = 1.
cutOffTolerance = 5.
cutOffType = "POTENTIAL"
if size == 1:
mu0 = [2.7e2, 20]
mutar = [3e2, 25]
murange = [[20., 10], [5.2e2, 30]]
elif size == 2:
mu0 = [2.7e2, 60]
mutar = [3e2, 75]
murange = [[20., 10], [5.2e2, 110]]
elif size == 3:
mu0 = [2.7e2, 160]
mutar = [3e2, 105]
murange = [[20., 10], [5.2e2, 310]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1]]]
n = 100
b = 10
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
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['R'] = SPivot
# params['S'] = 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")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params["S"] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
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:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', 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)))
try:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [1., 1], polesImTol = 2.)
except: pass
if show_norm:
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[1., 1.], relative = False,
nobeta = True)
print(1.j * approx.getPoles([None, 50.]))
diff --git a/examples/2d/pod/membrane_fracture_engine.py b/examples/2d/pod/membrane_fracture_engine.py
new file mode 120000
index 0000000..9e382c7
--- /dev/null
+++ b/examples/2d/pod/membrane_fracture_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/membrane_fracture_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/membrane_fracture_nodomain_engine.py b/examples/2d/pod/membrane_fracture_nodomain_engine.py
new file mode 120000
index 0000000..45ce4c8
--- /dev/null
+++ b/examples/2d/pod/membrane_fracture_nodomain_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/membrane_fracture_nodomain_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/square_pod.py b/examples/2d/pod/square_pod.py
index c089941..1999b22 100644
--- a/examples/2d/pod/square_pod.py
+++ b/examples/2d/pod/square_pod.py
@@ -1,190 +1,190 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_domain_problem_engine import \
HelmholtzSquareDomainProblemEngine as HSDPE
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 = 5
size = 6
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
MN = 6
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
SMarginal = 4
MMarginal = SMarginal - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples == "pivoted":
radialM = ""
radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
cutOffTolerance = 1. * np.inf
cutOffType = "MAGNITUDE"
matchingWeight = 10.
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [6.5 ** .5, .9 ** .5]
murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]]
aEff = 1.#1
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
if ignore_forcing: solver.nbs = 1
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
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"
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RIP
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
if samples in ["standard", "centered"]:
method = RB
elif samples == "pivoted":
raise
# params['S'] = SPivot
# params['SMarginal'] = SMarginal
# params['MMarginal'] = MMarginal
# params['polybasisMarginal'] = "MONOMIAL" + radialM
# params['matchingWeight'] = matchingWeight
# params['radialDirectionalWeightsMarginal'] = radialWeightM
# 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:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), 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":
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, 20, [2., 1.], clip = clip,
relative = True, nobeta = True)
diff --git a/examples/2d/pod/square_pod_hermite.py b/examples/2d/pod/square_pod_hermite.py
index e3d5851..9fc6517 100644
--- a/examples/2d/pod/square_pod_hermite.py
+++ b/examples/2d/pod/square_pod_hermite.py
@@ -1,95 +1,95 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_domain_problem_engine import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
RandomSampler as RS,
ManualSampler as MS)
verb = 0
size = 1
show_sample = False
show_norm = True
MN = 4
R = (MN + 2) * (MN + 1) // 2
MN0 = 2
R0 = (MN0 + 2) * (MN0 + 1) // 2
S0Tensorized = (MN0 + 1) ** 2
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [4.5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # large
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':R, 'POD':True}
params['polybasis'] = "CHEBYSHEV"
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
method = RB
if sampling == "quadrature":
sampler0 = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
S0 = S0Tensorized
else: # if sampling == "random":
sampler0 = RS(murange, "SOBOL", scalingExp = rescalingExp)
S0 = R0
params['sampler'] = MS(murange, points = sampler0.generatePoints(S0))
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), 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":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.])
if show_norm:
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 25, [2., 2.])
diff --git a/examples/2d/pod/square_simplified_pod.py b/examples/2d/pod/square_simplified_pod.py
index 0b3bde9..5dca5e5 100644
--- a/examples/2d/pod/square_simplified_pod.py
+++ b/examples/2d/pod/square_simplified_pod.py
@@ -1,120 +1,120 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_simplified_domain_problem_engine import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = False
show_norm = True
MN = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [8 ** .5, 1 ** .5]
murange = [[6 ** .5, .25 ** .5], [8 ** .5, 1.25 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
params['polybasis'] = "LEGENDRE"
params['polybasis'] = "MONOMIAL"
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
method = RB
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)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), 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":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.])
if show_norm:
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 25, [2., 2.])
diff --git a/examples/2d/pod/square_simplified_pod_hermite.py b/examples/2d/pod/square_simplified_pod_hermite.py
index 69b0bb4..5c6e302 100644
--- a/examples/2d/pod/square_simplified_pod_hermite.py
+++ b/examples/2d/pod/square_simplified_pod_hermite.py
@@ -1,95 +1,95 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
+from helmholtz_square_simplified_domain_problem_engine import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
RandomSampler as RS,
ManualSampler as MS)
verb = 0
size = 1
show_sample = False
show_norm = True
MN = 4
R = (MN + 2) * (MN + 1) // 2
MN0 = 2
R0 = (MN0 + 2) * (MN0 + 1) // 2
S0Tensorized = (MN0 + 1) ** 2
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [4.5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # large
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':R, 'POD':True}
params['polybasis'] = "CHEBYSHEV"
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
method = RB
if sampling == "quadrature":
sampler0 = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
S0 = S0Tensorized
else: # if sampling == "random":
sampler0 = RS(murange, "SOBOL", scalingExp = rescalingExp)
S0 = R0
params['sampler'] = MS(murange, points = sampler0.generatePoints(S0))
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), 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":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.])
if show_norm:
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50, [2., 2.])
diff --git a/examples/2d/pod/synthetic_bivariate_engine.py b/examples/2d/pod/synthetic_bivariate_engine.py
new file mode 120000
index 0000000..989399c
--- /dev/null
+++ b/examples/2d/pod/synthetic_bivariate_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/2d/base/synthetic_bivariate_engine.py
\ No newline at end of file
diff --git a/examples/2d/pod/synthetic_pod.py b/examples/2d/pod/synthetic_pod.py
index 074169d..92d6d3c 100644
--- a/examples/2d/pod/synthetic_pod.py
+++ b/examples/2d/pod/synthetic_pod.py
@@ -1,129 +1,128 @@
import numpy as np
-from rrompy.hfengines.linear_problem.bidimensional import \
- SyntheticBivariateEngine as SBE
+from synthetic_bivariate_engine import SyntheticBivariateEngine as SBE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
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
STensorized = (MN + 1) ** 2
PODTol = 1e-6
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
if samples == "standard":
radial = ""
radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
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)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':R,
'POD':True, 'PODTolerance':PODTol}
method = RB
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)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', 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:
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 2.], clip = clip, relative = False)
diff --git a/examples/3d/base/fracture3.py b/examples/3d/base/fracture3.py
index 13b6a32..7c96293 100644
--- a/examples/3d/base/fracture3.py
+++ b/examples/3d/base/fracture3.py
@@ -1,35 +1,34 @@
-from rrompy.hfengines.linear_problem.tridimensional import \
- MembraneFractureEngine3 as MFE
+from membrane_fracture_engine3 import MembraneFractureEngine3 as MFE
from fracture3_warping import fracture3_warping
verb = 100
mu0 = [45. ** .5, .7, .075]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
#solver.plotmesh(figsize = (7.5, 4.5))
#solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
#solver.plot(uh, what = 'REAL', figsize = (8, 5))
#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], filename = 'out', folder = True,
# forceNewFile = False)
warps = fracture3_warping(solver.V.mesh(), L, mu0[1], delta, mu0[2])
#solver.plotmesh(warps, figsize = (7.5, 4.5))
#solver.plot(u0, warps, what = 'REAL', figsize = (8, 5))
solver.plot(uh, warps, what = 'REAL', figsize = (8, 5))
#solver.plot(solver.residual(mu0, uh)[0], warps, name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], warps, filename = 'outW',
# folder = True, forceNewFile = False)
diff --git a/rrompy/hfengines/linear_problem/tridimensional/membrane_fracture_engine3.py b/examples/3d/base/membrane_fracture_engine3.py
similarity index 99%
rename from rrompy/hfengines/linear_problem/tridimensional/membrane_fracture_engine3.py
rename to examples/3d/base/membrane_fracture_engine3.py
index 6e5421d..f07eb34 100644
--- a/rrompy/hfengines/linear_problem/tridimensional/membrane_fracture_engine3.py
+++ b/examples/3d/base/membrane_fracture_engine3.py
@@ -1,264 +1,262 @@
# 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 fenics as fen
import ufl
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import (nextDerivativeIndices,
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['MembraneFractureEngine3']
-
def gen_mesh_fracture(W, delta, H, n):
n = 2 * (n // 2) + 1
xL = np.linspace(-.5 * W, - .5 * delta, n // 2 + 1)
xR = - xL[::-1]
nEff = int(np.ceil(delta / (xL[1] - xL[0])))
xC = np.linspace(-.5 * delta, .5 * delta, nEff + 1)[1:-1]
yA = np.linspace(-.5 * H, .5 * H, n)
yL = np.linspace(-.5 * H, 0., n // 2 + 1)
editor = fen.MeshEditor()
mesh = fen.Mesh()
editor.open(mesh, "triangle", 2, 2)
editor.init_vertices((2 * n + nEff - 1) * (n // 2 + 1))
editor.init_cells(2 * (2 * n - 2 + nEff) * (n // 2))
for j, y in enumerate(yA):
for i, x in enumerate(xL):
editor.add_vertex(i + j * (n // 2 + 1), np.array([x, y]))
for j in range(n - 1):
for i in range(n // 2):
editor.add_cell(2 * (i + j * (n // 2)),
np.array([i + j * (n // 2 + 1), i + 1 + j * (n // 2 + 1),
i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
editor.add_cell(2 * (i + j * (n // 2)) + 1,
np.array([i + 1 + j * (n // 2 + 1), i + 1 + (j + 1) * (n // 2 + 1),
i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
vBase1, cBase1 = n * (n // 2 + 1), 2 * (n - 1) * (n // 2)
for j, y in enumerate(yA):
for i, x in enumerate(xR):
editor.add_vertex(vBase1 + i + j * (n // 2 + 1), np.array([x, y]))
for j in range(n - 1):
for i in range(n // 2):
editor.add_cell(cBase1 + 2 * (i + j * (n // 2)),
np.array([vBase1 + i + j * (n // 2 + 1), vBase1 + i + 1 + j * (n // 2 + 1),
vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
editor.add_cell(cBase1 + 2 * (i + j * (n // 2)) + 1,
np.array([vBase1 + i + 1 + j * (n // 2 + 1),
vBase1 + i + 1 + (j + 1) * (n // 2 + 1),
vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
vBase2, cBase2 = 2 * n * (n // 2 + 1), 4 * (n - 1) * (n // 2)
for j, y in enumerate(yL):
for i, x in enumerate(xC):
editor.add_vertex(vBase2 + i + j * (nEff - 1), np.array([x, y]))
for j in range(n // 2):
for i in range(nEff - 2):
editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)),
np.array([vBase2 + i + j * (nEff - 1), vBase2 + i + 1 + j * (nEff - 1),
vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)) + 1,
np.array([vBase2 + i + 1 + j * (nEff - 1),
vBase2 + i + 1 + (j + 1) * (nEff - 1),
vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
if nEff == 1:
for j in range(n // 2):
editor.add_cell(cBase2 + 2 * j,
np.array([(j + 1) * (n // 2 + 1) - 1, vBase1 + j * (n // 2 + 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
editor.add_cell(cBase2 + 2 * j + 1,
np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
else:
cBase3 = 2 * (2 * n + nEff - 4) * (n // 2)
for j in range(n // 2):
editor.add_cell(cBase3 + 2 * j,
np.array([(j + 1) * (n // 2 + 1) - 1, vBase2 + j * (nEff - 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
editor.add_cell(cBase3 + 2 * j + 1,
np.array([vBase2 + j * (nEff - 1), vBase2 + (j + 1) * (nEff - 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
cBase4 = 2 * (2 * n + nEff - 3) * (n // 2)
for j in range(n // 2):
editor.add_cell(cBase4 + 2 * j,
np.array([vBase2 + (j + 1) * (nEff - 1) - 1, vBase1 + j * (n // 2 + 1),
vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
editor.add_cell(cBase4 + 2 * j + 1,
np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
editor.close()
return mesh
class MembraneFractureEngine3(HelmholtzProblemEngine):
def __init__(self, mu0 : paramVal = [20. ** .5, .6, .08], H : float = 1.,
L : float = .75, delta : float = .05, n : int = 50,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = False
self.nAs = 6
self.npar = 3
self._H = H
self._delta = delta
self._L = L
self._W = 2 * L + delta
self.rescalingExp = [2., 1., 1.]
mesh = gen_mesh_fracture(self._W, delta, H, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
and x[0] >= - .5 * delta and x[0] <= .5 * delta)
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
self._aboveIndicator = ufl.conditional(ufl.gt(y, 0.), fenONE, fenZERO)
self._belowCenterIndicator = ufl.conditional(
ufl.And(ufl.ge(x, - .5 * delta),
ufl.le(x, .5 * delta)),
fenONE, fenZERO)
self._belowSidesIndicator = (fenONE - self._aboveIndicator
- self._belowCenterIndicator)
self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
- .5 * ((x + .4 * L + .5 * delta) / (.1 * L)) ** 2.
) * self._belowSidesIndicator, fenZERO]
def buildA(self):
"""Build terms of operator of linear system."""
if self.As[1] is None or self.As[2] is None or self.As[5] is None:
thy2 = [(1.,), ('x', '()', 1, '*', 4., '-', 2. * self._H),
('x', '()', 1, '**', 2., '*', 6., '-',
('x', '()', 1, '*', 6. * self._H), '+', self._H ** 2.),
('x', '()', 1, '**', 2., '*', 4., '-',
('x', '()', 1, '*', 3. * self._H), '+', self._H ** 2.,
'*', ('x', '()', 1), '*', 2.),
('x', '()', 1, '*', -1., '+', self._H, '*', ('x', '()', 1),
'**', 2.)]
if self.As[3] is None or self.As[4] is None or self.As[5] is None:
thz2 = [(1.,), ('x', '()', 2, '*', 4., '-', 2. * self._W),
('x', '()', 2, '**', 2., '*', 6., '-',
('x', '()', 2, '*', 6. * self._W), '+', self._W ** 2.),
('x', '()', 2, '**', 2., '*', 4., '-',
('x', '()', 2, '*', 3. * self._W), '+', self._W ** 2.,
'*', ('x', '()', 2), '*', 2.),
('x', '()', 2, '*', -1., '+', self._W, '*', ('x', '()', 2),
'**', 2.)]
if self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fenZERO * fen.dot(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, fenZERO,
self.DirichletBoundary)
ajRe = ((self._aboveIndicator + self._belowSidesIndicator)
* 16. * self._L ** 2. * self.u.dx(0) * self.v.dx(0) * fen.dx)
self.As[1] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thz0 = [(1.,), ('x', '()', 2, '*', 2.), ('x', '()', 2, '**', 2.)]
self.thAs[1] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
for idx in idxs:
dy, dz = 4 - idx[1], 2 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[1] += [thy2[dy] + ('*') + (thz0[dz],)]
else:
self.thAs[1] += [(0.,)]
self.thAs[1] += [None]
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, fenZERO,
self.DirichletBoundary)
ajRe = (4. * self._delta ** 2. * self._belowCenterIndicator
* self.u.dx(0) * self.v.dx(0) * fen.dx)
self.As[2] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thz1 = [(1.,), ('x', '()', 2, '*', 2., '-', 2. * self._W),
('x', '()', 2, '*', -1., '+', self._W, '**', 2.)]
self.thAs[2] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
for idx in idxs:
dy, dz = 4 - idx[1], 2 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[2] += [thy2[dy] + ('*') + (thz1[dz],)]
else:
self.thAs[2] += [(0.,)]
self.thAs[2] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[3] is None:
vbMng(self, "INIT", "Assembling operator term A3.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = (self._H ** 2. * self._aboveIndicator
* self.u.dx(1) * self.v.dx(1) * fen.dx)
self.As[3] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thy1 = [(1.,), ('x', '()', 1, '*', 2., '-', 2. * self._H),
(self._H, '-', ('x', '()', 1), '**', 2.)]
self.thAs[3] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
for idx in idxs:
dy, dz = 2 - idx[1], 4 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[3] += [thy1[dy] + ('*') + (thz2[dz],)]
else:
self.thAs[3] += [(0.,)]
self.thAs[3] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[4] is None:
vbMng(self, "INIT", "Assembling operator term A4.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = ((self._belowSidesIndicator + self._belowCenterIndicator)
* self._H ** 2. * self.u.dx(1) * self.v.dx(1) * fen.dx)
self.As[4] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thy0 = [(1.,), ('x', '()', 1, '*', 2.), ('x', '()', 1, '**', 2.)]
self.thAs[4] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
for idx in idxs:
dy, dz = 2 - idx[1], 4 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[4] += [thy0[dy] + ('*') + (thz2[dz],)]
else:
self.thAs[4] += [(0.,)]
self.thAs[4] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[5] is None:
vbMng(self, "INIT", "Assembling operator term A5.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = - 4. * self.u * self.v * fen.dx
self.As[5] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thx = [(1.,), ('x', '()', 0)]
self.thAs[5] = []
idxs = nextDerivativeIndices([], 3, hashD([1, 4, 4]) + 1)
for idx in idxs:
dx, dy, dz = 1 - idx[0], 4 - idx[1], 4 - idx[2]
if dx >= 0 and dy >= 0 and dz >= 0:
self.thAs[5] += [thx[dx] + ('*') + (thy2[dy],)
+ ('*') + (thz2[dz],)]
else:
self.thAs[5] += [(0.,)]
self.thAs[5] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/examples/3d/base/scattering1.py b/examples/3d/base/scattering1.py
index f9a0c2d..b9805c4 100644
--- a/examples/3d/base/scattering1.py
+++ b/examples/3d/base/scattering1.py
@@ -1,29 +1,29 @@
-from rrompy.hfengines.linear_problem.tridimensional import Scattering1d as S1D
+from scattering_1d import Scattering1d as S1D
import fenics as fen
import numpy as np
verb = 100
mu0 = [4., np.pi, 0.]
mu = [4.5, np.pi, 1.]
n = 100
solver = S1D(mu0 = mu0, n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu)[0]
#solver.plotmesh(figsize = (12, 2))
#solver.plot(u0, what = ['REAL', 'IMAG'], figsize = (12, 4))
print(solver.norm(uh))
#solver.plot(uh, what = ['REAL', 'IMAG'], figsize = (12, 4))
#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = ['REAL', 'IMAG'], figsize = (12, 4))
x = fen.SpatialCoordinate(solver.V.mesh())
warps = [x * mu[1] / solver._L - x, x * solver._L / mu[1] - x]
#solver.plotmesh(warps, figsize = (12, 2))
#solver.plot(u0, warps, what = ['REAL', 'IMAG'], figsize = (12, 4))
solver.plot(uh, warps, what = ['REAL', 'IMAG'], figsize = (12, 4))
#solver.plot(solver.residual(mu0, uh)[0], warps, name = 'res',
# what = ['REAL', 'IMAG'], figsize = (12, 4))
diff --git a/rrompy/hfengines/linear_problem/tridimensional/scattering_1d.py b/examples/3d/base/scattering_1d.py
similarity index 99%
rename from rrompy/hfengines/linear_problem/tridimensional/scattering_1d.py
rename to examples/3d/base/scattering_1d.py
index c88e214..cd7f6ac 100644
--- a/rrompy/hfengines/linear_problem/tridimensional/scattering_1d.py
+++ b/examples/3d/base/scattering_1d.py
@@ -1,73 +1,71 @@
# 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 fenics as fen
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.scattering_problem_engine import (
ScatteringProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['Scattering1d']
-
class Scattering1d(ScatteringProblemEngine):
def __init__(self, mu0 : paramVal = [5.5, np.pi, 0.], n : int = 50,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = False
self.nAs = 3
self.npar = 3
self.rescalingExp = [1., 1., 1.]
self._L = np.real(mu0[1])
self.V = fen.FunctionSpace(fen.IntervalMesh(n, 0., self._L), "P", 1)
self.RobinBoundary = lambda x, on_b: (on_b and x[0] >= .5 * np.pi)
self.DirichletBoundary = "REST"
self.DirichletDatum = fenONE
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, fenZERO,
self.DirichletBoundary)
a0 = fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
self.As[0] = fenics2Sparse(a0, {}, DirichletBC0, 1)
self.thAs[0] = self.getMonomialSingleWeight([0, 0, 0])
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[1] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A1.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a1 = - self._L ** -1. * fen.dot(self.u, self.v) * self.ds(1)
self.As[1] = 1.j * fenics2Sparse(a1, {}, DirichletBC0, 0)
self.thAs[1] = self.getMonomialSingleWeight([1, 1, 1])
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, fenZERO,
self.DirichletBoundary)
a2 = - self._L ** -2. * fen.dot(self.u, self.v) * fen.dx
self.As[2] = fenics2Sparse(a2, {}, DirichletBC0, 0)
self.thAs[2] = self.getMonomialSingleWeight([2, 2, 0])
vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py
index a4e4b74..eadd736 100644
--- a/examples/3d/pod/fracture3_pod.py
+++ b/examples/3d/pod/fracture3_pod.py
@@ -1,241 +1,240 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
-from rrompy.hfengines.linear_problem.tridimensional import \
- MembraneFractureEngine3 as MFE
+from membrane_fracture_engine3 import MembraneFractureEngine3 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 = 70
size = 4
show_sample = False
show_norm = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 5
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
STensorized = (MN + 1) ** 3
PODTol = 1e-8
SPivot = MN + 1
SMarg1d = 5
MMarginal = SMarg1d - 1
SMarginal = (MMarginal + 2) * (MMarginal + 1) // 2
SMarginalTensorized = SMarg1d ** 2
matchingWeight = 10.
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 == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 3
if samples == "pivoted":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [47.5 ** .5, .4, .05]
mutar = [50 ** .5, .45, .07]
murange = [[40 ** .5, .3, .025], [55 ** .5, .5, .075]]
if size == 2:
mu0 = [50 ** .5, .3, .02]
mutar = [55 ** .5, .4, .03]
murange = [[40 ** .5, .1, .005], [60 ** .5, .5, .035]]
if size == 3:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .4, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
if size == 4:
mu0 = [45 ** .5, .4, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .5, .06]]
if size == 5:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
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[0][2] + bEff*murange[1][2]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1., 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
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['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RIP
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':R, 'POD':True, 'PODTolerance':PODTol}
if samples in ["standard", "centered"]:
method = RB
elif samples == "pivoted":
raise
# params['S'] = SPivot
# params['SMarginal'] = SMarginal
# params['MMarginal'] = MMarginal
# params['polybasisMarginal'] = "MONOMIAL" + radialM
# params['matchingWeight'] = matchingWeight
# params['radialDirectionalWeightsMarginal'] = radialWeightM
# 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")
params['SMarginal'] = SMarginalTensorized
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1:], murange[1][1:]], "CHEBYSHEV")
params['polybasisMarginal'] = "CHEBYSHEV" + radialM
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:
from fracture3_warping import fracture3_warping
warps = fracture3_warping(solver.V.mesh(), L, mutar[1], delta, mutar[2])
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', 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)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(np.real(approx.mus(0) ** 2.), np.real(approx.mus(1)),
np.real(approx.mus(2)), '.')
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
from plot_zero_set_3 import plotZeroSet3
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [2., 1., 1.],
# clip = clip)
muZeroScatter = plotZeroSet3(murange, murangeEff, approx, mu0, 50,
[None, None, None], [2., 1., 1.], clip = clip)
if show_norm:
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, .5, .05]) ** 2.)
diff --git a/examples/3d/pod/membrane_fracture_engine3.py b/examples/3d/pod/membrane_fracture_engine3.py
new file mode 120000
index 0000000..e970941
--- /dev/null
+++ b/examples/3d/pod/membrane_fracture_engine3.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/3d/base/membrane_fracture_engine3.py
\ No newline at end of file
diff --git a/examples/3d/pod/scattering1_pod.py b/examples/3d/pod/scattering1_pod.py
index 6f480fa..a4e1da8 100644
--- a/examples/3d/pod/scattering1_pod.py
+++ b/examples/3d/pod/scattering1_pod.py
@@ -1,173 +1,173 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
-from rrompy.hfengines.linear_problem.tridimensional import Scattering1d as S1D
+from scattering_1d import Scattering1d as S1D
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 50
size = 2
show_sample = True
show_norm = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 15
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
STensorized = (MN + 1) ** 3
PODTol = 1e-8
samples = "centered"
samples = "standard"
samples, nDer = "standard_MMM", 15
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 3
assert Delta <= 0
if size == 1:
mu0 = [4., np.pi, 0.]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, -.5], [6., 1.1 * np.pi, .5]]
if size == 2:
mu0 = [4., np.pi, .375]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, 0.], [6., 1.1 * np.pi, .75]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1],
aEff*murange[0][2] + bEff*murange[1][2]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
n = 100
solver = S1D(mu0 = mu0, n = n, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
else: #MMM
params['S'] = nDer * int(np.ceil(R / nDer))
method = RI
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':R, 'POD':True, 'PODTolerance':PODTol}
if samples == "MMM":
params['S'] = nDer * int(np.ceil(R / nDer))
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params['S'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples == "standard_MMM":
if sampling == "quadrature":
SdirEff = int(np.ceil(int(np.ceil(R / nDer)) ** (1. / 3))) ** 3
pts = QS(murange, "CHEBYSHEV").generatePoints(SdirEff)
elif sampling == "quadrature_total":
pts = QST(murange, "CHEBYSHEV").generatePoints(int(np.ceil(R / nDer)))
else: # if sampling == "random":
pts = RS(murange, "HALTON").generatePoints(int(np.ceil(R / nDer)))
params['sampler'] = MS(murange, points = pts)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
import fenics as fen
x = fen.SpatialCoordinate(solver.V.mesh())
warps = [x * mutar[1] / solver._L - x, x * solver._L / mutar[1] - x]
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', 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)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(approx.trainedModel.data.mus(0), approx.trainedModel.data.mus(1),
approx.trainedModel.data.mus(2), '.')
ax.set_xlim3d(murangeEff[0][0], murangeEff[1][0])
ax.set_ylim3d(murangeEff[0][1], murangeEff[1][1])
ax.set_zlim3d(murangeEff[0][2], murangeEff[1][2])
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
if algo == "rational" and approx.N > 0:
from plot_zero_set_3 import plotZeroSet3
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], None], [1., 1., 1.],
clip = clip)
plotZeroSet3(murange, murangeEff, approx, mu0, 100, [None, None, None],
[1., 1., 1.], clip = clip, imagTol = 1e-2)
if show_norm:
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, np.pi, 0.]))
diff --git a/examples/3d/pod/scattering_1d.py b/examples/3d/pod/scattering_1d.py
new file mode 120000
index 0000000..ed22763
--- /dev/null
+++ b/examples/3d/pod/scattering_1d.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/examples/3d/base/scattering_1d.py
\ No newline at end of file
diff --git a/examples/airfoil/airfoil_engine.py b/examples/airfoil/airfoil_engine.py
index f63b57d..a521257 100644
--- a/examples/airfoil/airfoil_engine.py
+++ b/examples/airfoil/airfoil_engine.py
@@ -1,50 +1,50 @@
import numpy as np
import fenics as fen
import ufl
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.solver.fenics import fenONE
PI = np.pi
class AirfoilScatteringEngine(HSP):
def __init__(self, kappa:float, theta:float,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(R = 5, kappa = kappa, theta = theta, n = 1,
degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
mesh = fen.Mesh('../data/mesh/airfoil2412_1.xml')
self.V = fen.FunctionSpace(mesh, "P", 1)
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
self.DirichletDatum = [u0R, u0I]
mu = 1.1
epsilon = .1
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
self.diffusivity = a
diff --git a/examples/airfoil/helmholtz_box_scattering_problem_engine.py b/examples/airfoil/helmholtz_box_scattering_problem_engine.py
new file mode 120000
index 0000000..265edd2
--- /dev/null
+++ b/examples/airfoil/helmholtz_box_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_box_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/base/helmholtz_box_scattering_problem_engine.py b/examples/base/helmholtz_box_scattering_problem_engine.py
new file mode 120000
index 0000000..265edd2
--- /dev/null
+++ b/examples/base/helmholtz_box_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_box_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/base/helmholtz_cavity_scattering_problem_engine.py b/examples/base/helmholtz_cavity_scattering_problem_engine.py
new file mode 120000
index 0000000..de6f867
--- /dev/null
+++ b/examples/base/helmholtz_cavity_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_cavity_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/base/helmholtz_square_bubble_problem_engine.py b/examples/base/helmholtz_square_bubble_problem_engine.py
new file mode 120000
index 0000000..29e8a1c
--- /dev/null
+++ b/examples/base/helmholtz_square_bubble_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_square_bubble_problem_engine.py
\ No newline at end of file
diff --git a/examples/base/helmholtz_square_transmission_problem_engine.py b/examples/base/helmholtz_square_transmission_problem_engine.py
new file mode 120000
index 0000000..61f47ea
--- /dev/null
+++ b/examples/base/helmholtz_square_transmission_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_square_transmission_problem_engine.py
\ No newline at end of file
diff --git a/examples/base/solver.py b/examples/base/solver.py
index 0868f79..b52905c 100644
--- a/examples/base/solver.py
+++ b/examples/base/solver.py
@@ -1,67 +1,67 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_transmission_problem_engine import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_cavity_scattering_problem_engine import \
HelmholtzCavityScatteringProblemEngine as HCSPE
testNo = 2
verb = 0
if testNo == 1:
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20,
verbosity = verb)
mu = 12.**.5
solver.setSolver("BICG", {"tol" : 1e-15})
uh = solver.solve(mu)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
solver.plot(solver.residual(mu, uh)[0], name = 'res')
###########
elif testNo == 2:
solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4.,
n = 50, verbosity = verb)
mu = 4.
uref = solver.liftDirichletData()
uh = solver.solve(mu)[0]
utot = uh - uref
print(solver.norm(uh))
print(solver.norm(uref))
solver.plot(uh)
solver.plot(uref, name = 'u_Dir')
solver.plot(utot, name = 'u_tot')
solver.plot(solver.residual(mu, uh)[0], name = 'res')
solver.plot(solver.residual(mu, utot)[0], name = 'res_tot')
###########
elif testNo == 3:
solver = HBSPE(R = 5, kappa = 12**.5, theta = - np.pi * 60 / 180, n = 30,
verbosity = verb)
mu = 12**.5
uref = solver.liftDirichletData()
uh = solver.solve(mu)[0]
utot = uh - uref
solver.plotmesh()
print(solver.norm(uh))
print(solver.norm(utot))
solver.plot(uh)
solver.plot(utot, name = 'u_tot')
solver.plot(solver.residual(mu, uh)[0], name = 'res')
solver.plot(solver.residual(mu, utot)[0], name = 'res_tot')
###########
elif testNo == 4:
solver = HCSPE(kappa = 5, n = 30, verbosity = verb)
mu = 10
uh = solver.solve(mu)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
solver.plot(solver.residual(mu, uh)[0], name = 'res')
diff --git a/examples/from_papers/greedy_scatteringAirfoil.py b/examples/from_papers/greedy_scatteringAirfoil.py
index 844d4d9..00dafc2 100644
--- a/examples/from_papers/greedy_scatteringAirfoil.py
+++ b/examples/from_papers/greedy_scatteringAirfoil.py
@@ -1,147 +1,147 @@
import numpy as np
import fenics as fen
import ufl
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RI
from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RB
from rrompy.solver.fenics import fenONE, L2NormMatrix
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
verb = 2
timed = False
algo = "RI"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
k0s = np.linspace(5, 20, 25)
k0 = np.mean(k0s)
kl, kr = min(k0s), max(k0s)
params = {'sampler':QS([kl, kr], "UNIFORM"), 'nTestPoints':500,
'greedyTol':1e-2, 'S':10, 'polybasis':polyBasis,
'errorEstimatorKind':'DIAGONAL'}
# 'errorEstimatorKind':'INTERPOLATORY'}
# 'errorEstimatorKind':'AFFINE'}
#########
PI = np.pi
R = 2
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
kappa = 10
theta = PI * - 45 / 180.
mu = 1.1
epsilon = .1
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
solver = HSP(R, np.real(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 2)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
#########
if algo == "RI":
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
params.pop("polybasis")
params.pop("errorEstimatorKind")
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus(0)),
1.05*np.max(norm)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus(0)),
4.*np.max(resApp)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/from_papers/helmholtz_box_scattering_problem_engine.py b/examples/from_papers/helmholtz_box_scattering_problem_engine.py
new file mode 120000
index 0000000..265edd2
--- /dev/null
+++ b/examples/from_papers/helmholtz_box_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_box_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/from_papers/pod_scatteringAirfoil.py b/examples/from_papers/pod_scatteringAirfoil.py
index 2849d33..2e05d7a 100644
--- a/examples/from_papers/pod_scatteringAirfoil.py
+++ b/examples/from_papers/pod_scatteringAirfoil.py
@@ -1,119 +1,119 @@
import numpy as np
import fenics as fen
import ufl
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
from rrompy.solver.fenics import fenONE
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
test = "solve"
test = "approx"
plotSamples = True
k0 = 10
kLeft, kRight = 8 + 0.j, 12 + 0.j
ktar = 11
ktars = np.linspace(8, 12, 21) + 0.j
PI = np.pi
R = 2
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
kappa = 10
theta = PI * - 45 / 180.
mu = 1.1
epsilon = .1
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
###
solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
###
if test == "solve":
uinc = solver.liftDirichletData()
uh = solver.solve(k0)[0]
uhtot = uh - uinc
print(solver.norm(uh))
print(solver.norm(uhtot))
solver.plot(fen.project(a, solver.V).vector(), what = 'Real',
name = 'a')
solver.plot(uinc, what = 'Real', name = 'u_inc')
solver.plot(uh, what = 'ABS')
solver.plot(uhtot, what = 'ABS', name = 'u + u_inc')
elif test == "approx":
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True,
'polybasis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parRI = subdict(params, ['N', 'M', 'S', 'POD', 'polybasis',
'sampler'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxRI = RI(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRI, verbosity = verb)
approxRB = RB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB, verbosity = verb)
approxRI.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxRI.plotSamples()
approxRI.plotHF(ktar, name = 'u_HF')
approxRI.plotApprox(ktar, name = 'u_RI''')
approxRI.plotErr(ktar, name = 'err_RI''')
approxRI.plotRes(ktar, name = 'res_RI''')
approxRB.plotApprox(ktar, name = 'u_RB')
approxRB.plotErr(ktar, name = 'err_RB')
approxRB.plotRes(ktar, name = 'res_RB')
HFNorm, RHSNorm = approxRI.normHF(ktar), approxRI.normRHS(ktar)
RIRes, RIErr = approxRI.normRes(ktar), approxRI.normErr(ktar)
RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar)
print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm))
print('RIRes:\t{}\nRIErr:\t{}'.format(RIRes, RIErr))
print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr))
print('\nPoles RI'':')
print(approxRI.getPoles())
diff --git a/examples/pod/PolesCentered.py b/examples/pod/PolesCentered.py
index 02ca328..9810a5b 100644
--- a/examples/pod/PolesCentered.py
+++ b/examples/pod/PolesCentered.py
@@ -1,70 +1,70 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.utilities.numerical import squareResonances
from rrompy.parameter.parameter_sampling import ManualSampler as MS
verb = 0
k0 = (12+0.j) ** .5
krange = [[9. ** .5], [15. ** .5]]
Nmin, Nmax = 2, 10
Nvals = np.arange(Nmin, Nmax + 1, 2)
params = {'N':Nmin, 'M':0, 'S':Nmin + 1, 'POD':True,
'sampler': MS(krange, points = [k0], scalingExp = 2.)}
#boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x)
# - np.real(z0))
#cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True}
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 25, verbosity = verb)
solver.omega = np.real(k0)
approxP = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approxR = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
rP, rE = [None] * len(Nvals), [None] * len(Nvals)
verbose = 1
for j, N in enumerate(Nvals):
if verbose > 0:
print('N = E = {}'.format(N))
approxP.approxParameters = {'N':N, 'S':N+1}
approxR.approxParameters = {'R':N, 'S':N+1}
if verbose > 1:
print(approxP.approxParameters)
print(approxR.approxParameters)
rP[j] = approxP.getPoles()
rE[j] = approxR.getPoles()
if verbose > 2:
print(rP)
print(rE)
from matplotlib import pyplot as plt
plotRows = int(np.ceil(len(Nvals) / 3))
fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows))
for j, N in enumerate(Nvals):
i1, i2 = int(np.floor(j / 3)), j % 3
axes[i1, i2].set_title('N = E = {}'.format(N))
axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb',
label="Pade'", markersize = 8)
axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), 'Pr',
label="RB", markersize = 8)
axes[i1, i2].axhline(linewidth=1, color='k')
xmin, xmax = axes[i1, i2].get_xlim()
height = (xmax - xmin) / 2.
res = np.power(squareResonances(xmin**2., xmax**2., False), .5)
axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4)
axes[i1, i2].plot(np.real(k0), np.imag(k0), 'om', markersize = 5)
axes[i1, i2].plot(np.real(k0) * np.ones(2),
1.5 * height * np.arange(-1, 3, 2), '--m')
axes[i1, i2].grid()
axes[i1, i2].set_xlim(xmin, xmax)
axes[i1, i2].set_ylim(- height, height)
p = axes[i1, i2].legend()
plt.tight_layout()
for j in range((len(Nvals) - 1) % 3 + 1, 3):
axes[plotRows - 1, j].axis('off')
diff --git a/examples/pod/PolesDistributed.py b/examples/pod/PolesDistributed.py
index 6dd37aa..d4bfc5e 100644
--- a/examples/pod/PolesDistributed.py
+++ b/examples/pod/PolesDistributed.py
@@ -1,45 +1,45 @@
from matplotlib import pyplot as plt
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
from rrompy.utilities.numerical import squareResonances
verb = 0
ks = [1., 46. ** .5]
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, verbosity = verb)
k0 = np.mean(np.power(ks, 2.)) ** .5
k0 = 3.46104724
solver.omega = np.real(k0)
nsets = 12
paramsPade = {'S':2, 'POD':True, 'polybasis':"LEGENDRE",
'sampler':QS(ks, "UNIFORM", rescaleExp = 2.)}
approx = RI(solver, mu0 = k0, approxParameters = paramsPade,
verbosity = verb)
poles = [None] * nsets
samples = [None] * nsets
polesexact = np.unique(np.power(
squareResonances(ks[0]**2., ks[1]**2., False), .5))
for i in range(1, nsets + 1):
print("N = {}".format(4 * i))
approx.approxParameters = {'N': 4 * i, 'M': 4 * i, 'S': 4 * i + 1}
approx.setupApprox()
poles[i - 1] = approx.getPoles()
samples[i - 1] = (approx.mus ** 2.).data
for i in range(1, nsets + 1):
plt.figure()
plt.plot(np.real(poles[i - 1]), np.imag(poles[i - 1]), 'kx')
plt.plot(polesexact, np.zeros_like(polesexact), 'm.')
plt.plot(np.real(samples[i - 1]), np.imag(samples[i - 1]), 'r*')
plt.xlim(ks)
plt.ylim((ks[0] - ks[1]) / 2., (ks[1] - ks[0]) / 2.)
plt.title("N = {}, Neff = {}".format(4 * i, len(poles[i - 1])))
plt.grid()
plt.show()
plt.close()
diff --git a/examples/pod/RationalHermiteInterpolant.py b/examples/pod/RationalHermiteInterpolant.py
index 3d20460..133e32e 100644
--- a/examples/pod/RationalHermiteInterpolant.py
+++ b/examples/pod/RationalHermiteInterpolant.py
@@ -1,131 +1,131 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_transmission_problem_engine import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
testNo = 3
verb = 100
polyBasis = "CHEBYSHEV"
#polyBasis = "LEGENDRE"
#polyBasis = "MONOMIAL"
rep = "REPEAT"
#rep = "TILE"
if testNo == 1:
k0s = np.power([10 + 0.j, 14 + 0.j], .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
samplerBase = QS(k0s, "CHEBYSHEV", 2.)
S = 8
if rep == "REPEAT":
points = np.repeat(samplerBase.generatePoints(2).data,
int(np.ceil(S / 2)))
else: # if rep == "TILE":
points = np.tile(samplerBase.generatePoints(2).data,
int(np.ceil(S / 2)))
params = {'N':7, 'M':6, 'S':S, 'POD':True, 'polybasis':polyBasis,
'sampler':MS(k0s, points = points, 2.)}
ktar = (11 + .5j) ** .5
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 2:
k0s = [3.85 + 0.j, 4.15 + 0.j]
k0 = np.mean(k0s)
ktar = 4 + 0.j
samplerBase = QS(k0s, "CHEBYSHEV", 2.)
S = 10
if rep == "REPEAT":
points = np.repeat(samplerBase.generatePoints(5).data,
int(np.ceil(S / 5)))
else: # if rep == "TILE":
points = np.tile(samplerBase.generatePoints(5).data,
int(np.ceil(S / 5)))
params = {'N':8, 'M':9, 'S':S, 'POD':True, 'polybasis':polyBasis,
'sampler':MS(k0s, points = points, 2.)}
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 3:
k0s = [2., 5.]
k0 = np.mean(k0s)
ktar = 4.5 - 0.j
samplerBase = QS(k0s, "CHEBYSHEV")
S = 15
if rep == "REPEAT":
points = np.repeat(samplerBase.generatePoints(5).data,
int(np.ceil(S / 5)))
else: # if rep == "TILE":
points = np.tile(samplerBase.generatePoints(5).data,
int(np.ceil(S / 5)))
params = {'N':14, 'M':14, 'S':S, 'POD':True, 'polybasis':polyBasis,
'sampler':MS(k0s, points = points)}
solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
diff --git a/examples/pod/RationalInterpolant.py b/examples/pod/RationalInterpolant.py
index f6e4a30..7e214c4 100644
--- a/examples/pod/RationalInterpolant.py
+++ b/examples/pod/RationalInterpolant.py
@@ -1,119 +1,119 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_transmission_problem_engine import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
testNo = 3
verb = 100
polyBasis = "CHEBYSHEV"
polyBasis = "LEGENDRE"
#polyBasis = "MONOMIAL"
loadName = "RationalInterpolantModel.pkl"
if testNo in [1, -1]:
if testNo > 0:
k0s = np.power([10 + 0.j, 14 + 0.j], .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
params = {'N':4, 'M':3, 'S':5, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0s, "CHEBYSHEV", 2.)}
ktar = (11 + .5j) ** .5
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
if testNo > 0:
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
else:
approx = RI(solver, mu0 = 0,
approxParameters = {'S':1, 'muBounds':[0, 1]},
verbosity = verb)
approx.loadTrainedModel(loadName)
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
if testNo > 0:
approx.storeTrainedModel("RationalInterpolantModel",
forceNewFile = False)
print(approx.trainedModel.data.__dict__)
############
elif testNo == 2:
k0s = [3.85 + 0.j, 4.15 + 0.j]
k0 = np.mean(k0s)
ktar = 4 + 0.j
params = {'N':8, 'M':9, 'S':10, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0s, "CHEBYSHEV", 2.)}
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 3:
k0s = [2., 5.]
k0 = np.mean(k0s)
ktar = 4.5 - 0.j
params = {'N':10, 'M':9, 'S':15, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0s, "CHEBYSHEV")}
solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RI(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
print('\nPoles Pade'':')
print(approx.getPoles())
diff --git a/examples/pod/ReducedBasis.py b/examples/pod/ReducedBasis.py
index ab42249..14a82b0 100644
--- a/examples/pod/ReducedBasis.py
+++ b/examples/pod/ReducedBasis.py
@@ -1,106 +1,106 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_transmission_problem_engine import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.linear_problem import \
+from helmholtz_box_scattering_problem_engine import \
HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
testNo = 1
verb = 100
loadName = "ReducedBasisModel.pkl"
if testNo in [1, -1]:
if testNo > 0:
k0s = np.power([10 + 0.j, 14 + 0.j], .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
params = {'S':5, 'R':4, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", 2.)}
ktar = (11 + .5j) ** .5
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
if testNo > 0:
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
else:
approx = RB(solver, mu0 = 0,
approxParameters = {'S':1, 'muBounds':[0, 1]},
verbosity = verb)
approx.loadTrainedModel(loadName)
approx.plotApprox(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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 testNo > 0:
approx.storeTrainedModel("ReducedBasisModel", forceNewFile = False)
print(approx.trainedModel.data.__dict__)
############
elif testNo == 2:
k0s = [3.85 + 0.j, 4.15 + 0.j]
k0 = np.mean(k0s)
ktar = 4 + .15j
params = {'S':10, 'R':9, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", 2.)}
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
############
elif testNo == 3:
k0s = [2., 5.]
k0 = np.mean(k0s)
ktar = 4.5 - 0.j
params = {'S':15, 'R':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")}
solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApprox(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
diff --git a/examples/pod/helmholtz_box_scattering_problem_engine.py b/examples/pod/helmholtz_box_scattering_problem_engine.py
new file mode 120000
index 0000000..265edd2
--- /dev/null
+++ b/examples/pod/helmholtz_box_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_box_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/pod/helmholtz_cavity_scattering_problem_engine.py b/examples/pod/helmholtz_cavity_scattering_problem_engine.py
new file mode 120000
index 0000000..de6f867
--- /dev/null
+++ b/examples/pod/helmholtz_cavity_scattering_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_cavity_scattering_problem_engine.py
\ No newline at end of file
diff --git a/examples/pod/helmholtz_square_bubble_problem_engine.py b/examples/pod/helmholtz_square_bubble_problem_engine.py
new file mode 120000
index 0000000..29e8a1c
--- /dev/null
+++ b/examples/pod/helmholtz_square_bubble_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_square_bubble_problem_engine.py
\ No newline at end of file
diff --git a/examples/pod/helmholtz_square_transmission_problem_engine.py b/examples/pod/helmholtz_square_transmission_problem_engine.py
new file mode 120000
index 0000000..61f47ea
--- /dev/null
+++ b/examples/pod/helmholtz_square_transmission_problem_engine.py
@@ -0,0 +1 @@
+/home/pradovera/Desktop/Repos/RROMPy/tests/hfengines/helmholtz_square_transmission_problem_engine.py
\ No newline at end of file
diff --git a/examples/pod/scatteringSquare.py b/examples/pod/scatteringSquare.py
index 163f778..a5b49c9 100644
--- a/examples/pod/scatteringSquare.py
+++ b/examples/pod/scatteringSquare.py
@@ -1,72 +1,72 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_cavity_scattering_problem_engine import \
HelmholtzCavityScatteringProblemEngine as CSPE
from rrompy.reduction_methods.standard import RationalInterpolant as PD
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
####################
test = "solve"
#test = "approx"
plotSamples = True
k0 = 10
kLeft, kRight = 9, 11
ktar = 9.5
ktars = np.linspace(8.5, 11.5, 125)
#ktars = np.array([k0])
kappa = 5
n = 50
solver = CSPE(kappa = kappa, n = n, verbosity = verb)
solver.omega = k0
if test == "solve":
uh = solver.solve(k0)[0]
print(solver.norm(uh))
solver.plot(uh, what = ['ABS', 'REAL'])
elif test == "approx":
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True,
'polybasis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'polybasis',
'sampler'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = PD(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parPade,
verbosity = verb)
approxRB = RB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB, verbosity = verb)
approxPade.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxPade.plotSamples()
PadeErr, solNorm = approxPade.normErr(ktar), approxPade.normHF(ktar)
RBErr = approxRB.normErr(ktar)
print(('SolNorm:\t{}\nErrPade:\t{}\nErrRelPade:\t{}\nErrRB:\t\t{}'
'\nErrRelRB:\t{}').format(solNorm, PadeErr,
np.divide(PadeErr, solNorm), RBErr,
np.divide(RBErr, solNorm)))
print('\nPoles Pade'':')
print(approxPade.getPoles())
print('\nPoles RB:')
print(approxRB.getPoles())
approxPade.plotHF(ktar, name = 'u_ex')
approxPade.plotApprox(ktar, name = 'u_Pade''')
approxRB.plotApprox(ktar, name = 'u_RB')
approxPade.plotErr(ktar, name = 'errPade''')
approxRB.plotErr(ktar, name = 'errRB')
diff --git a/examples/pod/with_error_plot.py b/examples/pod/with_error_plot.py
index c436d68..f95cc61 100644
--- a/examples/pod/with_error_plot.py
+++ b/examples/pod/with_error_plot.py
@@ -1,111 +1,111 @@
import numpy as np
-from rrompy.hfengines.linear_problem import \
+from helmholtz_square_bubble_problem_engine import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import RationalMovingLeastSquares as RIM
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
verb = 100
sol = "single"
sol = "sweep"
algo = "RI"
algo = "RIM"
#algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
radialBasis = "GAUSSIAN"
#radialBasis = "THINPLATE"
#radialBasis = "MULTIQUADRIC"
#radialBasis = "NEARESTNEIGHBOR"
radialBasisDen = "GAUSSIAN"
radialBasisDen = "THINPLATE"
radialBasisDen = "MULTIQUADRIC"
radialBasisDen = "NEARESTNEIGHBOR"
k0sC = np.power([7 + 0.j, 55 + 0.j], .5)
k0 = np.mean(k0sC ** 2.) ** .5
ktar = 14 ** .5
n = 20
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
params = {'N':1, 'M':1, 'S':30, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0sC, "CHEBYSHEV", 2.)}
if algo == "RI":
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
elif algo == "RIM":
params["radialBasis"] = radialBasis
params["radialDirectionalWeights"] = [.75 * params["S"]]
params["radialBasisDen"] = radialBasisDen
params["nNearestNeighborDen"] = params["N"] + 1
approx = RIM(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
params.pop("N")
params.pop("M")
params.pop("polybasis")
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.setupApprox()
if sol == "single":
approx.plotSamples(what = "REAL")
approx.plotApprox(ktar, what = "REAL", name = "uApp")
approx.plotHF(ktar, what = "REAL", name = "uHF")
approx.plotErr(ktar, what = "REAL", name = "err")
approx.plotRes(ktar, what = "REAL", name = "res")
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
poles = approx.getPoles()
try:
print('Poles:', poles ** 2.)
except: pass
if sol == "sweep":
z0s = np.real(np.linspace(k0sC[0] ** 2., k0sC[1] ** 2., 100))
k0s = z0s ** .5
zl, zr = min(z0s), max(z0s)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
# normRHS = approx.normRHS(k0s)
norm = approx.normHF(k0s)
normApp = approx.normApprox(k0s)
# res = approx.normRes(k0s) / normRHS
# err = approx.normErr(k0s) / norm
plt.figure()
plt.semilogy(z0s, norm)
plt.semilogy(z0s, normApp, '--')
plt.semilogy(np.real(approx.mus.data) ** 2.,
1.05*np.max(norm)*np.ones_like(approx.mus.data, dtype = float),
'rx')
plt.xlim([zl, zr])
plt.grid()
plt.show()
plt.close()
# plt.figure()
# plt.semilogy(z0s, res)
# plt.xlim([zl, zr])
# plt.grid()
# plt.show()
# plt.close()
#
# plt.figure()
# plt.semilogy(z0s, err)
# plt.xlim([zl, zr])
# plt.grid()
# plt.show()
# plt.close()
#for j, k in enumerate(k0s):
# print(k ** 2., approx.getPoles(mu = k) ** 2., norm[j], normApp[j])
diff --git a/rrompy/hfengines/linear_problem/__init__.py b/rrompy/hfengines/linear_problem/__init__.py
index 3f7fb0b..6dab090 100644
--- a/rrompy/hfengines/linear_problem/__init__.py
+++ b/rrompy/hfengines/linear_problem/__init__.py
@@ -1,49 +1,30 @@
# 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 .
#
from .laplace_base_problem_engine import LaplaceBaseProblemEngine
from .helmholtz_problem_engine import HelmholtzProblemEngine
from .scattering_problem_engine import ScatteringProblemEngine
-from .helmholtz_box_scattering_problem_engine import \
- HelmholtzBoxScatteringProblemEngine
-from .helmholtz_cavity_scattering_problem_engine import \
- HelmholtzCavityScatteringProblemEngine
-from .helmholtz_square_bubble_problem_engine import \
- HelmholtzSquareBubbleProblemEngine
-from .helmholtz_square_bubble_domain_problem_engine import \
- HelmholtzSquareBubbleDomainProblemEngine
-from .helmholtz_square_transmission_problem_engine import \
- HelmholtzSquareTransmissionProblemEngine
-from .laplace_disk_gaussian import LaplaceDiskGaussian
-from .membrane_fracture_engine_nodomain import MembraneFractureEngineNoDomain
__all__ = [
'LaplaceBaseProblemEngine',
'HelmholtzProblemEngine',
- 'ScatteringProblemEngine',
- 'HelmholtzBoxScatteringProblemEngine',
- 'HelmholtzCavityScatteringProblemEngine',
- 'HelmholtzSquareBubbleProblemEngine',
- 'HelmholtzSquareBubbleDomainProblemEngine',
- 'HelmholtzSquareTransmissionProblemEngine',
- 'LaplaceDiskGaussian',
- 'MembraneFractureEngineNoDomain'
+ 'ScatteringProblemEngine'
]
diff --git a/rrompy/hfengines/linear_problem/bidimensional/__init__.py b/rrompy/hfengines/linear_problem/bidimensional/__init__.py
deleted file mode 100644
index 54ca3ff..0000000
--- a/rrompy/hfengines/linear_problem/bidimensional/__init__.py
+++ /dev/null
@@ -1,38 +0,0 @@
-# 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 .
-#
-
-from .cookie_engine_single import CookieEngineSingle
-from .membrane_fracture_engine import MembraneFractureEngine
-from .laplace_disk_gaussian_2 import LaplaceDiskGaussian2
-from .helmholtz_square_simplified_domain_problem_engine import \
- HelmholtzSquareSimplifiedDomainProblemEngine
-from .helmholtz_square_domain_problem_engine import \
- HelmholtzSquareDomainProblemEngine
-from .synthetic_bivariate_engine import SyntheticBivariateEngine
-
-__all__ = [
- 'CookieEngineSingle',
- 'MembraneFractureEngine',
- 'LaplaceDiskGaussian2',
- 'HelmholtzSquareSimplifiedDomainProblemEngine',
- 'HelmholtzSquareDomainProblemEngine',
- 'SyntheticBivariateEngine'
- ]
-
-
-
diff --git a/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py b/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
deleted file mode 100644
index 5e0b8b7..0000000
--- a/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
+++ /dev/null
@@ -1,117 +0,0 @@
-# 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 fenics as fen
-import ufl
-from rrompy.utilities.base.types import ScOp, List, paramVal
-from rrompy.solver.fenics import fenONE, fenZERO
-from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
- HelmholtzProblemEngine)
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.numerical import hashDerivativeToIdx as hashD
-from rrompy.solver.fenics import fenics2Sparse
-
-__all__ = ['CookieEngineSingle']
-
-class CookieEngineSingle(HelmholtzProblemEngine):
-
- def __init__(self, kappa:float, theta:float, n:int, R : int = 1.,
- L : int = 2., nX : int = 1, nY : int = 1,
- mu0 : paramVal = [12. ** .5, 1.],
- degree_threshold : int = np.inf, verbosity : int = 10,
- timestamp : bool = True):
- super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
- verbosity = verbosity, timestamp = timestamp)
- self._affinePoly = False
- self.nAs = 3
- self.npar = 2
- self.rescalingExp = [2., 2.]
- mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(L * nX, L * nY),
- n * nX, n * nY)
- self.V = fen.FunctionSpace(mesh, "P", 1)
- x, y = fen.SpatialCoordinate(mesh)[:]
- cxs = np.linspace(0, L * nX, 2 * nX + 1)[1::2]
- cys = np.linspace(0, L * nY, 2 * nY + 1)[1::2]
- self.cookieIn = fenZERO
- for cx in cxs:
- for cy in cys:
- self.cookieIn += ufl.conditional(
- ufl.le((x-cx)**2. + (y-cy)**2., R**2.),
- fenONE, fenZERO)
- self.cookieOut = fenONE - self.cookieIn
- c, s = np.cos(theta), np.sin(theta)
- self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
- fen.sin(kappa * (c * x + s * y))]
-
- def buildA(self):
- """Build terms of operator of linear system."""
- if self.As[0] is None:
- self.autoSetDS()
- vbMng(self, "INIT", "Assembling operator term A0.", 20)
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- hRe, hIm = self.RobinDatumH
- termNames = ["diffusivity", "RobinDatumH"]
- parsRe = self.iterReduceQuadratureDegree(zip(
- [aRe, hRe],
- [x + "Real" for x in termNames]))
- parsIm = self.iterReduceQuadratureDegree(zip(
- [aIm, hIm],
- [x + "Imag" for x in termNames]))
- a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- + hRe * fen.dot(self.u, self.v) * self.ds(1))
- a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- + hIm * fen.dot(self.u, self.v) * self.ds(1))
- self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
- + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
- self.thAs[0] = self.getMonomialSingleWeight([0, 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, fenZERO,
- self.DirichletBoundary)
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
- ["refractionIndexSquaredReal"]))
- parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
- ["refractionIndexSquaredImag"]))
- a1Re = - n2Re * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
- a1Im = - n2Im * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
- self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
- + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
- self.thAs[1] = self.getMonomialSingleWeight([1, 0])
- 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, fenZERO,
- self.DirichletBoundary)
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
- ["refractionIndexSquaredReal"]))
- parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
- ["refractionIndexSquaredImag"]))
- a2Re = - n2Re * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
- a2Im = - n2Im * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
- self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
- + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
- self.thAs[2] = self.getMonomialSingleWeight([1, 1])
- vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
deleted file mode 100644
index 332d3d7..0000000
--- a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
+++ /dev/null
@@ -1,86 +0,0 @@
-# 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 numpy.polynomial.polynomial import polyfit as fit
-import fenics as fen
-from rrompy.utilities.base.types import paramVal
-from .laplace_base_problem_engine import LaplaceBaseProblemEngine
-from rrompy.solver.fenics import fenZERO
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.solver.fenics import fenics2Vector
-
-__all__ = ['LaplaceDiskGaussian']
-
-class LaplaceDiskGaussian(LaplaceBaseProblemEngine):
- """
- Solver for disk Laplace problems with parametric forcing term center.
- - \Delta u = C exp(-.5 * ||\cdot - (mu, 0)||^2) in \Omega = B(0, 5)
- u = 0 on \partial\Omega.
- """
-
- def __init__(self, n:int, mu0 : paramVal = [0.],
- degree_threshold : int = np.inf, verbosity : int = 10,
- timestamp : bool = True):
- super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
- verbosity = verbosity, timestamp = timestamp)
- self.nbs = 19
- import mshr
- mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), 3 * n)
- self.V = fen.FunctionSpace(mesh, "P", 1)
-
- def buildb(self):
- """Build terms of operator of linear system."""
- if self.thbs[0] is None:
- self.thbs = self.getMonomialWeights(self.nbs)
- bDEIMCoeffs = None
- for j in range(self.nbs):
- if self.bs[j] is None:
- vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
- 20)
- if bDEIMCoeffs is None:
- bDEIM = np.empty((self.nbs, self.spacedim),
- dtype = np.complex)
- muDEIM = 3. * np.linspace(0., 1., self.nbs // 2 + 1) ** 2.
- muDEIM = np.concatenate((-muDEIM[:0:-1], muDEIM))
- for jj, muD in enumerate(muDEIM):
- x, y = fen.SpatialCoordinate(self.V.mesh())[:]
- C = np.exp(-.5 * muD ** 2.)
- CR, CI = np.real(C), np.imag(C)
- f0 = ((2 * np.pi) ** -.5
- * fen.exp(-.5 * (x ** 2. + y ** 2.)))
- muR, muI = np.real(muD), np.imag(muD)
- f1R = fen.exp(muR * x) * fen.cos(muI * x)
- f1I = fen.exp(muR * x) * fen.sin(muI * x)
- fRe = f0 * (CR * f1R - CI * f1I) + fenZERO
- fIm = f0 * (CR * f1I + CI * f1R) + fenZERO
- parsRe = self.iterReduceQuadratureDegree(zip([fRe],
- ["forcingTerm{}Real".format(jj)]))
- parsIm = self.iterReduceQuadratureDegree(zip([fIm],
- ["forcingTerm{}Imag".format(jj)]))
- LR = fen.dot(fRe, self.v) * fen.dx
- LI = fen.dot(fIm, self.v) * fen.dx
- DBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1)
- + 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
- bDEIMCoeffs = (fit(muDEIM / 3., bDEIM, self.nbs - 1).T
- * np.power(3., - np.arange(self.nbs))).T
- self.bs[j] = bDEIMCoeffs[j]
- vbMng(self, "DEL", "Done assembling forcing term.", 20)
-
diff --git a/rrompy/hfengines/linear_problem/tridimensional/__init__.py b/rrompy/hfengines/linear_problem/tridimensional/__init__.py
deleted file mode 100644
index f0493e4..0000000
--- a/rrompy/hfengines/linear_problem/tridimensional/__init__.py
+++ /dev/null
@@ -1,30 +0,0 @@
-# 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 .
-#
-
-from .membrane_fracture_engine3 import MembraneFractureEngine3
-from .scattering_1d import Scattering1d
-from .matrix_dynamical_passive import MatrixDynamicalPassive
-
-__all__ = [
- 'MembraneFractureEngine3',
- 'Scattering1d',
- 'MatrixDynamicalPassive'
- ]
-
-
-
diff --git a/rrompy/hfengines/vector_linear_problem/__init__.py b/rrompy/hfengines/vector_linear_problem/__init__.py
index 805fc8b..689207b 100644
--- a/rrompy/hfengines/vector_linear_problem/__init__.py
+++ b/rrompy/hfengines/vector_linear_problem/__init__.py
@@ -1,32 +1,30 @@
# 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 .
#
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from .linear_elasticity_helmholtz_problem_engine import LinearElasticityHelmholtzProblemEngine
from .linear_elasticity_helmholtz_problem_engine_damped import LinearElasticityHelmholtzProblemEngineDamped
-from .linear_elasticity_beam_poisson_ratio import LinearElasticityBeamPoissonRatio
__all__ = [
'LinearElasticityProblemEngine',
'LinearElasticityHelmholtzProblemEngine',
- 'LinearElasticityHelmholtzProblemEngineDamped',
- 'LinearElasticityBeamPoissonRatio'
+ 'LinearElasticityHelmholtzProblemEngineDamped'
]
diff --git a/rrompy/hfengines/vector_linear_problem/bidimensional/__init__.py b/rrompy/hfengines/vector_linear_problem/bidimensional/__init__.py
deleted file mode 100644
index 3311e14..0000000
--- a/rrompy/hfengines/vector_linear_problem/bidimensional/__init__.py
+++ /dev/null
@@ -1,26 +0,0 @@
-# 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 .
-#
-
-from .linear_elasticity_beam_elasticity_constants import LinearElasticityBeamElasticityConstants
-
-__all__ = [
- 'LinearElasticityBeamElasticityConstants'
- ]
-
-
-
diff --git a/rrompy/hfengines/vector_linear_problem/bidimensional/linear_elasticity_beam_elasticity_constants.py b/rrompy/hfengines/vector_linear_problem/bidimensional/linear_elasticity_beam_elasticity_constants.py
deleted file mode 100644
index 5ca7e55..0000000
--- a/rrompy/hfengines/vector_linear_problem/bidimensional/linear_elasticity_beam_elasticity_constants.py
+++ /dev/null
@@ -1,119 +0,0 @@
-# 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 fenics as fen
-from rrompy.hfengines.vector_linear_problem.\
- linear_elasticity_beam_poisson_ratio import LinearElasticityBeamPoissonRatio
-from rrompy.solver.fenics import fenZERO, fenZEROS
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
-
-__all__ = ['LinearElasticityBeamElasticityConstants']
-
-class LinearElasticityBeamElasticityConstants(
- LinearElasticityBeamPoissonRatio):
- """
- Solver for linear elasticity problem of a beam subject to its own weight,
- with parametric Joung modulus and 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, E0:float, nu0:float,
- length:float, degree_threshold : int = np.inf,
- verbosity : int = 10, timestamp : bool = True):
- super().__init__(mu0 = [nu0, E0], degree_threshold = degree_threshold,
- verbosity = verbosity, timestamp = timestamp)
- self._affinePoly = False
- self.nAs, self.nbs = 3, 2
- self.npar = 2
-
- 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)), 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, 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 = fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
- self.As[1] = fenics2Sparse(a1Re, {}, DirichletBC0, 0)
- self.thAs[1] = [('x', '()', 0, '*', -2., '+', 1.,
- '*', ('x', '()', 1)),
- ('x', '()', 1, '*', -2.),
- ('x', '()', 0, '*', -2.),
- (0.,), (-2.,), None]
- 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)
- a4Re = fen.div(self.u) * fen.div(self.v) * fen.dx
- self.As[2] = fenics2Sparse(a4Re, {}, DirichletBC0, 0)
- self.thAs[2] = self.getMonomialSingleWeight([1, 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)
- L0 = 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(L0, {}, DBCR, 1)
- + 1.j * fenics2Vector(L0, {}, DBCI, 1))
- self.thbs[0] = self.getMonomialSingleWeight([0, 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.),
- (0.,), (-2.0,), None]
- vbMng(self, "DEL", "Done assembling forcing term.", 20)
-
diff --git a/tests/hfengines/fracture.py b/tests/hfengines/fracture.py
index 576f130..2a96b1f 100644
--- a/tests/hfengines/fracture.py
+++ b/tests/hfengines/fracture.py
@@ -1,58 +1,57 @@
# 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 \
+from membrane_fracture_engine import MembraneFractureEngine as MFE
+from membrane_fracture_nodomain_engine import MembraneFractureNoDomainEngine \
as MFEND
from rrompy.solver.fenics import affine_warping
@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(mu0[0], uh1), dual = True)[0],
solver2.norm(solver2.residual(mu0, uh2), dual = True)[0],
atol = 1e-5)
assert np.isclose(solver1.norm(uh1 - uh2), 0., atol = 1e-6)
diff --git a/rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py b/tests/hfengines/helmholtz_box_scattering_problem_engine.py
similarity index 95%
rename from rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py
rename to tests/hfengines/helmholtz_box_scattering_problem_engine.py
index add93e4..7b57010 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py
+++ b/tests/hfengines/helmholtz_box_scattering_problem_engine.py
@@ -1,57 +1,55 @@
# 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 fenics as fen
-from .scattering_problem_engine import ScatteringProblemEngine
-
-__all__ = ['HelmholtzBoxScatteringProblemEngine']
+from rrompy.hfengines.linear_problem import ScatteringProblemEngine
class HelmholtzBoxScatteringProblemEngine(ScatteringProblemEngine):
"""
Solver for scattering problem outside a box with parametric wavenumber.
- \Delta u - omega^2 * n^2 * u = 0 in \Omega
u = 0 on \Gamma_D
\partial_nu - i k u = 0 on \Gamma_R
with exact solution a transmitted plane wave.
"""
def __init__(self, R:float, kappa:float, theta:float, n:int,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = [kappa], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
import mshr
scatterer = mshr.Polygon([fen.Point(-1, -.5), fen.Point(1, -.5),
fen.Point(1, .5), fen.Point(.8, .5),
fen.Point(.8, -.3), fen.Point(-.8, -.3),
fen.Point(-.8, .5), fen.Point(-1, .5),])
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer,
3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.DirichletBoundary = (lambda x, on_boundary:
on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R)
self.RobinBoundary = "REST"
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
self.DirichletDatum = [u0R, u0I]
diff --git a/rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py b/tests/hfengines/helmholtz_cavity_scattering_problem_engine.py
similarity index 95%
rename from rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py
rename to tests/hfengines/helmholtz_cavity_scattering_problem_engine.py
index b59cf5b..62e288a 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py
+++ b/tests/hfengines/helmholtz_cavity_scattering_problem_engine.py
@@ -1,57 +1,55 @@
# 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 fenics as fen
-from .scattering_problem_engine import ScatteringProblemEngine
-
-__all__ = ['HelmholtzCavityScatteringProblemEngine']
+from rrompy.hfengines.linear_problem import ScatteringProblemEngine
class HelmholtzCavityScatteringProblemEngine(ScatteringProblemEngine):
"""
Solver for scattering problem inside a cavity with parametric wavenumber.
- \Delta u - omega^2 * n^2 * u = 0 in \Omega
u = 0 on \Gamma_D
\partial_nu - i k u = 0 on \Gamma_R
with exact solution a transmitted plane wave.
"""
def __init__(self, kappa:float, n:int, gamma : float = 0.,
signR : int = -1, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(mu0 = [kappa], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.signR = signR
pi = np.pi
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.RobinBoundary = (lambda x, on_boundary: on_boundary
and np.isclose(x[1], np.pi))
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
C = 4. / pi ** 4.
bR = C * (2 * (x * (pi - x) + y * (2 * pi - y))
+ (kappa * gamma) ** 2. * x * (pi - x) * y * (2 * pi - y))
bI = C * signR * 2 * kappa * (gamma * (pi - 2 * x) * y * (pi - y)
+ 2 * x * (pi - x) * (pi - y))
wR = fen.cos(kappa * signR * (gamma * x + y))
wI = fen.sin(kappa * signR * (gamma * x + y))
self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI]
diff --git a/tests/hfengines/helmholtz_external.py b/tests/hfengines/helmholtz_external.py
index 13f6ba6..9a4493a 100644
--- a/tests/hfengines/helmholtz_external.py
+++ b/tests/hfengines/helmholtz_external.py
@@ -1,68 +1,70 @@
# 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)
+from helmholtz_box_scattering_problem_engine import \
+ HelmholtzBoxScatteringProblemEngine
+from helmholtz_cavity_scattering_problem_engine import \
+ HelmholtzCavityScatteringProblemEngine
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(mu, uh)[0], dual = True),
4.25056407e-13, 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.setthAs(solver1.thAs)
solver2.setbs(solver1.bs)
solver2.setthbs(solver1.thbs)
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(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(mu, uh)[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 ea30bf2..0a13eae 100644
--- a/tests/hfengines/helmholtz_internal.py
+++ b/tests/hfengines/helmholtz_internal.py
@@ -1,93 +1,96 @@
# 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 os, shutil
import numpy as np
-from rrompy.hfengines.linear_problem import (
- HelmholtzSquareBubbleDomainProblemEngine, HelmholtzSquareBubbleProblemEngine,
- HelmholtzSquareTransmissionProblemEngine)
+from helmholtz_square_bubble_domain_problem_engine import \
+ HelmholtzSquareBubbleDomainProblemEngine
+from helmholtz_square_bubble_problem_engine import \
+ HelmholtzSquareBubbleProblemEngine
+from helmholtz_square_transmission_problem_engine import \
+ HelmholtzSquareTransmissionProblemEngine
def test_helmholtz_square_io():
solver = HelmholtzSquareBubbleProblemEngine(kappa = 4, theta = 1., n = 20,
verbosity = 0)
mu = 5
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 145.0115, rtol = 1e-3)
assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
1.19934e-11, 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 = 20, verbosity = 0)
mu = 5.
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 138.6609, rtol = 1e-2)
assert np.isclose(solver.norm(solver.residual(mu, uh)[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")
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), 10.07843, rtol = 1e-2)
assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
6.14454989e-13, rtol = 1e-1)
diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py b/tests/hfengines/helmholtz_square_bubble_domain_problem_engine.py
similarity index 95%
rename from rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
rename to tests/hfengines/helmholtz_square_bubble_domain_problem_engine.py
index e3b5098..0ac331d 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
+++ b/tests/hfengines/helmholtz_square_bubble_domain_problem_engine.py
@@ -1,131 +1,131 @@
# 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 numpy.polynomial.polynomial import polyfit as fit
import fenics as fen
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO
-from .helmholtz_problem_engine import HelmholtzProblemEngine
+from rrompy.hfengines.linear_problem import HelmholtzProblemEngine
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
-__all__ = ['HelmholtzSquareBubbleDomainProblemEngine']
-
class HelmholtzSquareBubbleDomainProblemEngine(HelmholtzProblemEngine):
"""
Solver for square bubble Helmholtz problems with parametric domain heigth.
- \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi]
u = 0 on \Gamma_mu = \partial\Omega_mu
with exact solution square bubble times plane wave.
"""
def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [1.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = False
self.nAs, self.nbs = 2, 15
self.kappa = kappa
self.theta = theta
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(np.pi, np.pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.rescalingExp = [1.]
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, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * 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 A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2)
k2n2Re = k2Re * n2Re - k2Im * n2Im
k2n2Im = k2Re * n2Im + k2Im * n2Re
parsRe = self.iterReduceQuadratureDegree(zip([k2n2Re],
["kappaSquaredRefractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([k2n2Im],
["kappaSquaredRefractionIndexSquaredImag"]))
a2Re = (fen.dot(self.u.dx(0), self.v.dx(0))
- k2n2Re * fen.dot(self.u, self.v)) * fen.dx
a2Im = - k2n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
self.thAs[1] = self.getMonomialSingleWeight([2])
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 = self.getMonomialWeights(self.nbs)
bDEIMCoeffs = None
for j in range(self.nbs):
if self.bs[j] is None:
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
if bDEIMCoeffs is None:
- bDEIM = np.empty((self.nbs, self.spacedim),
- dtype = np.complex)
muDEIM = np.linspace(.5, 4., self.nbs)
for jj, muD in enumerate(muDEIM):
pi = np.pi
c, s = np.cos(self.theta), np.sin(self.theta)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
muR, muI = np.real(muD), np.imag(muD)
mu2R, mu2I = np.real(muD ** 2.), np.imag(muD ** 2.)
C = 16. / pi ** 4.
bR = C * (2 * (x * (pi - x) + y * (pi - y))
+ (self.kappa * s) ** 2. * (mu2R - 1.)
* x * (pi - x) * y * (pi - y))
bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y
* (pi - y) + s * x * (pi - x)
* (pi - 2 * y))
+ (self.kappa * s) ** 2. * mu2I
* x * (pi - x) * y * (pi - y))
wR = (fen.cos(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
wI = (fen.sin(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
fRe, fIm = bR * wR + bI * wI, bI * wR - bR * wI
fRe = mu2R * fRe - mu2I * fIm + fenZERO
fIm = mu2R * fIm + mu2I * fRe + fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTerm{}Real".format(jj)]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTerm{}Imag".format(jj)]))
LR = fen.dot(fRe, self.v) * fen.dx
LI = fen.dot(fIm, self.v) * fen.dx
DBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1)
+ bjj = (fenics2Vector(LR, parsRe, DBC0, 1)
+ 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
+ if jj == 0:
+ bDEIM = np.empty((self.nbs, len(bjj)),
+ dtype = np.complex)
+ bDEIM[jj] = bjj
bDEIMCoeffs = fit(muDEIM, bDEIM, self.nbs - 1)
self.bs[j] = bDEIMCoeffs[j]
vbMng(self, "DEL", "Done assembling forcing term.", 20)
diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py b/tests/hfengines/helmholtz_square_bubble_problem_engine.py
similarity index 94%
rename from rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py
rename to tests/hfengines/helmholtz_square_bubble_problem_engine.py
index 81675d2..f0a6fdd 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py
+++ b/tests/hfengines/helmholtz_square_bubble_problem_engine.py
@@ -1,51 +1,49 @@
# 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 fenics as fen
-from .helmholtz_problem_engine import HelmholtzProblemEngine
-
-__all__ = ['HelmholtzSquareBubbleProblemEngine']
+from rrompy.hfengines.linear_problem import HelmholtzProblemEngine
class HelmholtzSquareBubbleProblemEngine(HelmholtzProblemEngine):
"""
Solver for square bubble Helmholtz problems with parametric wavenumber.
- \Delta u - omega^2 * u = f in \Omega
u = 0 on \Gamma_D
with exact solution square bubble times plane wave.
"""
def __init__(self, kappa:float, theta:float, n:int,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = [kappa], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
pi = np.pi
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
C = 16. / pi ** 4.
bR = C * 2 * (x * (pi - x) + y * (pi - y))
bI = C * 2 * kappa * (c * (pi - 2 * x) * y * (pi - y)
+ s * x * (pi - x) * (pi - 2 * y))
wR = fen.cos(kappa * (c * x + s * y))
wI = fen.sin(kappa * (c * x + s * y))
self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI]
diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py b/tests/hfengines/helmholtz_square_transmission_problem_engine.py
similarity index 96%
rename from rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py
rename to tests/hfengines/helmholtz_square_transmission_problem_engine.py
index ea45b45..5460ef0 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py
+++ b/tests/hfengines/helmholtz_square_transmission_problem_engine.py
@@ -1,72 +1,70 @@
# 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 fenics as fen
import ufl
-from .helmholtz_problem_engine import HelmholtzProblemEngine
-
-__all__ = ['HelmholtzSquareTransmissionProblemEngine']
+from rrompy.hfengines.linear_problem import HelmholtzProblemEngine
class HelmholtzSquareTransmissionProblemEngine(HelmholtzProblemEngine):
"""
Solver for square transmission Helmholtz problems with parametric
wavenumber.
- \Delta u - omega^2 * n^2 * u = 0 in \Omega
u = 0 on \Gamma_D
with exact solution a transmitted plane wave.
"""
def __init__(self, nT:float, nB:float, kappa:float, theta:float, n:int,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = [kappa], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
mesh = fen.RectangleMesh(fen.Point(-np.pi/2, -np.pi/2),
fen.Point(np.pi/2, np.pi/2), 3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
dx, dy = np.cos(theta), np.sin(theta)
Kx = kappa * nB * dx
Ky = kappa * (nT**2. - (nB * dx)**2. + 0.j)**.5
T = 2 * kappa * nB * dy / (Ky + kappa * nB * dy)
x, y = fen.SpatialCoordinate(mesh)[:]
TR, TI = np.real(T), np.imag(T)
if np.isclose(np.imag(Ky), 0.):
u0RT = (TR * fen.cos(Kx * x + np.real(Ky) * y)
- TI * fen.sin(Kx * x + np.real(Ky) * y))
u0IT = (TR * fen.sin(Kx * x + np.real(Ky) * y)
+ TI * fen.cos(Kx * x + np.real(Ky) * y))
else:
u0RT = fen.exp(- np.imag(Ky) * y) * (TR * fen.cos(Kx * x)
- TI * fen.sin(Kx * x))
u0IT = fen.exp(- np.imag(Ky) * y) * (TR * fen.sin(Kx * x)
+ TI * fen.cos(Kx * x))
u0RB = (fen.cos(kappa * nB * (dx * x + dy * y))
+ (TR - 1) * fen.cos(kappa * nB * (dx*x - dy*y))
- TI * fen.sin(kappa * nB * (dx*x - dy*y)))
u0IB = (fen.sin(kappa * nB * (dx * x + dy * y))
+ (TR - 1) * fen.sin(kappa * nB * (dx*x - dy*y))
+ TI * fen.cos(kappa * nB * (dx*x - dy*y)))
u0R = ufl.conditional(ufl.ge(y, 0.), u0RT, u0RB)
u0I = ufl.conditional(ufl.ge(y, 0.), u0IT, u0IB)
self.refractionIndex = ufl.conditional(ufl.ge(y, 0.),
fen.Constant(nT),
fen.Constant(nB))
self.DirichletDatum = [u0R, u0I]
diff --git a/tests/hfengines/laplace.py b/tests/hfengines/laplace.py
index f13bd4f..f6f9fef 100644
--- a/tests/hfengines/laplace.py
+++ b/tests/hfengines/laplace.py
@@ -1,39 +1,38 @@
# 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
+from laplace_disk_gaussian import LaplaceDiskGaussian, 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(mu, uh)[0], dual = True),
4.66728e-10, atol = 1e-7)
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(mu, uh)[0], dual = True),
2.40043363e-08, atol = 1e-7)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py b/tests/hfengines/laplace_disk_gaussian.py
similarity index 56%
rename from rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py
rename to tests/hfengines/laplace_disk_gaussian.py
index 4b240bc..bb39579 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py
+++ b/tests/hfengines/laplace_disk_gaussian.py
@@ -1,97 +1,156 @@
# 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 numpy.polynomial.polynomial import polyfit as fit
import fenics as fen
from rrompy.utilities.base.types import paramVal
-from rrompy.hfengines.linear_problem.laplace_disk_gaussian import (
- LaplaceDiskGaussian)
+from rrompy.hfengines.linear_problem import LaplaceBaseProblemEngine
from rrompy.solver.fenics import fenZERO
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
from rrompy.solver.fenics import fenics2Vector
-__all__ = ['LaplaceDiskGaussian2']
+class LaplaceDiskGaussian(LaplaceBaseProblemEngine):
+ """
+ Solver for disk Laplace problems with parametric forcing term center.
+ - \Delta u = C exp(-.5 * ||\cdot - (mu, 0)||^2) in \Omega = B(0, 5)
+ u = 0 on \partial\Omega.
+ """
+
+ def __init__(self, n:int, mu0 : paramVal = [0.],
+ degree_threshold : int = np.inf, verbosity : int = 10,
+ timestamp : bool = True):
+ super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
+ verbosity = verbosity, timestamp = timestamp)
+ self.nbs = 19
+ import mshr
+ mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), 3 * n)
+ self.V = fen.FunctionSpace(mesh, "P", 1)
+
+ def buildb(self):
+ """Build terms of operator of linear system."""
+ if self.thbs[0] is None:
+ self.thbs = self.getMonomialWeights(self.nbs)
+ bDEIMCoeffs = None
+ for j in range(self.nbs):
+ if self.bs[j] is None:
+ vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
+ 20)
+ if bDEIMCoeffs is None:
+ muDEIM = 3. * np.linspace(0., 1., self.nbs // 2 + 1) ** 2.
+ muDEIM = np.concatenate((-muDEIM[:0:-1], muDEIM))
+ for jj, muD in enumerate(muDEIM):
+ x, y = fen.SpatialCoordinate(self.V.mesh())[:]
+ C = np.exp(-.5 * muD ** 2.)
+ CR, CI = np.real(C), np.imag(C)
+ f0 = ((2 * np.pi) ** -.5
+ * fen.exp(-.5 * (x ** 2. + y ** 2.)))
+ muR, muI = np.real(muD), np.imag(muD)
+ f1R = fen.exp(muR * x) * fen.cos(muI * x)
+ f1I = fen.exp(muR * x) * fen.sin(muI * x)
+ fRe = f0 * (CR * f1R - CI * f1I) + fenZERO
+ fIm = f0 * (CR * f1I + CI * f1R) + fenZERO
+ parsRe = self.iterReduceQuadratureDegree(zip([fRe],
+ ["forcingTerm{}Real".format(jj)]))
+ parsIm = self.iterReduceQuadratureDegree(zip([fIm],
+ ["forcingTerm{}Imag".format(jj)]))
+ LR = fen.dot(fRe, self.v) * fen.dx
+ LI = fen.dot(fIm, self.v) * fen.dx
+ DBC0 = fen.DirichletBC(self.V, fenZERO,
+ self.DirichletBoundary)
+ bjj = (fenics2Vector(LR, parsRe, DBC0, 1)
+ + 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
+ if jj == 0:
+ bDEIM = np.empty((self.nbs, len(bjj)),
+ dtype = np.complex)
+ bDEIM[jj] = bjj
+ bDEIMCoeffs = (fit(muDEIM / 3., bDEIM, self.nbs - 1).T
+ * np.power(3., - np.arange(self.nbs))).T
+ self.bs[j] = bDEIMCoeffs[j]
+ vbMng(self, "DEL", "Done assembling forcing term.", 20)
class LaplaceDiskGaussian2(LaplaceDiskGaussian):
"""
Solver for disk Laplace problems with parametric forcing term center.
- \Delta u = C exp(-.5 * ||\cdot - (mu1, mu2)||^2) in \Omega = B(0, 5)
u = 0 on \partial\Omega.
"""
def __init__(self, n:int, mu0 : paramVal = [0., 0.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(n = n, mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.nbs = 16
self.npar = 2
def buildb(self):
"""Build terms of operator of linear system."""
if self.thbs[0] is None:
self.thbs = [None] * self.nbs
bDEIMCoeffs = None
for j in range(self.nbs):
j1, j2 = j % int(self.nbs ** .5), j // int(self.nbs ** .5)
if self.bs[j] is None:
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
if bDEIMCoeffs is None:
- bDEIM = np.empty((self.nbs, self.spacedim),
- dtype = np.complex)
muD1 = np.linspace(-2., 2., int(self.nbs ** .5))
muDEIM = np.empty((self.nbs, 2))
muDEIM[:, 0] = np.repeat(muD1, int(self.nbs ** .5))
muDEIM[:, 1] = np.tile(muD1, int(self.nbs ** .5))
for jj, muD in enumerate(muDEIM):
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
C = np.exp(-.5 * (muD[0] ** 2. + muD[1] ** 2.))
CR, CI = np.real(C), np.imag(C)
f0 = ((2 * np.pi) ** -.5
* fen.exp(-.5 * (x ** 2. + y ** 2.)))
muxR, muxI = np.real(muD[0]), np.imag(muD[0])
muyR, muyI = np.real(muD[1]), np.imag(muD[1])
f1R = (fen.exp(muxR * x + muyR * y)
* fen.cos(muxI * x + muyI * y))
f1I = (fen.exp(muxR * x + muyR * y)
* fen.sin(muxI * x + muyI * y))
fRe = f0 * (CR * f1R - CI * f1I) + fenZERO
fIm = f0 * (CR * f1I + CI * f1R) + fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTerm{}Real".format(jj)]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTerm{}Imag".format(jj)]))
LR = fen.dot(fRe, self.v) * fen.dx
LI = fen.dot(fIm, self.v) * fen.dx
DBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1)
+ bjj = (fenics2Vector(LR, parsRe, DBC0, 1)
+ 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
+ if jj == 0:
+ bDEIM = np.empty((self.nbs, len(bjj)),
+ dtype = np.complex)
+ bDEIM[jj] = bjj
p = PI()
wellCond, _ = p.setupByInterpolation(muDEIM, bDEIM,
int(self.nbs ** .5) - 1,
"MONOMIAL", False, False)
bDEIMCoeffs = p.coeffs
self.bs[j1 + int(self.nbs ** .5) * j2] = bDEIMCoeffs[j1, j2]
self.thbs[j1 + int(self.nbs ** .5) * j2] = (
self.getMonomialSingleWeight([j1, j2]))
vbMng(self, "DEL", "Done assembling forcing term.", 20)
diff --git a/tests/hfengines/linear_elasticity.py b/tests/hfengines/linear_elasticity.py
index c79bddb..85578ef 100644
--- a/tests/hfengines/linear_elasticity.py
+++ b/tests/hfengines/linear_elasticity.py
@@ -1,38 +1,38 @@
# 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 linear_elasticity_beam_poisson_ratio_engine import \
+ LinearElasticityBeamPoissonRatioEngine
from rod_3d import rod3Dsolver
def test_elastic_beam():
- solver = LinearElasticityBeamPoissonRatio(n = 10, rho_ = 1e3, g = 3,
+ solver = LinearElasticityBeamPoissonRatioEngine(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), 9.33957e-08, rtol = 1e-1)
assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
8.4545952e-13, 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], dual = True),
1.2210129e-07, rtol = 1e-1)
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py b/tests/hfengines/linear_elasticity_beam_poisson_ratio_engine.py
similarity index 96%
rename from rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
rename to tests/hfengines/linear_elasticity_beam_poisson_ratio_engine.py
index 624ab2e..1a6ed59 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
+++ b/tests/hfengines/linear_elasticity_beam_poisson_ratio_engine.py
@@ -1,116 +1,115 @@
# 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 fenics as fen
-from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
+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
-__all__ = ['LinearElasticityBeamPoissonRatio']
-
-class LinearElasticityBeamPoissonRatio(LinearElasticityProblemEngine):
+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, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(mu0 = [nu0], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
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)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py b/tests/hfengines/membrane_fracture_engine.py
similarity index 97%
rename from rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py
rename to tests/hfengines/membrane_fracture_engine.py
index b3ab7d6..59e03f3 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py
+++ b/tests/hfengines/membrane_fracture_engine.py
@@ -1,145 +1,142 @@
# 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 fenics as fen
import mshr, ufl
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO, fenONE
-from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
- HelmholtzProblemEngine)
+from rrompy.hfengines.linear_problem import HelmholtzProblemEngine
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['MembraneFractureEngine']
-
class MembraneFractureEngine(HelmholtzProblemEngine):
def __init__(self, mu0 : paramVal = [20. ** .5, .6], H : float = 1.,
L : float = .75, delta : float = .05, n : int = 50,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = False
self.nAs = 5
self.npar = 2
self.H = H
self.rescalingExp = [2., 1.]
domain = (mshr.Rectangle(fen.Point(0., - H / 2.),
fen.Point(2. * L + delta, H / 2.))
- mshr.Rectangle(fen.Point(L, 0.),
fen.Point(L + delta, H / 2.)))
mesh = mshr.generate_mesh(domain, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
and x[0] >= L
and x[0] <= L + delta)
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
self._belowIndicator = ufl.conditional(ufl.le(y, 0.), fenONE, fenZERO)
self._aboveIndicator = fenONE - self._belowIndicator
self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
- .5 * ((x - .6 * L) / (.1 * L)) ** 2.
) * self._belowIndicator, fenZERO]
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, fenZERO,
self.DirichletBoundary)
a0Re = fenZERO * fen.dot(self.u, self.v) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
self.thAs[0] = self.getMonomialSingleWeight([0, 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, fenZERO,
self.DirichletBoundary)
a1Re = (self.H ** 3 / 4. * self._aboveIndicator
* fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx)
self.As[1] = fenics2Sparse(a1Re, {}, DirichletBC0, 0)
self.thAs[1] = [('x', '()', 1, '*', -2., '+', self.H),
(0.,), (-2.,), None]
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, fenZERO,
self.DirichletBoundary)
a2Re = self.H ** 2 * fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
self.thAs[2] = [('x', '()', 1, '*', -1., '+', self.H, '*',
('x', '()', 1), '**', 2.),
(0.,), ('x', '()', 1, '**', 2., '*', 4., '-',
('x', '()', 1, '*', 6. * self.H), '+',
2. * self.H ** 2., '*', ('x', '()', 1)),
(0.,), (0.,), ('x', '()', 1, '**', 2., '*', 6., '-',
('x', '()', 1, '*', 6. * self.H), '+',
self.H ** 2.),
(0.,), (0.,), (0.,), ('x', '()', 1, '*', 4.,
'-', 2. * self.H),
(0.,), (0.,), (0.,), (0.,), (1.,), None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[3] is None:
vbMng(self, "INIT", "Assembling operator term A3.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a3Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a3Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[3] = (fenics2Sparse(a3Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a3Im, parsIm, DirichletBC0, 0))
self.thAs[3] = [('x', '()', 1, '*', -1., '+', self.H, '*',
('x', '()', 1), '**', 2., '*', ('x', '()', 0)),
('x', '()', 1, '*', -1., '+', self.H, '*',
('x', '()', 1), '**', 2.),
(2. * self.H ** 2., '-',
('x', '()', 1, '*', 6. * self.H), '+',
('x', '()', 1, '**', 2., '*', 4.), '*',
('x', '()', 1), '*', ('x', '()', 0)),
(0.,), ('x', '()', 1, '**', 2., '*', 4., '-',
('x', '()', 1, '*', 6. * self.H), '+',
2. * self.H ** 2., '*', ('x', '()', 1)),
('x', '()', 1, '**', 2., '*', 6., '-',
('x', '()', 1, '*', 6. * self.H), '+',
self.H ** 2., '*', ('x', '()', 0)), (0.,), (0.,),
('x', '()', 1, '**', 2., '*', 6., '-',
('x', '()', 1, '*', 6. * self.H), '+',
self.H ** 2.),
('x', '()', 1, '*', 4., '-', 2. * self.H, '*',
('x', '()', 0)), (0.,), (0.,), (0.,),
('x', '()', 1, '*', 4., '-', 2. * self.H),
('x', '()', 0), (0.,), (0.,), (0.,), (0.,),
(1.,), None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[4] is None:
vbMng(self, "INIT", "Assembling operator term A4.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a4Re = .25 * self.H ** 2 * fen.dot(self.u.dx(1),
self.v.dx(1)) * fen.dx
self.As[4] = fenics2Sparse(a4Re, {}, DirichletBC0, 0)
self.thAs[4] = self.getMonomialSingleWeight([0, 2])
vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py b/tests/hfengines/membrane_fracture_nodomain_engine.py
similarity index 97%
rename from rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py
rename to tests/hfengines/membrane_fracture_nodomain_engine.py
index 5c17226..c471691 100644
--- a/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py
+++ b/tests/hfengines/membrane_fracture_nodomain_engine.py
@@ -1,91 +1,89 @@
# 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 fenics as fen
import mshr, ufl
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse
-__all__ = ['MembraneFractureEngineNoDomain']
-
-class MembraneFractureEngineNoDomain(HelmholtzProblemEngine):
+class MembraneFractureNoDomainEngine(HelmholtzProblemEngine):
def __init__(self, mu0 : paramVal = [20. ** .5, .6], H : float = 1.,
L : float = .75, delta : float = .05, n : int = 50,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0[0], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.npar = 1
self.lFrac = mu0[1]
self.H = H
self.rescalingExp = [2.]
domain = (mshr.Rectangle(fen.Point(0., - H / 2.),
fen.Point(2. * L + delta, H / 2.))
- mshr.Rectangle(fen.Point(L, 0.),
fen.Point(L + delta, H / 2.)))
mesh = mshr.generate_mesh(domain, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
and x[0] >= L
and x[0] <= L + delta)
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
self._belowIndicator = ufl.conditional(ufl.le(y, 0.), fenONE, fenZERO)
self._aboveIndicator = fenONE - self._belowIndicator
self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
- .5 * ((x - .6 * L) / (.1 * L)) ** 2.
) * self._belowIndicator, fenZERO]
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = (fen.dot(self.u.dx(0), self.v.dx(0))
+ self.H ** 4 / 4. * (self.lFrac ** -2. * self._aboveIndicator
+ (self.H - self.lFrac) ** -2. * self._belowIndicator)
* fen.dot(self.u.dx(1), self.v.dx(1))
) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
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, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)