Page MenuHomec4science
Files F61143163

linear_elasticity_beam_elasticity_constants.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 19:45

linear_elasticity_beam_elasticity_constants.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.\
linear_elasticity_beam_poisson_ratio import LinearElasticityBeamPoissonRatio
from rrompy.solver.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LinearElasticityBeamElasticityConstants']
class LinearElasticityBeamElasticityConstants(
LinearElasticityBeamPoissonRatio):
"""
Solver for linear elasticity problem of a beam subject to its own weight,
with parametric Joung modulus and 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, E0:float, nu0:float,
length:float, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(mu0 = [nu0, E0], degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.nAs, self.nbs = 5, 4
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)), fenZEROS(2)]
self.DirichletBoundary = lambda x, on_b: on_b and fen.near(x[0], 0.)
self.NeumannBoundary = "REST"
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
for j in [1, 3]:
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
if derI <= 0 and self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
a0Re = fen.inner(fenZEROS(2), self.v) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 4 and self.As[2] is None:
vbMng(self, "INIT", "Assembling operator term A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
a2Re = 2. * fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 4 and self.As[4] is None:
vbMng(self, "INIT", "Assembling operator term A4.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
a4Re = fen.div(self.u) * fen.div(self.v) * fen.dx
self.As[4] = (fenics2Sparse(a4Re, {}, DirichletBC0, 0)
- 2. * self.As[2])
vbMng(self, "DEL", "Done assembling operator term.", 20)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
RROMPyAssert(homogeneized, False, "Homogeneized")
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if derI <= 0 and self.bs[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling forcing term b0.", 20)
fRe, fIm = self.forcingTerm
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTermReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTermImag"]))
L0Re = fen.inner(fRe, self.v) * fen.dx
L0Im = fen.inner(fIm, 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, parsRe, DBCR, 1)
+ 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if derI <= 3 and 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))
vbMng(self, "DEL", "Done assembling forcing term.", 20)
if derI <= 2 and self.bs[2] is None:
self.bs[2] = self.checkbInBounds(-1)
if derI <= 3 and self.bs[3] is None:
vbMng(self, "INIT", "Assembling forcing term b3.", 20)
self.bs[3] = 2. * self.bs[1]
vbMng(self, "DEL", "Done assembling forcing term.", 20)
return self._assembleb(mu, der, derI, homogeneized)

Event Timeline

c4science · Help