Page MenuHomec4science
Files F60484316

solver.py
No OneTemporary
Actions

File Metadata

Created
Tue, Apr 30, 13:55

solver.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
verb = 2
dampingEta = 0 * 1e4 / 2. / np.pi
k = 7773.051993943557 # [Hz]
theta = 20. * np.pi / 180.
phi = 10. * np.pi / 180.
mesh = fen.Mesh("./diapason_1.xml")
subdomains = fen.MeshFunction("size_t", mesh,
"./diapason_1_physical_region.xml")
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
c = 3.e2
rho = 8e3 * (2. * np.pi) ** 2.
E = 1.93e11
nu = .3
T = 1e6
lambda_ = E * nu / (1. + nu) / (1. - 2. * nu)
mu_ = E / (1. + nu)
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(k) / 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_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))
###
if dampingEta > 0:
solver = LEHPED(degree_threshold = 8, verbosity = verb)
solver.eta = dampingEta
else:
solver = LEHPE(degree_threshold = 8, verbosity = verb)
solver.omega = np.real(k)
solver.lambda_ = lambda_eff
solver.mu_ = mu_eff
solver.rho_ = rho_eff
solver.V = fen.VectorFunctionSpace(mesh, "P", 1)
solver.DirichletBoundary = lambda x, on_b: on_b and x[1] < Hball
solver.NeumannBoundary = "REST"
solver.forcingTerm = neumannDatum
uh = solver.solve(k)
solver.outParaviewTimeDomain(uh, omega = 2. * np.pi * k,
filename = "out/outT{}_{}_".format(k, dampingEta),
forceNewFile = False)

Event Timeline

c4science · Help