diff --git a/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py b/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
index 5afb153..49e4acd 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
@@ -1,167 +1,137 @@
# 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 scsp
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 verbosityDepth
from rrompy.utilities.poly_fitting.polynomial 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.nAs = 5
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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
self.As[2] = self.checkAInBounds(-1)
if derI <= 3 and self.As[3] is None:
self.As[3] = self.checkAInBounds(-1)
if derI <= 4 and self.As[4] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A4.",
timestamp = self.timestamp)
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"]))
a4Re = - n2Re * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
a4Im = - n2Im * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
- A4Re = fen.assemble(a4Re, form_compiler_parameters = parsRe)
- A4Im = fen.assemble(a4Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A4Re)
- DirichletBC0.zero(A4Im)
- A4ReMat = fen.as_backend_type(A4Re).mat()
- A4ImMat = fen.as_backend_type(A4Im).mat()
- A4Rer, A4Rec, A4Rev = A4ReMat.getValuesCSR()
- A4Imr, A4Imc, A4Imv = A4ImMat.getValuesCSR()
- self.As[4] = (scsp.csr_matrix((A4Rev, A4Rec, A4Rer),
- shape = A4ReMat.size)
- + 1.j * scsp.csr_matrix((A4Imv, A4Imc, A4Imr),
- shape = A4ImMat.size))
+ self.As[4] = (fenics2Sparse(a4Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a4Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py
index eafce7b..cf05c29 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py
@@ -1,228 +1,203 @@
# 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 scsp
import fenics as fen
from rrompy.utilities.base.types import (ScOp, List, Tuple, paramVal, Np1D,
FenExpr)
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
+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.nAs, self.nbs = 6, 11 * 12 // 2
self.npar = 2
self.rescalingExp = [2., 1.]
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 / self.mu0(0, 1) * y)),
fen.sin(kappa * (c * x + s / self.mu0(0, 1) * y))]
self.forcingTermDer = kappa * s * y
def getExtraFactorB(self, mu : paramVal = [],
derI : int = 0) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
mu = self.checkParameter(mu)
def getPowMinusj(x, power):
powR = x ** power
powI = fenZERO
if power % 2 == 1:
powR, powI = powI, powR
if power % 4 > 1:
powR, powI = - powR, - powI
return powR, powI
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
if derI == 0: return fenONE, fenZERO
coeffs = np.zeros(derI + 1)
coeffs[1] = - 2.
for j in range(2, derI + 1):
coeffs[1 :] = (-2. / j * coeffs[1 :]
- (3 - (1 + 2 * np.arange(derI)) / j) * coeffs[: -1])
for j in range(derI):
powR, powI = getPowMinusj(self.forcingTermDer, derI - j)
mupBase = coeffs[j + 1] * mu(0, 1) ** (- 3 * derI + 2 * j)
mupR, mupI = np.real(mupBase), np.imag(mupBase)
if j == 0:
exprR = mupR * powR - mupI * powI
exprI = mupI * powR + mupR * powI
else:
exprR += mupR * powR - mupI * powI
exprI += mupI * powR + mupR * powI
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return exprR, exprI
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
for j in range(2, 5):
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 5 and self.As[5] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A5.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a5Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- A5Re = fen.assemble(a5Re)
- DirichletBC0.zero(A5Re)
- A5ReMat = fen.as_backend_type(A5Re).mat()
- A5Rer, A5Rec, A5Rev = A5ReMat.getValuesCSR()
- self.As[5] = scsp.csr_matrix((A5Rev, A5Rec, A5Rer),
- shape = A5ReMat.size,
- dtype = np.complex)
+ self.As[5] = fenics2Sparse(a5Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
derH = hashI(j, self.npar)
if bs[j] is None:
if np.sum(derH) != derH[-1]:
if homogeneized:
self.bsH[j] = self.checkbInBounds(-1)
else:
self.bs[j] = self.checkbInBounds(-1)
continue
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j == 0:
u0Re, u0Im = self.DirichletDatum
else:
u0Re, u0Im = fenZERO, fenZERO
if j < self.nbs:
fRe, fIm = self.forcingTerm
cRe, cIm = self.getExtraFactorB(self.mu0, derH[-1])
cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(j)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(j)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
- DBCR.apply(b0Re)
- DBCI.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if homogeneized:
Ader = self.A(self.mu0, derH)
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py
index d310906..68a0f00 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_simplified_domain_problem_engine.py
@@ -1,129 +1,106 @@
# 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 scsp
import fenics as fen
from rrompy.utilities.base.types import ScOp, List, paramVal
from rrompy.solver.fenics import fenZERO
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+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.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)[:]
- C = 16. / pi ** 4.
self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
fen.sin(kappa * (c * x + s * y))]
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a2Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- A2Re = fen.assemble(a2Re)
- DirichletBC0.zero(A2Re)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- self.As[2] = scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size,
- dtype = np.complex)
+ self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py b/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py
index ef85fa7..cd3ec26 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py
@@ -1,125 +1,123 @@
# 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 Np1D, Tuple, List, FenExpr, paramVal
from rrompy.hfengines.linear_problem.laplace_disk_gaussian import (
LaplaceDiskGaussian)
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.solver.fenics import fenics2Vector
__all__ = ['LaplaceDiskGaussian2']
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 = 1
self.npar = 2
def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
mu = self.checkParameter(mu)
if mu != self.forcingTermMu:
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling base expression for "
"forcing term."),
timestamp = self.timestamp)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
C = np.exp(-.5 * (mu(0, 0) ** 2. + mu(0, 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(mu(0, 0)), np.imag(mu(0, 0))
muyR, muyI = np.real(mu(0, 1)), np.imag(mu(0, 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)
self.forcingTerm = [f0 * (CR * f1R - CI * f1I),
f0 * (CR * f1I + CI * f1R)]
self.forcingTermMu = mu
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling base expression.",
timestamp = self.timestamp)
return self.forcingTerm
def computebsFactors(self):
pass
def getExtraFactorB(self, mu : paramVal = [],
derI : int = 0) -> Tuple[FenExpr, FenExpr]:
if derI == 0: return [fenONE, fenZERO]
raise RROMPyException("Not implemented.")
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
if bs[j] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j < self.nbs:
fRe, fIm = self.getForcingTerm(mu)
cRe, cIm = self.getExtraFactorB(mu, j)
cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(j)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(j)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- DirichletBC0.apply(b0Re)
- DirichletBC0.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1))
if homogeneized:
Ader = self.A(self.mu0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py b/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py
index 6e03ad0..98c1119 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/membrane_fracture_engine.py
@@ -1,248 +1,188 @@
# 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 scsp
import fenics as fen
import mshr, ufl
from rrompy.utilities.base.types import ScOp, List, paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+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.nAs = 20
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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
for j in [1, 3, 4, 6, 7, 10, 11, 12, 15, 16, 17, 18]:
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = (self.H ** 4 / 4. * self._aboveIndicator
* fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx)
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a2Re = (- self.H ** 3 / 2. * self._aboveIndicator
* fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx)
- A2Re = fen.assemble(a2Re)
- DirichletBC0.zero(A2Re)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- self.As[2] = scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size,
- dtype = np.complex)
+ self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 5 and self.As[5] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A6.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a5Re = self.H ** 2 * (fen.dot(self.u.dx(0), self.v.dx(0))
+ .25 * fen.dot(self.u.dx(1), self.v.dx(1))) * fen.dx
- A5Re = fen.assemble(a5Re)
- DirichletBC0.zero(A5Re)
- A5ReMat = fen.as_backend_type(A5Re).mat()
- A5Rer, A5Rec, A5Rev = A5ReMat.getValuesCSR()
- self.As[5] = scsp.csr_matrix((A5Rev, A5Rec, A5Rer),
- shape = A5ReMat.size,
- dtype = np.complex)
+ self.As[5] = fenics2Sparse(a5Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 8 and self.As[8] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A8.",
timestamp = self.timestamp)
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"]))
a8Re = - self.H ** 2. * n2Re * fen.dot(self.u, self.v) * fen.dx
a8Im = - self.H ** 2. * n2Im * fen.dot(self.u, self.v) * fen.dx
- A8Re = fen.assemble(a8Re, form_compiler_parameters = parsRe)
- A8Im = fen.assemble(a8Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A8Re)
- DirichletBC0.zero(A8Im)
- A8ReMat = fen.as_backend_type(A8Re).mat()
- A8ImMat = fen.as_backend_type(A8Im).mat()
- A8Rer, A8Rec, A8Rev = A8ReMat.getValuesCSR()
- A8Imr, A8Imc, A8Imv = A8ImMat.getValuesCSR()
- self.As[8] = (scsp.csr_matrix((A8Rev, A8Rec, A8Rer),
- shape = A8ReMat.size)
- + 1.j * scsp.csr_matrix((A8Imv, A8Imc, A8Imr),
- shape = A8ImMat.size))
+ self.As[8] = (fenics2Sparse(a8Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a8Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 9 and self.As[9] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A9.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a9Re = - 2. * self.H * fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
- A9Re = fen.assemble(a9Re)
- DirichletBC0.zero(A9Re)
- A9ReMat = fen.as_backend_type(A9Re).mat()
- A9Rer, A9Rec, A9Rev = A9ReMat.getValuesCSR()
- self.As[9] = scsp.csr_matrix((A9Rev, A9Rec, A9Rer),
- shape = A9ReMat.size,
- dtype = np.complex)
+ self.As[9] = fenics2Sparse(a9Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 13 and self.As[13] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A13.",
timestamp = self.timestamp)
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"]))
a13Re = 2. * self.H * n2Re * fen.dot(self.u, self.v) * fen.dx
a13Im = 2. * self.H * n2Im * fen.dot(self.u, self.v) * fen.dx
- A13Re = fen.assemble(a13Re, form_compiler_parameters = parsRe)
- A13Im = fen.assemble(a13Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A13Re)
- DirichletBC0.zero(A13Im)
- A13ReMat = fen.as_backend_type(A13Re).mat()
- A13ImMat = fen.as_backend_type(A13Im).mat()
- A13Rer, A13Rec, A13Rev = A13ReMat.getValuesCSR()
- A13Imr, A13Imc, A13Imv = A13ImMat.getValuesCSR()
- self.As[13] = (scsp.csr_matrix((A13Rev, A13Rec, A13Rer),
- shape = A13ReMat.size)
- + 1.j * scsp.csr_matrix((A13Imv, A13Imc, A13Imr),
- shape = A13ImMat.size))
+ self.As[13] = (fenics2Sparse(a13Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a13Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 14 and self.As[14] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A14.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a14Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
- A14Re = fen.assemble(a14Re)
- DirichletBC0.zero(A14Re)
- A14ReMat = fen.as_backend_type(A14Re).mat()
- A14Rer, A14Rec, A14Rev = A14ReMat.getValuesCSR()
- self.As[14] = scsp.csr_matrix((A14Rev, A14Rec, A14Rer),
- shape = A14ReMat.size,
- dtype = np.complex)
+ self.As[14] = fenics2Sparse(a14Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 19 and self.As[19] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A19.",
timestamp = self.timestamp)
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"]))
a19Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a19Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
- A19Re = fen.assemble(a19Re, form_compiler_parameters = parsRe)
- A19Im = fen.assemble(a19Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A19Re)
- DirichletBC0.zero(A19Im)
- A19ReMat = fen.as_backend_type(A19Re).mat()
- A19ImMat = fen.as_backend_type(A19Im).mat()
- A19Rer, A19Rec, A19Rev = A19ReMat.getValuesCSR()
- A19Imr, A19Imc, A19Imv = A19ImMat.getValuesCSR()
- self.As[19] = (scsp.csr_matrix((A19Rev, A19Rec, A19Rer),
- shape = A19ReMat.size)
- + 1.j * scsp.csr_matrix((A19Imv, A19Imc, A19Imr),
- shape = A19ImMat.size))
+ self.As[19] = (fenics2Sparse(a19Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a19Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py b/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
index 21f2b87..2fb7ca5 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
@@ -1,153 +1,123 @@
# 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 scsp
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 verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+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.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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
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
- A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
- A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A2Re)
- DirichletBC0.zero(A2Im)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2ImMat = fen.as_backend_type(A2Im).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
- self.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size)
- + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
- shape = A2ImMat.size))
+ self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
index a77dd46..5b04852 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
@@ -1,160 +1,140 @@
# 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 scsp
import fenics as fen
from .laplace_base_problem_engine import LaplaceBaseProblemEngine
from rrompy.utilities.base.types import List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+from rrompy.solver.fenics import fenics2Sparse
__all__ = ['HelmholtzProblemEngine']
class HelmholtzProblemEngine(LaplaceBaseProblemEngine):
"""
Solver for generic Helmholtz problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
refractionIndex: Value of n.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, 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.nAs = 2
self.rescalingExp = [2.]
self.refractionIndex = fenONE
@property
def refractionIndex(self):
"""Value of n."""
return self._refractionIndex
@refractionIndex.setter
def refractionIndex(self, refractionIndex):
self.resetAs()
if not isinstance(refractionIndex, (list, tuple,)):
refractionIndex = [refractionIndex, fenZERO]
self._refractionIndex = refractionIndex
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.As[0] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
index 6babb7c..9e27a74 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
@@ -1,246 +1,227 @@
# 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 scsp
import fenics as fen
from rrompy.utilities.base.types import (Np1D, ScOp, Tuple, List, FenExpr,
paramVal)
from rrompy.solver.fenics import fenZERO
from .helmholtz_problem_engine import HelmholtzProblemEngine
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
+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.nAs, self.nbs = 3, 20
self.kappa = kappa
self.theta = theta
self.forcingTermMu = np.nan
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 buildEnergyNormForm(self): # H1
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
mudxM = np.abs(self.mu0(0, 0)) * (fen.dot(self.u.dx(0), self.v.dx(0))
+ fen.dot(self.u, self.v))
imudy = 1. / np.abs(self.mu0(0, 0)) * fen.dot(self.u.dx(1),
self.v.dx(1))
normMatFen = fen.assemble((mudxM + imudy) * fen.dx)
normMat = fen.as_backend_type(normMatFen).mat()
self.energyNormMatrix = scsp.csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
mu = self.checkParameter(mu)
if mu != self.forcingTermMu:
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling base expression for "
"forcing term."),
timestamp = self.timestamp)
pi = np.pi
c, s = np.cos(self.theta), np.sin(self.theta)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0))
mu2R, mu2I = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 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))
self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI]
self.forcingTermMu = mu
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling base expression.",
timestamp = self.timestamp)
return self.forcingTerm
def getExtraFactorB(self, mu : paramVal = [],
derI : int = 0) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
mu = self.checkParameter(mu)
def getPowMinusj(x, power):
powR = x ** power
powI = fenZERO
if power % 2 == 1:
powR, powI = powI, powR
if (power + 3) % 4 < 2:
powR, powI = - powR, - powI
return powR, powI
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
from scipy.special import factorial as fact
y = fen.SpatialCoordinate(self.V.mesh())[1]
powR, powI = [(self.kappa * np.sin(self.theta)) ** derI * k\
for k in getPowMinusj(y, derI)]
mu2R, mu2I = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.)
exprR = mu2R * powR - mu2I * powI
exprI = mu2I * powR + mu2R * powI
if derI >= 1:
muR, muI = np.real(2. * mu(0, 0)), np.imag(2. * mu(0, 0))
powR, powI = [(self.kappa * np.sin(self.theta)) ** (derI - 1) * k\
* derI for k in getPowMinusj(y, derI - 1)]
exprR += muR * powR - muI * powI
exprI += muI * powR + muR * powI
if derI >= 2:
powR, powI = [(self.kappa * np.sin(self.theta)) ** (derI - 2) * k\
* derI * (derI - 1) for k in getPowMinusj(y, derI - 2)]
exprR += powR
exprI += powI
fac = fact(derI)
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return [exprR / fac, exprI / fac]
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
self.As[1] = self.checkAInBounds(-1)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
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
- A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
- A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A2Re)
- DirichletBC0.zero(A2Im)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2ImMat = fen.as_backend_type(A2Im).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
- self.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size)
- + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
- shape = A2ImMat.size))
+ self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
if bs[j] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j < self.nbs:
fRe, fIm = self.getForcingTerm(self.mu0)
cRe, cIm = self.getExtraFactorB(self.mu0, j)
cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(j)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(j)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
-
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- DirichletBC0.apply(b0Re)
- DirichletBC0.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1))
if homogeneized:
Ader = self.A(self.mu0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
index 1cecb1a..f152950 100644
--- a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
@@ -1,320 +1,307 @@
# 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 scsp
import fenics as fen
from rrompy.hfengines.base.problem_engine_base import ProblemEngineBase
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO, fenONE, H1NormMatrix
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.parameter import checkParameter
+from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LaplaceBaseProblemEngine']
class LaplaceBaseProblemEngine(ProblemEngineBase):
"""
Solver for generic Laplace problems.
- \nabla \cdot (a \nabla u) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.mu0 = checkParameter(mu0)
self.npar = self.mu0.shape[1]
self.omega = np.abs(self.mu0(0, 0)) if self.npar > 0 else 0.
self.diffusivity = fenONE
self.forcingTerm = fenZERO
self.DirichletDatum = fenZERO
self.NeumannDatum = fenZERO
self.RobinDatumG = fenZERO
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
ProblemEngineBase.V.fset(self, V)
self.dsToBeSet = True
@property
def diffusivity(self):
"""Value of a."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
self.resetAs()
if not isinstance(diffusivity, (list, tuple,)):
diffusivity = [diffusivity, fenZERO]
self._diffusivity = diffusivity
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm, fenZERO]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
if not isinstance(DirichletDatum, (list, tuple,)):
DirichletDatum = [DirichletDatum, fenZERO]
self._DirichletDatum = DirichletDatum
@property
def NeumannDatum(self):
"""Value of g1."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
self.resetbs()
if not isinstance(NeumannDatum, (list, tuple,)):
NeumannDatum = [NeumannDatum, fenZERO]
self._NeumannDatum = NeumannDatum
@property
def RobinDatumG(self):
"""Value of g2."""
return self._RobinDatumG
@RobinDatumG.setter
def RobinDatumG(self, RobinDatumG):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG, fenZERO]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
if self.verbosity >= 20:
verbosityDepth("INIT", "Initializing boundary measures.",
timestamp = self.timestamp)
mesh = self.V.mesh()
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = mesh,
subdomain_data = boundary_markers)
self.dsToBeSet = False
if self.verbosity >= 20:
verbosityDepth("DEL", "Done initializing boundary measures.",
timestamp = self.timestamp)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = H1NormMatrix(self.V, np.abs(self.omega)**2)
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.As[0] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
if bs[j] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
termNames, terms = [], []
if j == 0:
u0Re, u0Im = self.DirichletDatum
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
termNames += ["forcingTerm", "NeumannDatum", "RobinDatumG"]
terms += [[fRe, fIm], [g1Re, g1Im], [g2Re, g2Im]]
else:
u0Re, u0Im = fenZERO, fenZERO
fRe, fIm = fenZERO, fenZERO
g1Re, g1Im = fenZERO, fenZERO
g2Re, g2Im = fenZERO, fenZERO
if len(termNames) > 0:
parsRe = self.iterReduceQuadratureDegree(zip(
[term[0] for term in terms],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[term[1] for term in terms],
[x + "Imag" for x in termNames]))
else:
parsRe, parsIm = {}, {}
L0Re = (fen.dot(fRe, self.v) * fen.dx
+ fen.dot(g1Re, self.v) * self.ds(0)
+ fen.dot(g2Re, self.v) * self.ds(1))
L0Im = (fen.dot(fIm, self.v) * fen.dx
+ fen.dot(g1Im, self.v) * self.ds(0)
+ fen.dot(g2Im, self.v) * self.ds(1))
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
- DBCR.apply(b0Re)
- DBCI.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if homogeneized:
Ader = self.A(self.mu0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
index 4da30a2..227c446 100644
--- a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
+++ b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
@@ -1,160 +1,158 @@
# 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 Np1D, Tuple, FenExpr, paramVal
from .laplace_base_problem_engine import LaplaceBaseProblemEngine
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
+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 = 20
self.computebsFactors()
self.forcingTermMu = np.nan
import mshr
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
mu = self.checkParameter(mu)
if mu != self.forcingTermMu:
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling base expression for "
"forcing term."),
timestamp = self.timestamp)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
C = np.exp(-.5 * mu(0, 0) ** 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(mu(0, 0)), np.imag(mu(0, 0))
f1R = fen.exp(muR * x) * fen.cos(muI * x)
f1I = fen.exp(muR * x) * fen.sin(muI * x)
self.forcingTerm = [f0 * (CR * f1R - CI * f1I),
f0 * (CR * f1I + CI * f1R)]
self.forcingTermMu = mu
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling base expression.",
timestamp = self.timestamp)
return self.forcingTerm
def computebsFactors(self):
self.bsFactors = np.zeros((self.nbs, self.nbs), dtype = float)
self.bsFactors[0, 0] = 1.
self.bsFactors[1, 1] = 1.
for j in range(2, self.nbs):
l = (j + 1) % 2 + 1
J = np.arange(l, j + 1, 2)
self.bsFactors[j, J] = self.bsFactors[j - 1, J - 1]
if l == 2:
l = 0
J = np.arange(l, j, 2)
self.bsFactors[j, J] += np.multiply(- 1 - J,
self.bsFactors[j - 1, J + 1])
self.bsFactors[j, l : j + 2 : 2] /= j
def getExtraFactorB(self, mu : paramVal = [],
derI : int = 0) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
mu = self.checkParameter(mu)
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0))
x = fen.SpatialCoordinate(self.V.mesh())[0]
l = derI % 2
if l == 0:
powR, powI = fenONE, fenZERO
else:
powR, powI = x - muR, fen.Constant(muI)
exprR, exprI = [self.bsFactors[derI, l] * k for k in [powR, powI]]
for j in range(l + 2, derI + 1, 2):
for _ in range(2):
powR, powI = (powR * (x - muR) - powI * muI,
powR * muI + powI * (x - muR))
exprR += self.bsFactors[derI, j] * powR
exprI += self.bsFactors[derI, j] * powI
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return[exprR, exprI]
def b(self, mu : paramVal = [], der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
if bs[j] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j < self.nbs:
fRe, fIm = self.getForcingTerm(self.mu0)
cRe, cIm = self.getExtraFactorB(self.mu0, j)
cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(j)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(j)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
- DirichletBC0.apply(b0Re)
- DirichletBC0.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1))
if homogeneized:
Ader = self.A(self.mu0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py b/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py
index 0586874..a7ad7bb 100644
--- a/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py
+++ b/rrompy/hfengines/linear_problem/membrane_fracture_engine_nodomain.py
@@ -1,120 +1,104 @@
# 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 scsp
import fenics as fen
import mshr, ufl
from rrompy.utilities.base.types import ScOp, List, paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+from rrompy.solver.fenics import fenics2Sparse
__all__ = ['MembraneFractureEngineNoDomain']
class MembraneFractureEngineNoDomain(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.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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.As[0] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/linear_problem/scattering_problem_engine.py b/rrompy/hfengines/linear_problem/scattering_problem_engine.py
index 2b62cea..aa99e7c 100644
--- a/rrompy/hfengines/linear_problem/scattering_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/scattering_problem_engine.py
@@ -1,174 +1,150 @@
# 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 numpy import inf
-import scipy.sparse as scsp
import fenics as fen
from rrompy.utilities.base.types import List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO
from rrompy.utilities.base import verbosityDepth
from .helmholtz_problem_engine import HelmholtzProblemEngine
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
from rrompy.utilities.exception_manager import RROMPyWarning
+from rrompy.solver.fenics import fenics2Sparse
__all__ = ['ScatteringProblemEngine']
class ScatteringProblemEngine(HelmholtzProblemEngine):
"""
Solver for scattering problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu +- i omega u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
signR: Sign in ABC.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
A2: Scipy sparse array representation (in CSC format) of A2.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
signR = - 1.
def __init__(self, mu0 : paramVal = [0.], degree_threshold : int = inf,
verbosity : int = 10, timestamp : bool = True):
self.silenceWarnings = True
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
del self.silenceWarnings
self.nAs = 3
self.rescalingExp = [1.]
@property
def RobinDatumH(self):
"""Value of h."""
return self.signR * self.omega
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
if not hasattr(self, "silenceWarnings"):
RROMPyWarning(("Scattering problems do not allow changes of h. "
"Ignoring assignment."))
return
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
parsRe = self.iterReduceQuadratureDegree(zip([aRe],
["diffusivityReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([aIm],
["diffusivityImag"]))
a0Re = aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
a0Im = aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a1 = fen.dot(self.u, self.v) * self.ds(1)
- A1 = fen.assemble(a1)
- DirichletBC0.zero(A1)
- A1Mat = fen.as_backend_type(A1).mat()
- A1r, A1c, A1v = A1Mat.getValuesCSR()
- self.As[1] = self.signR * 1.j * scsp.csr_matrix((A1v, A1c, A1r),
- shape = A1Mat.size)
+ self.As[1] = (self.signR * 1.j
+ * fenics2Sparse(a1, {}, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
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 * fen.dot(self.u, self.v) * fen.dx
a2Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
- A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
- A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A2Re)
- DirichletBC0.zero(A2Im)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2ImMat = fen.as_backend_type(A2Im).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
- self.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size)
- + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
- shape = A2ImMat.size))
+ self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
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
index 4c3aa5d..80c1a2b 100644
--- 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
@@ -1,162 +1,157 @@
# 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 scsp
import fenics as fen
from rrompy.hfengines.vector_linear_problem.\
linear_elasticity_beam_poisson_ratio import LinearElasticityBeamPoissonRatio
from rrompy.solver.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.utilities.base import verbosityDepth
-from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.utilities.exception_manager import RROMPyAssert
+from rrompy.utilities.poly_fitting.polynomial import (
+ hashDerivativeToIdx as hashD)
+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.nAs, self.nbs = 5, 4
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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
- for j in [0, 1, 3]:
+ for j in [1, 3]:
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
+ if derI <= 0 and self.As[0] is None:
+ if self.verbosity >= 20:
+ verbosityDepth("INIT", "Assembling operator term A0.",
+ timestamp = self.timestamp)
+ DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
+ self.DirichletBoundary)
+ a0Re = fen.inner(fenZEROS(2), self.v) * fen.dx
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
+ if self.verbosity >= 20:
+ verbosityDepth("DEL", "Done assembling operator term.",
+ timestamp = self.timestamp)
if derI <= 4 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
- a0Re = fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[2] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ a2Re = 2. * fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
+ self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 4 and self.As[4] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A4.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
- a1Re = fen.div(self.u) * fen.div(self.v) * fen.dx
- A1Re = fen.assemble(a1Re)
- DirichletBC0.zero(A1Re)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- self.As[4] = 2. * (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size,
- dtype = np.complex)
- - self.As[2])
+ a4Re = fen.div(self.u) * fen.div(self.v) * fen.dx
+ self.As[4] = (fenics2Sparse(a4Re, {}, DirichletBC0, 0)
+ - 2. * self.As[2])
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
RROMPyAssert(homogeneized, False, "Homogeneized")
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
- if derI <= 3 and self.bs[0] is None:
+ if derI <= 0 and self.bs[0] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b0.",
timestamp = self.timestamp)
fRe, fIm = self.forcingTerm
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTermReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTermImag"]))
L0Re = fen.inner(fRe, self.v) * fen.dx
L0Im = fen.inner(fIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
- DBCR.apply(b0Re)
- DBCI.apply(b0Im)
- self.bs[0] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ self.bs[0] = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if derI <= 3 and self.bs[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b1.",
timestamp = self.timestamp)
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
- b1Re = fen.assemble(L1Re, form_compiler_parameters = parsRe)
- b1Im = fen.assemble(L1Im, form_compiler_parameters = parsIm)
DBC0 = fen.DirichletBC(self.V, fenZEROS(2), self.DirichletBoundary)
- DBC0.apply(b1Re)
- DBC1.apply(b1Im)
- self.bs[1] = np.array(b1Re[:] + 1.j * b1Im[:], dtype = np.complex)
+ self.bs[1] = (fenics2Vector(L1Re, parsRe, DBC0, 1)
+ + 1.j * fenics2Vector(L1Im, parsIm, DBC0, 1))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
if derI <= 2 and self.bs[2] is None:
self.bs[2] = self.checkbInBounds(-1)
if derI <= 3 and self.bs[3] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b3.",
timestamp = self.timestamp)
self.bs[3] = 2. * self.bs[1]
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized)
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
index 0b9515e..096c4ce 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
@@ -1,161 +1,145 @@
# 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 scsp
import fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.solver.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LinearElasticityBeamPoissonRatio']
class LinearElasticityBeamPoissonRatio(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.nAs, self.nbs = 2, 3
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 A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
- if derI <= 1 and self.As[0] is None:
+ if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
- a0Re = self.E_ * fen.inner(epsilon(self.u),
- epsilon(self.v)) * fen.dx
- A0Re = fen.assemble(a0Re)
- DirichletBC0.apply(A0Re)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size,
- dtype = np.complex)
+ a0Re = 2. * self.E_ * fen.inner(epsilon(self.u),
+ epsilon(self.v)) * fen.dx
+ self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
- a1Re = self.E_ * fen.div(self.u) * fen.div(self.v) * fen.dx
- A1Re = fen.assemble(a1Re)
- DirichletBC0.zero(A1Re)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- self.As[1] = 2. * (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size,
- dtype = np.complex)
- - self.As[0])
+ epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
+ a1Re = self.E_ * (fen.div(self.u) * fen.div(self.v)
+ - 4. * fen.inner(epsilon(self.u),
+ epsilon(self.v))) * fen.dx
+ self.As[1] = fenics2Sparse(a1Re, {}, DirichletBC0, 0)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
RROMPyAssert(homogeneized, False, "Homogeneized")
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.bs[0] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b0.",
timestamp = self.timestamp)
fRe, fIm = self.forcingTerm
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTermReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTermImag"]))
L0Re = fen.inner(fRe, self.v) * fen.dx
L0Im = fen.inner(fIm, self.v) * fen.dx
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
- DBCR.apply(b0Re)
- DBCI.apply(b0Im)
- self.bs[0] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ self.bs[0] = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if derI <= 2 and self.bs[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b1.",
timestamp = self.timestamp)
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
- b1Re = fen.assemble(L1Re, form_compiler_parameters = parsRe)
- b1Im = fen.assemble(L1Im, form_compiler_parameters = parsIm)
DBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
- DBC0.apply(b1Re)
- DBC0.apply(b1Im)
- self.bs[1] = np.array(b1Re[:] + 1.j * b1Im[:], dtype = np.complex)
+ self.bs[1] = (fenics2Vector(L1Re, parsRe, DBC0, 1)
+ + 1.j * fenics2Vector(L1Im, parsIm, DBC0, 1))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
if derI <= 2 and self.bs[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b2.",
timestamp = self.timestamp)
self.bs[2] = 2. * self.bs[1]
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized)
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
index 0290795..aee428e 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
@@ -1,180 +1,160 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.utilities.base.types import List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO, fenZEROS, fenONE, elasticNormMatrix
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+from rrompy.solver.fenics import fenics2Sparse
__all__ = ['LinearElasticityHelmholtzProblemEngine']
class LinearElasticityHelmholtzProblemEngine(LinearElasticityProblemEngine):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * mu^2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, 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.nAs = 2
self.omega = np.abs(self.mu0(0, 0))
self.rho_ = fenONE
self.rescalingExp = [2.]
@property
def rho_(self):
"""Value of rho_."""
return self._rho_
@rho_.setter
def rho_(self, rho_):
self.resetAs()
if not isinstance(rho_, (list, tuple,)):
rho_ = [rho_, fenZERO]
self._rho_ = rho_
def buildEnergyNormForm(self): # energy + omega norm
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
lambda_Re, _ = self.lambda_
mu_Re, _ = self.mu_
r_Re, _ = self.rho_
self.energyNormMatrix = elasticNormMatrix(self.V, lambda_Re, mu_Re,
np.abs(self.omega)**2 * r_Re)
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re],
["rho_Real"]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im],
["rho_Imag"]))
a1Re = - rho_Re * fen.inner(self.u, self.v) * fen.dx
a1Im = - rho_Im * fen.inner(self.u, self.v) * fen.dx
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
index 898e32d..bc24df2 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
@@ -1,205 +1,175 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from .linear_elasticity_helmholtz_problem_engine import \
LinearElasticityHelmholtzProblemEngine
from rrompy.utilities.base.types import List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO, fenZEROS
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
+from rrompy.solver.fenics import fenics2Sparse
__all__ = ['LinearElasticityHelmholtzProblemEngineDamped']
class LinearElasticityHelmholtzProblemEngineDamped(
LinearElasticityHelmholtzProblemEngine):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * (mu^2 - i * eta * mu) * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
eta: Value of eta.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, 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.nAs = 3
self.eta = fenZERO
self.rescalingExp = [1.]
@property
def eta(self):
"""Value of eta."""
return self._eta
@eta.setter
def eta(self, eta):
self.resetAs()
if not isinstance(eta, (list, tuple,)):
eta = [eta, fenZERO]
self._eta = eta
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
eta_Re, eta_Im = self.eta
termNames = ["rho_", "eta"]
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re, eta_Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im, eta_Im],
[x + "Imag" for x in termNames]))
a1Re = - ((eta_Re * rho_Im + eta_Im * rho_Re)
* fen.inner(self.u, self.v)) * fen.dx
a1Im = ((eta_Re * rho_Re - eta_Im * rho_Im)
* fen.inner(self.u, self.v)) * fen.dx
- A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
- A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.As[1] = (csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
+ self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re],
["rho_Real"]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im],
["rho_Imag"]))
a2Re = - rho_Re * fen.inner(self.u, self.v) * fen.dx
a2Im = - rho_Im * fen.inner(self.u, self.v) * fen.dx
- A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
- A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
- DirichletBC0.zero(A2Re)
- DirichletBC0.zero(A2Im)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2ImMat = fen.as_backend_type(A2Im).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
- self.As[2] = (csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size)
- + 1.j * csr_matrix((A2Imv, A2Imc, A2Imr),
- shape = A2ImMat.size))
+ self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
index 5fe95b0..a265875 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
@@ -1,346 +1,333 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from rrompy.hfengines.base.vector_problem_engine_base import \
VectorProblemEngineBase
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.solver.fenics import fenZERO, fenZEROS, fenONE, elasticNormMatrix
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.parameter import checkParameter
+from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LinearElasticityProblemEngine']
class LinearElasticityProblemEngine(VectorProblemEngineBase):
"""
Solver for generic linear elasticity problems.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.lambda_ = fenONE
self.mu_ = fenONE
self.mu0 = checkParameter(mu0)
self.npar = self.mu0.shape[1]
self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
VectorProblemEngineBase.V.fset(self, V)
self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.dsToBeSet = True
@property
def lambda_(self):
"""Value of lambda_."""
return self._lambda_
@lambda_.setter
def lambda_(self, lambda_):
self.resetAs()
if not isinstance(lambda_, (list, tuple,)):
lambda_ = [lambda_, fenZERO]
self._lambda_ = lambda_
@property
def mu_(self):
"""Value of mu_."""
return self._mu_
@mu_.setter
def mu_(self, mu_):
self.resetAs()
if not isinstance(mu_, (list, tuple,)):
mu_ = [mu_, fenZERO]
self._mu_ = mu_
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm,
fenZEROS(self.V.mesh().topology().dim())]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
if not isinstance(DirichletDatum, (list, tuple,)):
DirichletDatum = [DirichletDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._DirichletDatum = DirichletDatum
@property
def NeumannDatum(self):
"""Value of g1."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
self.resetbs()
if not isinstance(NeumannDatum, (list, tuple,)):
NeumannDatum = [NeumannDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._NeumannDatum = NeumannDatum
@property
def RobinDatumG(self):
"""Value of g2."""
return self._RobinDatumG
@RobinDatumG.setter
def RobinDatumG(self, RobinDatumG):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG,
fenZEROS(self.V.mesh().topology().dim())]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
if self.verbosity >= 20:
verbosityDepth("INIT", "Initializing boundary measures.",
timestamp = self.timestamp)
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", self.V.mesh(),
self.V.mesh().topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = self.V.mesh(),
subdomain_data = boundary_markers)
self.dsToBeSet = False
if self.verbosity >= 20:
verbosityDepth("DEL", "Done initializing boundary measures.",
timestamp = self.timestamp)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
lambda_Re, _ = self.lambda_
mu_Re, _ = self.mu_
self.energyNormMatrix = elasticNormMatrix(self.V, lambda_Re, mu_Re)
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
- A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
- A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
+ self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu != self.mu0BC:
self.liftDirichletData(self.mu)
fenZEROSEff = fenZEROS(self.V.mesh().topology().dim())
for j in range(derI, nbsTot):
if bs[j] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j == 0:
u0Re, u0Im = self.DirichletDatum
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
else:
u0Re, u0Im = fenZEROSEff, fenZEROSEff
fRe, fIm = fenZEROSEff, fenZEROSEff
g1Re, g1Im = fenZEROSEff, fenZEROSEff
g2Re, g2Im = fenZEROSEff, fenZEROSEff
termNames = ["forcingTerm", "NeumannDatum", "RobinDatumG"]
parsRe = self.iterReduceQuadratureDegree(zip(
[fRe, g1Re, g2Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[fIm, g1Im, g2Im],
[x + "Imag" for x in termNames]))
L0Re = (fen.inner(fRe, self.v) * fen.dx
+ fen.inner(g1Re, self.v) * self.ds(0)
+ fen.inner(g2Re, self.v) * self.ds(1))
L0Im = (fen.inner(fIm, self.v) * fen.dx
+ fen.inner(g1Im, self.v) * self.ds(0)
+ fen.inner(g2Im, self.v) * self.ds(1))
- b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
- b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
- DBCR.apply(b0Re)
- DBCI.apply(b0Im)
- b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
+ b = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if homogeneized:
Ader = self.A(self.mu0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)
diff --git a/rrompy/solver/fenics/__init__.py b/rrompy/solver/fenics/__init__.py
index db6dc54..8e5e891 100644
--- a/rrompy/solver/fenics/__init__.py
+++ b/rrompy/solver/fenics/__init__.py
@@ -1,43 +1,46 @@
# 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 .fenics_constants import (fenZERO, fenZEROS, fenONE, fenONES,
bdrTrue, bdrFalse)
+from .fenics_la import fenics2Sparse, fenics2Vector
from .fenics_norms import (L2NormMatrix, H1NormMatrix, Hminus1NormMatrix,
elasticNormMatrix, elasticDualNormMatrix)
from .fenics_plotting import fenplot, affine_warping
from .fenics_projecting import interp_project
__all__ = [
'fenZERO',
'fenZEROS',
'fenONE',
'fenONES',
'bdrTrue',
'bdrFalse',
+ 'fenics2Sparse',
+ 'fenics2Vector',
'L2NormMatrix',
'H1NormMatrix',
'Hminus1NormMatrix',
'elasticNormMatrix',
'elasticDualNormMatrix',
'fenplot',
'affine_warping',
'interp_project'
]
diff --git a/rrompy/solver/fenics/fenics_matrices.py b/rrompy/solver/fenics/fenics_la.py
similarity index 64%
rename from rrompy/solver/fenics/fenics_matrices.py
rename to rrompy/solver/fenics/fenics_la.py
index 9f38974..d2bace5 100644
--- a/rrompy/solver/fenics/fenics_matrices.py
+++ b/rrompy/solver/fenics/fenics_la.py
@@ -1,35 +1,45 @@
# 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 numpy import array, complex
from scipy.sparse import csr_matrix
import fenics as fen
from rrompy.utilities.base.types import Np2D, FenBC, FenExpr
-__all__ = ['fenics2Sparse']
+__all__ = ['fenics2Sparse', 'fenics2Vector']
def fenics2Sparse(expr:FenExpr, formCompPars : dict = {}, DBC : FenBC = None,
- zero : bool = False) -> Np2D:
+ BCType : int = -1, dtype = complex) -> Np2D:
assembled = fen.assemble(expr, form_compiler_parameters = formCompPars)
- if zero:
+ if BCType == 0:
DBC.zero(assembled)
- else:
+ elif BCType >= 1:
DBC.apply(assembled)
emat = fen.as_backend_type(assembled).mat()
er, ec, ev = emat.getValuesCSR()
- return csr_matrix((ev, ec, er), shape = emat.size)
+ return csr_matrix((ev, ec, er), shape = emat.size, dtype = dtype)
+
+def fenics2Vector(expr:FenExpr, formCompPars : dict = {}, DBC : FenBC = None,
+ BCType : int = -1, dtype = complex) -> Np2D:
+ assembled = fen.assemble(expr, form_compiler_parameters = formCompPars)
+ if BCType == 0:
+ DBC.zero(assembled)
+ elif BCType >= 1:
+ DBC.apply(assembled)
+ return array(assembled[:], dtype = dtype)
diff --git a/rrompy/solver/fenics/fenics_norms.py b/rrompy/solver/fenics/fenics_norms.py
index a4e36c8..96f2c24 100644
--- a/rrompy/solver/fenics/fenics_norms.py
+++ b/rrompy/solver/fenics/fenics_norms.py
@@ -1,76 +1,72 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from rrompy.utilities.base.types import Np2D, FenFunc, DictAny, FenFuncSpace
from rrompy.solver import Np2DLikeInv, Np2DLikeInvLowRank
+from .fenics_la import fenics2Sparse
__all__ = ['L2NormMatrix', 'H1NormMatrix', 'Hminus1NormMatrix',
'elasticNormMatrix', 'elasticDualNormMatrix']
-def _fen2sp(expr):
- matFen = fen.as_backend_type(fen.assemble(expr)).mat()
- return csr_matrix(matFen.getValuesCSR()[::-1], shape = matFen.size)
-
def L2NormMatrix(V:FenFuncSpace, r_ : FenFunc = 1.) -> Np2D:
u = fen.TrialFunction(V)
v = fen.TestFunction(V)
- return _fen2sp(r_ * fen.dot(u, v) * fen.dx)
+ return fenics2Sparse(r_ * fen.dot(u, v) * fen.dx)
def H1NormMatrix(V:FenFuncSpace, w : float = 0., r_ : FenFunc = 1.,
a_ : FenFunc = 1.) -> Np2D:
u = fen.TrialFunction(V)
v = fen.TestFunction(V)
- return _fen2sp((w * r_ * fen.dot(u, v)
- + fen.dot(a_ * fen.grad(u), fen.grad(v))) * fen.dx)
+ return fenics2Sparse((w * r_ * fen.dot(u, v)
+ + fen.dot(a_ * fen.grad(u), fen.grad(v))) * fen.dx)
def Hminus1NormMatrix(V:FenFuncSpace, w : float = 0., r_ : FenFunc = 1.,
a_ : FenFunc = 1., solverType : str = "SPSOLVE",
solverArgs : DictAny = {}, compressRank : int = None,
compressOversampling : int = 10,
compressSeed : int = 420) -> Np2D:
if compressRank is None:
return Np2DLikeInv(H1NormMatrix(V, w, r_, a_), L2NormMatrix(V, r_),
solverType, solverArgs)
return Np2DLikeInvLowRank(H1NormMatrix(V, w, r_, a_), L2NormMatrix(V, r_),
solverType, solverArgs, compressRank,
compressOversampling, compressSeed)
def elasticNormMatrix(V:FenFuncSpace, l_:FenFunc, m_:FenFunc,
w : float = 0., r_ : FenFunc = 1.) -> Np2D:
u = fen.TrialFunction(V)
v = fen.TestFunction(V)
epsilon = lambda f: 0.5 * (fen.grad(f) + fen.nabla_grad(f))
sigma = (l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
- return _fen2sp((w * r_ * fen.dot(u, v)
- + fen.inner(sigma, epsilon(v))) * fen.dx)
+ return fenics2Sparse((w * r_ * fen.dot(u, v)
+ + fen.inner(sigma, epsilon(v))) * fen.dx)
def elasticDualNormMatrix(V:FenFuncSpace, l_:FenFunc, m_:FenFunc,
w : float = 0., solverType : str = "SPSOLVE",
solverArgs : DictAny = {}, r_ : FenFunc = 1.,
compressRank : int = None,
compressOversampling : int = 10,
compressSeed : int = 420) -> Np2D:
if compressRank is None:
return Np2DLikeInv(elasticNormMatrix(V, l_, m_, w, r_),
L2NormMatrix(V, r_), solverType, solverArgs)
return Np2DLikeInvLowRank(elasticNormMatrix(V, l_, m_, w, r_),
L2NormMatrix(V, r_), solverType, solverArgs,
compressRank, compressOversampling, compressSeed)
diff --git a/tests/hfengines/linear_elasticity.py b/tests/hfengines/linear_elasticity.py
index 0ef3355..92e10da 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 rod_3d import rod3Dsolver
def test_elastic_beam():
solver = LinearElasticityBeamPoissonRatio(n = 10, rho_ = 1e3, g = 3,
E = 1e6, nu0 = .45, length = 5, verbosity = 0)
mu = .45
uh = solver.solve(mu)[0]
- assert np.isclose(solver.norm(uh), 7.637337113310191e-08, rtol = 1e-1)
+ assert np.isclose(solver.norm(uh), 5.866888537228743e-08, rtol = 1e-1)
assert np.isclose(solver.norm(solver.residual(uh, mu)[0]), 8.4545952e-13,
rtol = 1e-1)
def test_elastic_rod():
solver = rod3Dsolver()
uh = solver.solve()[0]
assert np.isclose(solver.norm(uh), 0.15563476339534466, rtol = 1e-5)
assert np.isclose(solver.norm(solver.residual(uh)[0]), 5.708389944e-08,
rtol = 1e-1)