Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py
c6af74908c2cmaster

Last Change
Raw File

helmholtz_base_problem_engine.py
View Options

#!/usr/bin/python
# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
import scipy.sparse as scsp
import fenics as fen
from rrompy.hfengines.fenics.generic_problem_engine import GenericProblemEngine
from rrompy.utilities.types import Np1D, Np2D
from rrompy.utilities.fenics import fenZERO, fenONE, bdrTrue, bdrFalse
__all__ = ['HelmholtzBaseProblemEngine']
class HelmholtzBaseProblemEngine(GenericProblemEngine):
"""
ABSTRACT
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:
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.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
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.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
"""
Aterms = 2
omega = 1.
_diffusivity = [fenONE, fenZERO]
_refractionIndex = [fenONE, fenZERO]
_forcingTerm = [fenZERO, fenZERO]
_DirichletDatum = [fenZERO, fenZERO]
_NeumannDatum = [fenZERO, fenZERO]
_RobinDatumG = [fenZERO, fenZERO]
def __init__(self):
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.DirichletBoundary = "ALL"
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
if hasattr(self, "b0"): del self.b0
if not type(V).__name__ == 'FunctionSpace':
raise Exception("V type not recognized.")
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
self.dsToBeSet = True
@property
def diffusivity(self):
"""Value of a."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
if hasattr(self, "A0"): del self.A0
if not isinstance(diffusivity, (list, tuple,)):
diffusivity = [diffusivity, fenZERO]
self._diffusivity = diffusivity
@property
def refractionIndex(self):
"""Value of n."""
return self._refractionIndex
@refractionIndex.setter
def refractionIndex(self, refractionIndex):
if hasattr(self, "A1"): del self.A1
if not isinstance(refractionIndex, (list, tuple,)):
refractionIndex = [refractionIndex, fenZERO]
self._refractionIndex = refractionIndex
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
if hasattr(self, "b0"): del self.b0
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm, fenZERO]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""
Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
DirichletBCIm.
"""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
if hasattr(self, "b0"): del self.b0
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):
if hasattr(self, "b0"): del self.b0
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):
if hasattr(self, "b0"): del self.b0
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG, fenZERO]
self._RobinDatumG = RobinDatumG
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self._DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
if isinstance(DirichletBoundary, (str,)):
if DirichletBoundary.upper() == "NONE":
self._DirichletBoundary = bdrFalse
self.DREST = 0
elif DirichletBoundary.upper() == "ALL":
self._DirichletBoundary = bdrTrue
self.NeumannBoundary = "NONE"
self.RobinBoundary = "NONE"
self.DREST = 0
elif DirichletBoundary.upper() == "REST":
if self.NREST + self.RREST > 0:
raise Exception("Only 1 'REST' wildcard can be specified.")
self._DirichletBoundary = lambda x, on_boundary : (on_boundary
and not self.NeumannBoundary(x, on_boundary)
and not self.RobinBoundary(x, on_boundary))
self.DREST = 1
else:
raise Exception("DirichletBoundary label not recognized.")
elif callable(DirichletBoundary):
self._DirichletBoundary = DirichletBoundary
self.DREST = 0
else:
raise Exception("DirichletBoundary type not recognized.")
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self._NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
if hasattr(self, "A1"): del self.A1
self.dsToBeSet = True
if isinstance(NeumannBoundary, (str,)):
if NeumannBoundary.upper() == "NONE":
self._NeumannBoundary = bdrFalse
self.NREST = 0
elif NeumannBoundary.upper() == "ALL":
self._NeumannBoundary = bdrTrue
self.DirichletBoundary = "NONE"
self.RobinBoundary = "NONE"
self.NREST = 0
elif NeumannBoundary.upper() == "REST":
if self.DREST + self.RREST > 0:
raise Exception("Only 1 'REST' wildcard can be specified.")
self._NeumannBoundary = lambda x, on_boundary : (on_boundary
and not self.DirichletBoundary(x, on_boundary)
and not self.RobinBoundary(x, on_boundary))
self.NREST = 1
else:
raise Exception("NeumannBoundary label not recognized.")
elif callable(NeumannBoundary):
self._NeumannBoundary = NeumannBoundary
self.NREST = 0
else:
raise Exception("DirichletBoundary type not recognized.")
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self._RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
self.dsToBeSet = True
if isinstance(RobinBoundary, (str,)):
if RobinBoundary.upper() == "NONE":
self._RobinBoundary = bdrFalse
self.RREST = 0
elif RobinBoundary.upper() == "ALL":
self._RobinBoundary = bdrTrue
self.DirichletBoundary = "NONE"
self.NeumannBoundary = "NONE"
self.RREST = 0
elif RobinBoundary.upper() == "REST":
if self.DREST + self.NREST > 0:
raise Exception("Only 1 'REST' wildcard can be specified.")
self._RobinBoundary = lambda x, on_boundary : (on_boundary
and not self.DirichletBoundary(x, on_boundary)
and not self.NeumannBoundary(x, on_boundary))
self.RREST = 1
else:
raise Exception("RobinBoundary label not recognized.")
return
elif callable(RobinBoundary):
self._RobinBoundary = RobinBoundary
self.RREST = 0
else:
raise Exception("RobinBoundary type not recognized.")
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
NB = self.NeumannBoundary
RB = self.RobinBoundary
class NBoundary(fen.SubDomain):
def inside(self, x, on_boundary):
return NB(x, on_boundary)
class RBoundary(fen.SubDomain):
def inside(self, x, on_boundary):
return RB(x, on_boundary)
boundary_markers = fen.MeshFunction("size_t", self.V.mesh(),
self.V.mesh().topology().dim() - 1)
NBoundary().mark(boundary_markers, 0)
RBoundary().mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = self.V.mesh(),
subdomain_data = boundary_markers)
self.dsToBeSet = False
def energyNormMatrix(self) -> Np2D:
"""Sparse matrix (in CSR format) representative of scalar product."""
normMatFen = fen.assemble((fen.dot(fen.grad(self.u), fen.grad(self.v))
+ np.abs(self.omega)**2 * fen.dot(self.u, self.v)) * fen.dx)
normMat = fen.as_backend_type(normMatFen).mat()
return scsp.csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
def liftDirichletData(self) -> Np1D:
"""Lift Dirichlet datum."""
solLRe = fen.interpolate(self.DirichletDatum[0], self.V)
solLIm = fen.interpolate(self.DirichletDatum[1], self.V)
return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector())
def assembleb(self):
"""Assemble matrix blocks of linear system."""
if not hasattr(self, "b0"):
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
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)
b0Im = fen.assemble(L0Im)
fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary).apply(b0Re)
fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary).apply(b0Im)
self.b0 = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
self.bs = [self.b0]
c4science · Help