Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_elasticity_constants.py
2638d26e1058master

Last Change
Raw File

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 scipy.sparse as scsp
import fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.solver.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, ScOp, paramVal
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['LinearElasticityBeamElasticityConstants']
class LinearElasticityBeamElasticityConstants(LinearElasticityProblemEngine):
"""
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
"""
npar = 2
nAs, nbs = 1, 1
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__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.lambda_ = E0 * nu0 / (1. + nu0) / (1. - 2 * nu0)
self.mu_ = E0 / (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)), 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 : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
Anull = self.checkAInBounds(der)
if der > 0:
raise RROMPyException("Not implemented.")
if Anull is not None: return Anull
self.autoSetDS()
if self.Asmu != mu:
self.Asmu = mu
self.resetAs()
A = self.As[der]
if A is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
m_ = mu(0) / (1. + mu(1))
l_ = m_ * mu(1) / (1. - 2 * mu(1))
lambda_Re, lambda_Im = np.real(l_), np.imag(l_)
mu_Re, mu_Im = np.real(m_), np.imag(m_)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
a0Re = (mu_Re * fen.inner(epsilon(self.u), epsilon(self.v))
+ lambda_Re * fen.div(self.u) * fen.div(self.v)) * fen.dx
a0Im = (mu_Im * fen.inner(epsilon(self.u), epsilon(self.v))
+ lambda_Im * fen.div(self.u) * fen.div(self.v)) * fen.dx
A0Re = fen.assemble(a0Re)
A0Im = fen.assemble(a0Im)
DirichletBC0.apply(A0Re)
DirichletBC0.apply(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()
A = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), shape = A0ImMat.size))
self.As[0] = A
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return A
c4science · Help