Page MenuHomec4science
Files F77790593

linear_elasticity_helmholtz_archway_frequency_broken.py
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 16, 10:29

linear_elasticity_helmholtz_archway_frequency_broken.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 .linear_elasticity_helmholtz_problem_engine import \
LinearElasticityHelmholtzProblemEngine
from rrompy.utilities.base.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityHelmholtzArchwayFrequency']
class LinearElasticityHelmholtzArchwayFrequency(
LinearElasticityHelmholtzProblemEngine):
"""
Solver for archway linear elasticity Helmholtz problem with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * omega^2 * u = rho_ * g / omega in \Omega
u = 0 on \Gamma_D
\partial_nu = 0 on \Gamma_N
"""
nAs = 2
nbs = 20
def __init__(self, kappa:float, n:int, rho_:float, g:float, lambda_:float,
mu_:float, R:float, r:float, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = kappa
self.lambda_ = lambda_
self.mu_ = mu_
self.rho_ = rho_
import mshr
domain = (mshr.Circle(fen.Point(0, 0), R)
- mshr.Circle(fen.Point(0, 0), r)
- mshr.Rectangle(fen.Point(-1.05*R, -1.05*R),
fen.Point(1.05*R, 0)))
mesh = mshr.generate_mesh(domain, n)
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[1], 0.)
self.NeumannBoundary = "REST"
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),
timestamp = self.timestamp)
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)
if np.isclose(mu, 0.):
if der == 0:
self.bs[homogeneized][der] = self.bsBase
else:
self.bs[homogeneized][der] = np.zeros_like(self.bsBase)
else:
self.bs[homogeneized][der] = (mu * np.power(- mu, - der)
* self.bsBase)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self.bs[homogeneized][der]

Event Timeline

c4science · Help