Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/tests/hfengines/linear_elasticity_beam_poisson_ratio_engine.py
0fc95802a615master

Last Change
Raw File

linear_elasticity_beam_poisson_ratio_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 rrompy.hfengines.vector_linear_problem import \
LinearElasticityProblemEngine
from rrompy.solver.fenics import fenZERO, fenZEROS
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
class LinearElasticityBeamPoissonRatioEngine(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._affinePoly = False
self.nAs, self.nbs = 3, 2
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 buildA(self):
"""Build terms of operator of linear system."""
if self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
a0Re = fenZERO * fen.inner(self.u, self.v) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
self.thAs[0] = self.getMonomialSingleWeight([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(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
a1Re = self.E_ * fen.inner(epsilon(self.u),
epsilon(self.v)) * fen.dx
self.As[1] = fenics2Sparse(a1Re, {}, DirichletBC0, 0)
self.thAs[1] = [('x', '()', 0, '*', -2., '+', 1.),
(-2.0,), None]
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
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(2),
self.DirichletBoundary)
a2Re = self.E_ * fen.div(self.u) * fen.div(self.v) * fen.dx
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
self.thAs[2] = self.getMonomialSingleWeight([1])
vbMng(self, "DEL", "Done assembling operator term.", 20)
def buildb(self):
"""Build terms of operator of linear system."""
if self.thbs[0] is None:
self.thbs = [None] * self.nbs
if self.bs[0] is None:
vbMng(self, "INIT", "Assembling forcing term b0.", 20)
L0Re = fen.inner(fenZEROS(2), self.v) * fen.dx
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
self.bs[0] = (fenics2Vector(L0Re, {}, DBCR, 1)
+ 1.j * fenics2Vector(L0Re, {}, DBCI, 1))
self.thbs[0] = self.getMonomialSingleWeight([0])
if self.bs[1] is None:
vbMng(self, "INIT", "Assembling forcing term b1.", 20)
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
DBC0 = fen.DirichletBC(self.V, fenZEROS(2), self.DirichletBoundary)
self.bs[1] = (fenics2Vector(L1Re, parsRe, DBC0, 1)
+ 1.j * fenics2Vector(L1Im, parsIm, DBC0, 1))
self.thbs[1] = [('x', '()', 0, '**', 2., '*', -2., '-',
('x', '()', 0), '+', 1.),
('x', '()', 0, '*', -4., '-', 1.),
(-2.0,), None]
vbMng(self, "DEL", "Done assembling forcing term.", 20)
c4science · Help