Page MenuHomec4science
Files F60126809

linear_elasticity_problem_engine.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 17:18

linear_elasticity_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
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, ScOp
from rrompy.utilities.base.fenics import fenZERO, fenZEROS, fenONE
from rrompy.utilities.base import verbosityDepth
__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.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
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, 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.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): # energy norm
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling energy matrix.",
timestamp = self.timestamp)
lambda_Re, _ = self.lambda_
mu_Re, _ = self.mu_
termNames = ["lambda_", "mu_"]
pars = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re],
[x + "Real" 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))
normMatFen = fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
normMatFen = fen.assemble(normMatFen, form_compiler_parameters = pars)
normMat = fen.as_backend_type(normMatFen).mat()
self.energyNormMatrix = csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling energy matrix.",
timestamp = self.timestamp)
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if 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))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self.As[0]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if homogeneized and not np.isclose(self.mu0BC, mu):
self.u0BC = self.liftDirichletData(mu)
if (max(self.nbs, self.nAs * homogeneized) > 1
and not np.isclose(self.bsmu, mu)):
self.bsmu = mu
self.resetbs()
if self.bs[homogeneized][der] is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(der),
timestamp = self.timestamp)
if der < self.nbs:
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
else:
fRe = fenZEROS(self.V.mesh().topology().dim())
fIm = fenZEROS(self.V.mesh().topology().dim())
g1Re = fenZEROS(self.V.mesh().topology().dim())
g1Im = fenZEROS(self.V.mesh().topology().dim())
g2Re = fenZEROS(self.V.mesh().topology().dim())
g2Im = fenZEROS(self.V.mesh().topology().dim())
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)
if homogeneized:
Ader = self.A(mu, der)
b0Re[:] -= np.real(Ader.dot(self.u0BC))
b0Im[:] -= np.imag(Ader.dot(self.u0BC))
DBCR = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
else:
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[homogeneized][der] = np.array(b0Re[:]
+ 1.j * b0Im[:], dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self.bs[homogeneized][der]

Event Timeline

c4science · Help