Page MenuHomec4science
Files F61483330

linear_elasticity_helmholtz_problem_engine.py
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 23:20

linear_elasticity_helmholtz_problem_engine.py
View Options

# 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 fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import (fenZERO, fenZEROS, fenONE, fenics2Sparse,
elasticNormMatrix, elasticDualNormMatrix)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.parameter import parameterMap as pMap
__all__ = ['LinearElasticityHelmholtzProblemEngine',
'LinearElasticityHelmholtzProblemEngineDamped']
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.
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, 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._affinePoly = True
self.nAs = 2
self.omega = np.abs(self.mu0(0, 0))
self.rho_ = fenONE
self.parameterMap = pMap([2.] + [1.] * (self.npar - 1))
@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.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self.energyNormMatrix = elasticNormMatrix(
self.V, self.lambda_[0], self.mu_[0],
np.abs(self.omega)**2 * self.rho_[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 = elasticDualNormMatrix(
self.V, self.lambda_[0], self.mu_[0],
np.abs(self.omega)**2 * self.rho_[0],
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
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))
self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[1] is None:
vbMng(self, "INIT", "Assembling operator term A1.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
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
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
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.
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._affinePoly = True
self.nAs = 3
self.eta = fenZERO
self.parameterMap = pMap(1., self.npar)
@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 buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
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))
self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[1] is None:
vbMng(self, "INIT", "Assembling operator term A1.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
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
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[2] is None:
vbMng(self, "INIT", "Assembling operator term A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
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
self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)

Event Timeline

c4science · Help