I have the following UFL based on the UFL defined for plasticity in Fig. 26.4 of fenics book. The purpose is 2D elasto-plastic wave propagation. Before I added the derivative of acceleration term in Bilinear form 'a', the UFL compiled but now it won't. Can someone please advice if this is the right way to express the derivative?
elementA = VectorElement("Lagrange", triangle, 2)
elementT = VectorElement("Quadrature", triangle, 2, 9)
elementS = VectorElement("Quadrature", triangle, 2, 3)
DG0 = FiniteElement("Discontinuous Lagrange", triangle, 0)
w = TestFunction(elementA)
du = TrialFunction(elementA)
f = Coefficient(elementA) # External forces
t = Coefficient(elementT) # consistent tangent moduli
s = Coefficient(elementS) # stress
# Fields from previous time step
u0 = Coefficient(elementA) # displacement at previous step
u = Coefficient(elementA) # current displacement
v0 = Coefficient(elementA) # velocity
a0 = Coefficient(elementA) # acceleration
rho = Coefficient(DG0) # mass density
# Time stepping parameters (Newmark Beta Scheme)
beta = Coefficient(DG0)
gamma = Coefficient(DG0)
dt = Coefficient(DG0)
def eps(u): # eps_xx, eps_yy, gam_xy for 2D plane strain
return as_vector([u[i].dx(i) for i in range(2)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1)]])
def sigma(s): # Stress tensor xx, xy, yx, yy
return as_matrix([[s[0], s[2]], [s[2], s[1]]])
def tangent(t): # Elastoplastic Consistent Tangent
return as_matrix([[t[i*3 + j] for j in range(3)] for i in range(3)])
def ddot_u(u): # Acceleration at t_n+1 from Newmark scheme
return ((u - u0 - dt*v0)/(beta*dt*dt) - a0*(1 - 2*beta)/(2*beta))
RoA = (u - u0 - dt*v0)/(beta*dt*dt) - a0*(1 - 2*beta)/(2*beta)
tripledot_u = derivative(RoA,u,du) # Rate of Acceleration
a = inner(eps(w), dot(tangent(t), eps(u)) )*dx + rho*inner(tripledot_u, w)*dx
L = inner(grad(w), sigma(s))*dx - inner(w, f)*dx + rho*inner(ddot_u(u), w)*dx
