Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/diapason/diapason_engine.py
fb0a24fbc16ffeature_greedy_coarseni…

Last Change
Raw File

diapason_engine.py
View Options

import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.vector_linear_problem import \
LinearElasticityHelmholtzProblemEngine as LEHPE
from rrompy.hfengines.vector_linear_problem import \
LinearElasticityHelmholtzProblemEngineDamped as LEHPED
class DiapasonEngine(LEHPE):
def __init__(self, kappa:float, c:float, rho:float, E:float, nu:float,
T:float, theta:float, phi:float, meshNo : int = 1,
deg : int = 1, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = kappa
mesh = fen.Mesh("../data/mesh/diapason_{}.xml".format(meshNo))
subdomains = fen.MeshFunction("size_t", mesh,
("../data/mesh/diapason_{}_physical_"
"region.xml").format(meshNo))
meshBall = fen.SubMesh(mesh, subdomains, 2)
meshFork = fen.SubMesh(mesh, subdomains, 1)
Hball = np.max(meshBall.coordinates()[:, 1]) #.00257
Ltot = np.max(mesh.coordinates()[:, 1]) #.1022
Lhandle = np.max(meshFork.coordinates()[:, 1]) #.026
Lrod = Ltot - Lhandle #.0762
L2var = (Lrod / 4.) ** 2.
Ehandle_ratio = 3.
rhohandle_ratio = 1.5
kWave = (np.cos(theta) * np.cos(phi), np.sin(phi),
np.sin(theta) * np.cos(phi))
x, y, z = fen.SpatialCoordinate(mesh)[:]
yCorr = y - Ltot
compPlane = kWave[0] * x + kWave[1] * y + kWave[2] * z
xPlane, yPlane, zPlane = (xx - compPlane * xx for xx in (x, y, z))
xOrtho, yOrtho, zOrtho = (compPlane * xx for xx in (x, y, z))
forcingBase = (T / (2. * np.pi * L2var)**.5
* fen.exp(- (xPlane**2. + yPlane**2. + zPlane**2.)
/ (2.*L2var)))
forcingWeight = np.real(kappa) / c * (xOrtho + yOrtho + zOrtho)
neumannDatum = [ufl.as_vector(tuple(forcingBase
* fen.cos(forcingWeight)
* kWavedir for kWavedir in kWave)),
ufl.as_vector(tuple(forcingBase
* fen.sin(forcingWeight)
* kWavedir for kWavedir in kWave))]
lambda_ = E * nu / (1. + nu) / (1. - 2. * nu)
mu_ = E / (1. + nu)
lambda_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(lambda_),
fen.Constant(Ehandle_ratio * lambda_))
mu_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(mu_),
fen.Constant(Ehandle_ratio * mu_))
rho_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(rho),
fen.Constant(rhohandle_ratio * rho))
self.lambda_ = lambda_eff
self.mu_ = mu_eff
self.rho_ = rho_eff
self.V = fen.VectorFunctionSpace(mesh, "P", deg)
self.DirichletBoundary = lambda x, on_b: on_b and x[1] < Hball
self.NeumannBoundary = "REST"
self.forcingTerm = neumannDatum
class DiapasonEngineDamped(LEHPED):
def __init__(self, kappa:float, c:float, rho:float, E:float, nu:float,
T:float, theta:float, phi:float, dampingEta:float,
meshNo : int = 1, deg : int = 1,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.eta = dampingEta
self.omega = kappa
mesh = fen.Mesh("../data/mesh/diapason_{}.xml".format(meshNo))
subdomains = fen.MeshFunction("size_t", mesh,
("../data/mesh/diapason_{}_physical_"
"region.xml").format(meshNo))
meshBall = fen.SubMesh(mesh, subdomains, 2)
meshFork = fen.SubMesh(mesh, subdomains, 1)
Hball = np.max(meshBall.coordinates()[:, 1]) #.00257
Ltot = np.max(mesh.coordinates()[:, 1]) #.1022
Lhandle = np.max(meshFork.coordinates()[:, 1]) #.026
Lrod = Ltot - Lhandle #.0762
L2var = (Lrod / 4.) ** 2.
Ehandle_ratio = 3.
rhohandle_ratio = 1.5
kWave = (np.cos(theta) * np.cos(phi), np.sin(phi),
np.sin(theta) * np.cos(phi))
x, y, z = fen.SpatialCoordinate(mesh)[:]
yCorr = y - Ltot
compPlane = kWave[0] * x + kWave[1] * y + kWave[2] * z
xPlane, yPlane, zPlane = (xx - compPlane * xx for xx in (x, y, z))
xOrtho, yOrtho, zOrtho = (compPlane * xx for xx in (x, y, z))
forcingBase = (T / (2. * np.pi * L2var)**.5
* fen.exp(- (xPlane**2. + yPlane**2. + zPlane**2.)
/ (2.*L2var)))
forcingWeight = np.real(kappa) / c * (xOrtho + yOrtho + zOrtho)
neumannDatum = [ufl.as_vector(tuple(forcingBase
* fen.cos(forcingWeight)
* kWavedir for kWavedir in kWave)),
ufl.as_vector(tuple(forcingBase
* fen.sin(forcingWeight)
* kWavedir for kWavedir in kWave))]
lambda_ = E * nu / (1. + nu) / (1. - 2. * nu)
mu_ = E / (1. + nu)
lambda_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(lambda_),
fen.Constant(Ehandle_ratio * lambda_))
mu_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(mu_),
fen.Constant(Ehandle_ratio * mu_))
rho_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(rho),
fen.Constant(rhohandle_ratio * rho))
self.lambda_ = lambda_eff
self.mu_ = mu_eff
self.rho_ = rho_eff
self.V = fen.VectorFunctionSpace(mesh, "P", deg)
self.DirichletBoundary = lambda x, on_b: on_b and x[1] < Hball
self.NeumannBoundary = "REST"
self.forcingTerm = neumannDatum
c4science · Help