Error with interpolation of initial conditions (function space)

Question

Hello,

I have been trying to run a linear elastic wave simulation code (newmark elastodynamics) by imposing an initial displacement (Gaussian pulse) instead of a force loading. However, I am not able to interpolate the initial displacement on the Vectorfunctionspace. If I do that I get an error.

*** Error:   Unable to interpolate function into function space.
*** Reason:  Rank of function (0) does not match rank of function space (1).

When I change the VectorFunctionSpace to just FunctionSpace, the interpolation works but the stress tensor gives me an error.

Symmetric part of tensor with rank != 2 is undefined.

I tried a projection technique given as solution on fenics forum - that solved the interpolation and stress tensor problem but gave me an error further ahead in the newmark scheme where I never had an error before. I am trying to separate the stress in tensorfunctionspace now. Meanwhile, if someone can suggest a quick fix, the minimal code is shown below:

# Define function space
V = VectorFunctionSpace(mesh, "CG", 2)

# Test and trial functions
u1, w = TrialFunction(V), TestFunction(V)

# Fields from previous time step (displacement, velocity, acceleration)
u0, v0, a0 = Function(V), Function(V), Function(V)

# Velocity and acceleration at t_(n+1)
v1 = (gamma/(beta*dt))*(u1 - u0) - (gamma/beta - 1.0)*v0 - dt*(gamma/(2.0*beta) - 1.0)*a0
a1 = (1.0/(beta*dt**2))*(u1 - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0

#Initial conditions
zeros = Expression('0.0')
ui  = Expression('-40.*sin(pi*x[0])*sin(pi*x[1])*exp(-1000.*(pow(x[0]-0.5,2)+pow(x[1]-0.5,2) - c*t))', c=c, t=0.0)
f  = Constant((0.0, 0.0))
u0 = interpolate(ui, V)
#plot(u0, "displacement", interactive = True)
v0 = interpolate(zeros, V)
a0 = interpolate(zeros, V)

# Stress tensor with damping term
def sigma(u, v):
    return 2.0*mu*sym(grad(u)) + (lmbda*tr(grad(u)) + eta*tr(grad(v)))*Identity(mesh.geometry().dim())

# Governing equation
F = (rho*dot(a1, w) + inner(sigma(u1, v1), sym(grad(w))))*dx - dot(f, w)*dx

Answer 1

Hi,
your first problem is that your defined expressions for the initial conditions are scalars. Consider to use:

#Initial conditions
zeros = Expression(("0.0","0.0"))
ui  = Expression(("0.0","-40.*sin(pi*x[0])*sin(pi*x[1])*exp(-1000.*(pow(x[0]-0.5,2)+pow(x[1]-    0.5,2) - c*t))"), c=c, t=0.0)

Whereas your second problem is related to the dimension of grad(u). If u belongs to a scalar function space then grad(u) is a vectorial function, therefore doesn't exist its symmetric part.

Regards.

Comments

Thanks for your reply Hernan. If I modify the expression like that while sticking to a vectorfunctionspace, it gives me a completely different projection which does not propagate anything. You can see that using

plot(u0, "displacement", interactive = True)

in the code. The output of code doesn't show anything.

---

Change the lines

# Velocity and acceleration at t_(n+1)
v1 = (gamma/(beta*dt))*(u1 - u0) - (gamma/beta - 1.0)*v0 - dt*(gamma/(2.0*beta) - 1.0)*a0
a1 = (1.0/(beta*dt**2))*(u1 - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0

# Stress tensor with damping term
def sigma(u, v):
    return 2.0*mu*sym(grad(u)) + (lmbda*tr(grad(u)) + eta*tr(grad(v)))*Identity(mesh.geometry().dim())

# Governing equation
F = (rho*dot(a1, w) + inner(sigma(u1, v1), sym(grad(w))))*dx - dot(f, w)*dx

to

# Velocity and acceleration at t_(n+1)
def dot_u(u):
  return (gamma/(beta*dt))*(u - u0) - (gamma/beta - 1.0)*v0 - dt*(gamma/(2.0*beta) - 1.0)*a0

def ddot_u(u):
  return (1.0/(beta*dt**2))*(u - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0

def sigma(u):
    return 2.0*mu*sym(grad(u)) + (lmbda*tr(grad(u)) + eta*tr(grad(dot_u(u))))*Identity(mesh.geometry().dim())

# Governing equation
F = (rho*inner(ddot_u(u1), w) + inner(sigma(u1), sym(grad(w))))*dx #- dot(f, w)*dx

(the above indications have worked for me).

---

Thanks a lot Hernan!

---

Hi,

This is a very useful post for me. Can you please tell me how can I know that the damping is working? When I set damping zero, i.e., eta=0, and damping 20%, i.e., eta =0.2, I don't see any visible difference in Paraview. Should I try with very high eta and then compare? or should I run both cases for very long time so the waves in damped model should disappear?

Please tell me, how to know the difference between paraview output of eta = 0 and eta = 0.2?

Thank you.