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)