Page MenuHomec4science
Files F60883743

linear_elasticity_beam_poisson_ratio.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 03:50

linear_elasticity_beam_poisson_ratio.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 scipy.sparse as scsp
import fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.utilities.base.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, ScOp
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityBeamPoissonRatio']
class LinearElasticityBeamPoissonRatio(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
"""
nAs = 2
nbs = 3
def __init__(self, n:int, rho_:float, g:float, E:float, nu0:float,
length:float, degree_threshold : int = np.inf,
verbosity : int = 10):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity)
self.lambda_ = E * nu0 / (1. + nu0) / (1. - 2 * nu0)
self.mu_ = E / (1. + nu0)
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 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 der <= 1 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.")
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
a0Re = 2 * fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
A0Re = fen.assemble(a0Re)
DirichletBC0.apply(A0Re)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size,
dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.")
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.")
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
a1Re = fen.div(self.u) * fen.div(self.v) * fen.dx
A1Re = fen.assemble(a1Re)
DirichletBC0.apply(A1Re)
A1ReMat = fen.as_backend_type(A1Re).mat()
A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size,
dtype = np.complex)
- 2. * self.As[0])
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.")
if der == 0:
return self.As[0] + mu * self.As[1]
return self.As[1]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
homogeneized = False
bnull = self.checkbInBounds(der)
if bnull is not None: return bnull
if (self.nbs > 1 and not np.isclose(self.bsmu, mu)):
self.bsmu = mu
self.resetbs()
if self.bs[homogeneized][der] is None:
self.autoSetDS()
if self.bs[homogeneized][0] is None and der > 0: self.b(mu, 0)
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b{}.".format(
der))
if der == 0:
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
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
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.bsBase = np.array(b0Re[:] + 1.j * b0Im[:],
dtype = np.complex)
self.bs[homogeneized][0] = (1 - mu - 2 * mu ** 2) * self.bsBase
elif der == 1:
self.bs[homogeneized][der] = (- 1 - 4 * mu) * self.bsBase
elif der == 2:
self.bs[homogeneized][der] = - 2. * self.bsBase
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.")
return self.bs[homogeneized][der]

Event Timeline

c4science · Help