Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/hfengines/fenics_engines/linear_elasticity_helmholtz_problem_engine_augmented.py
d6451a36626adevelop

Last Change
Raw File

linear_elasticity_helmholtz_problem_engine_augmented.py
View Options

# Copyright (C) 2018-2020 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
from scipy.sparse import eye, bmat, block_diag
from .linear_elasticity_helmholtz_problem_engine import (
LinearElasticityHelmholtzProblemEngine,
LinearElasticityHelmholtzProblemEngineDamped)
from rrompy.solver.fenics import (augmentedElasticNormMatrix,
augmentedElasticDualNormMatrix)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import parameterMap as pMap
__all__ = ['LinearElasticityHelmholtzProblemEngineAugmented',
'LinearElasticityHelmholtzProblemEngineDampedAugmented']
class LinearElasticityHelmholtzProblemEngineAugmented(
LinearElasticityHelmholtzProblemEngine):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * mu * v = f in \Omega
mu * u = v in \overline{\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.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
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).
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.parameterMap = pMap(1., self.npar)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self._energyNormMatrix = augmentedElasticNormMatrix(self.V,
self.lambda_[0], self.mu_[0])
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self._energyNormDualMatrix = augmentedElasticDualNormMatrix(
self.V, self.lambda_[0], self.mu_[0],
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
ANone = any([A is None for A in self.As])
if not ANone: return
self.nAs = 2
super().buildA()
I = eye(self.spacedim // 2)
self.As[0] = block_diag((self.As[0], I), format = "csr")
self.As[1] = bmat([[None, self.As[1]], [- I, None]], format = "csr")
def buildb(self):
"""Build terms of operator of linear system."""
bNone = any([b is None for b in self.bs])
if not bNone: return
self.nbs = 1
dim = self.spacedim // 2
super().buildb()
self.bs[0] = np.pad(self.bs[0], (0, dim), "constant")
def plot(self, u, warping = None, is_state = False, name = "u",
save = None, what = 'all', forceNewFile = True,
saveFormat = "eps", saveDPI = 100, show = True, colorMap = "jet",
fenplotArgs = {}, **figspecs):
uh = u[: self.spacedim // 2] if is_state or self.isCEye else u
return super().plot(uh, warping, is_state, name, save, what,
forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
def outParaview(self, u, warping = None, is_state = False, name = "u",
filename = "out", time = 0., what = 'all',
forceNewFile = True, folder = False, filePW = None):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaview(u[: self.spacedim // 2], warping, is_state,
name, filename, time, what, forceNewFile,
folder, filePW)
def outParaviewTimeDomain(self, u, omega, warping = None, is_state = False,
timeFinal = None, periodResolution = 20,
name = "u", filename = "out",
forceNewFile = True, folder = False):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaviewTimeDomain(u[: self.spacedim // 2], omega,
warping, is_state, timeFinal,
periodResolution, name, filename,
forceNewFile, folder)
class LinearElasticityHelmholtzProblemEngineDampedAugmented(
LinearElasticityHelmholtzProblemEngineDamped):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * (mu - i * eta) * v = f in \Omega
mu * u = v in \overline{\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.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
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).
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.nAs = 2
self._weight0 = 1.
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self._energyNormMatrix = augmentedElasticNormMatrix(self.V,
self.lambda_[0], self.mu_[0])
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self._energyNormDualMatrix = augmentedElasticDualNormMatrix(
self.V, self.lambda_[0], self.mu_[0],
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
ANone = any([A is None for A in self.As])
if not ANone: return
self.nAs = 3
super().buildA()
self._nAs = 2
I = eye(self.spacedim // 2)
self.As[0] = bmat([[self.As[0], self._weight0 * self.As[1]],
[None, I]], format = "csr")
self.As[1] = bmat([[(1. - self._weight0) * self.As[1], self.As[2]],
[- I, None]], format = "csr")
self.thAs.pop()
self.As.pop()
def buildb(self):
"""Build terms of operator of linear system."""
bNone = any([b is None for b in self.bs])
if not bNone: return
self.nbs = 1
dim = self.spacedim // 2
super().buildb()
self.bs[0] = np.pad(self.bs[0], (0, dim), "constant")
def plot(self, u, warping = None, is_state = False, name = "u",
save = None, what = 'all', forceNewFile = True,
saveFormat = "eps", saveDPI = 100, show = True, colorMap = "jet",
fenplotArgs = {}, **figspecs):
uh = u[: self.spacedim // 2] if is_state or self.isCEye else u
return super().plot(uh, warping, is_state, name, save, what,
forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
def outParaview(self, u, warping = None, is_state = False, name = "u",
filename = "out", time = 0., what = 'all',
forceNewFile = True, folder = False, filePW = None):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaview(u[: self.spacedim // 2], warping, is_state,
name, filename, time, what, forceNewFile,
folder, filePW)
def outParaviewTimeDomain(self, u, omega, warping = None, is_state = False,
timeFinal = None, periodResolution = 20,
name = "u", filename = "out",
forceNewFile = True, folder = False):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaviewTimeDomain(u[: self.spacedim // 2], omega,
warping, is_state, timeFinal,
periodResolution, name, filename,
forceNewFile, folder)
c4science · Help