Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_archway_frequency.py
4c9bfd2ccc95feature_greedy_coarseni…

Last Change
Raw File

linear_elasticity_helmholtz_archway_frequency.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
__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
"""
def __init__(self, kappa:float, n:int, rho_:float, T: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)
import ufl
x, y = fen.SpatialCoordinate(mesh)[:]
NeumannNonZero = ufl.And(ufl.gt(y, r),
ufl.And(ufl.ge(x, -.25 * R), ufl.le(x, .25 * R)))
self.NeumannDatum = [ufl.as_vector((0.,
ufl.conditional(NeumannNonZero,
fen.Constant(T),
0.))),
fenZEROS(2)]
self.DirichletBoundary = lambda x, on_b: on_b and fen.near(x[1], 0.)
self.NeumannBoundary = "REST"
c4science · Help