Hi,
I'm trying to integrate elastodynamic equation using leapfrog algorythm, but my solution with initial velocity is accelerating, instead of "slowing down" due to stifness. Here is my code:
# wave equation in 2D
from dolfin import *
import math
import numpy as np
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 20.0)
mesh = RectangleMesh(Point(0.0, 0.0), Point(20.0, 5.0), 60, 15)
polyDeg = 1
V = VectorFunctionSpace(mesh, "CG", polyDeg)
rho = 10.0
E, nu = 1.0, 0.3
mu, lmbd = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()
I = Identity(d)
def L(q):
return lmbd * tr(q) * I + mu * (q + q.T)
k = inner(L(grad(u)),grad(v))*dx
m = rho * inner(u, v) * dx
f = Expression(('0.0','0.0'))
b = inner(f,v) * dx
T = 10
dt = 0.001
nts = int(math.ceil(T/dt))
u = Function(V)
v = Function(V)
u0 = interpolate(Expression(('0','0')), V)
v0 = interpolate(Expression(('0','0')), V)
u0m = u0.vector().array()
v0m = v0.vector().array()
for initC in range(1700, 1952):
v0m[initC] = 0.1
um1 = u0m
vm1 = v0m
K = assemble(k)
M = assemble(m)
B = assemble(b)
K = -K.array()
B = -B.array()
M = M.array()
Mi = np.linalg.inv(M)
am1 = np.dot(Mi, np.dot(K, u0m))
for ts in range(1, nts) :
t = ts * dt
K = assemble(k)
M = assemble(m)
B = assemble(b)
K = -K.array()
B = -B.array()
M = M.array()
Mi = np.linalg.inv(M)
um = um1 + vm1 * dt + 0.5 * am1 * dt*dt
for dirBC in range(0, 40):
um[dirBC] = 0
u.vector()[:] = um
plot(u, mode = "color", key = "displ", title = str(t))
am = np.dot(Mi, np.dot(K, um))
vm = vm1 + 0.5 * (am1 + am) * dt
print(max(vm))
v.vector()[:] = vm
plot(v, mode = "color", key = "velo", range_max=1.0, range_min=0.0, title = str(t))
um1 = um
am1 = am
vm1 = vm
interactive()