Thank you very much, Prof. Dokken!
Here is the code. It takes from cases library. I want to have a try to using it to conduct dynamic analysis of three-dimensional structures as an exercise and learn more about FENICS by the way. It seems quite challenging now.
from fenics import *
import dolfinx
import matplotlib.pyplot as plt
import meshio
import h5py
from mpi4py import MPI
import numpy as np
import dolfinx.io
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, 'Ldam.xdmf', "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
V = VectorFunctionSpace(mesh, "CG", 2)
# Create classes for defining parts of the boundaries of the domain
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -454.598)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 464.200)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], -250.0)
class Front(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -277.267)
class Back(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 343.795)
# # Initialize sub-domain instances
left = Left()
right = Right()
front = Front()
bottom = Bottom()
back = Back()
# # Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
front.mark(boundaries, 3)
bottom.mark(boundaries, 4)
back.mark(boundaries, 5)
dss = ds(subdomain_data=boundaries)
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, (0.0,0.0,0.0), boundaries, 1),
DirichletBC(V, (0.0,0.0,0.0), boundaries, 2),
DirichletBC(V, (0.0,0.0,0.0), boundaries, 3),
DirichletBC(V, (0.0,0.0,0.0), boundaries, 4),
DirichletBC(V, (0.0,0.0,0.0), boundaries, 5)]
# Parameters
E = 24e9
nu = 0.25
mu = Constant(E / (2.0*(1.0 + nu)))
lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu)))
rho = Constant(2400)
eta_m = Constant(0.1)
eta_k = Constant(0.1)
alpha_m = Constant(0.2)
alpha_f = Constant(0.4)
gamma = Constant(0.5+alpha_f-alpha_m)
beta = Constant((gamma+0.5)**2/4.)
# time
T = 4.0
Nsteps = 200
dt = Constant(T/Nsteps)
time = np.linspace(0, T, Nsteps+1)
p0 = 1.
cutoff_Tc = T/10
p = Expression(("0", "t <= tc ? p0*t/tc : 0", "0"), t=0, tc=cutoff_Tc, p0=p0, degree=0)
u0 = Constant((0.0, 0.0, 0.0))
b = Constant((0.0, 0.0, 0.0))
def epsilon(u_):
return sym(grad(u_))
def sigma(u_):
return 2.0*mu*epsilon(u_) + lmbda*tr(epsilon(u_))*Identity(len(u_))
def m(u_,v_):
return rho*inner(u_, v_)*dx
def k(u_,v_):
return inner(sigma(u_), epsilon(v_))*dx
def c(u_,v_):
return eta_m*m(u_,v_) + eta_k*k(u_,v_)
def f(v_):
return rho*dot(b,v_)*dx+dot(p,v_)*dss(3)
def update_a(u_, u_old_, v_old_, a_old_, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u_-u_old_-dt_*v_old_)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old_
def update_v(a_, u_old_, v_old_, a_old_, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old_ + dt_*((1-gamma_)*a_old_ + gamma_*a_)
def update_fields(u_, u_old_, v_old_, a_old_):
u_vec, u0_vec = u_.vector(), u_old_.vector()
v0_vec, a0_vec = v_old_.vector(), a_old_.vector()
#
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
# (u_old <- u)
v_old_.vector()[:], a_old_.vector()[:] = v_vec, a_vec
u_old_.vector()[:] = u_.vector()
def avg(x_old, x_new, alpha):
return (1-alpha)*x_new + alpha*x_old
# Definition question
u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_old = Function(V)
v_old = Function(V)
a_old = Function(V)
a_new = update_a(u, u_old, v_old, a_old, ufl=True)
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True)
F = m(avg(a_old, a_new, alpha_m),v)+c(avg(v_old, v_new, alpha_f),v)+k(avg(u_old, u, alpha_f),v) - f(v)
a = lhs(F)
L = rhs(F)
#Calculating
vtkfile = File("elasdynamics.pvd")
K, res = assemble_system(a, L, bcs)
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = True
for (i, dt) in enumerate(np.diff(time)):
t = time[i+1]
print("Time: ", t)
p.t = t-float(alpha_f*dt)
# solve
#solve(a == L,u_,bc)
res = assemble(L)
bcs.apply(res)
solver.solve(K, u_.vector(), res)
# save
vtkfile << (u_, t)
# updating
update_fields(u_, u_old, v_old, a_old)