Page MenuHomec4science
Files F60851906

linear_elasticity_helmholtz_archway_frequency_3d.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 23:10

linear_elasticity_helmholtz_archway_frequency_3d.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.Sphere(fen.Point(0, 0, 0), R)
- mshr.Sphere(fen.Point(0, 0, 0), r)
- mshr.Box(fen.Point(-1.05*R, -1.05*R, -1.05*R),
fen.Point(1.05*R, 1.05*R, 0))
- mshr.Box(fen.Point(-1.05*R, -1.05*R, -1.05*R),
fen.Point(1.05*R, -.05*R, 1.05*R))
- mshr.Box(fen.Point(1.05*R, 1.05*R, 1.05*R),
fen.Point(-1.05*R, .05*R, -1.05*R)))
mesh = mshr.generate_mesh(domain, n)
self.V = fen.VectorFunctionSpace(mesh, "P", 1)
import ufl
x, y, z = fen.SpatialCoordinate(mesh)[:]
NeumannNonZero = ufl.And(ufl.gt(z, r),
ufl.And(ufl.ge(x, -.25 * R), ufl.le(x, .25 * R)))
self.NeumannDatum = [ufl.as_vector((0., 0.,fen.Constant(T))),
# ufl.conditional(NeumannNonZero,
# fen.Constant(T),
# 0.))),
fenZEROS(3)]
self.DirichletBoundary = lambda x, on_b: on_b and fen.near(x[2], 0.)
self.NeumannBoundary = "REST"

Event Timeline

c4science · Help