Avoiding assembly with linear elasticity

Question

Hi, I'm trying to avoid assembly in a time loop with linear elasticity. Code demo is below.
I have the displacement function u, and another function c (calculated elsewhere) and define a python function eps=eps(u) which implements voigt notation. Then I want to solve:
eps(u) = c*eps0
where eps0 is a constant vector.

Can anyone please point out where I am going wrong? If I don't use the assemble( ) function line, the answer is wrong.

Thanks!

from dolfin import *

def eps(u):
    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)]])

fcn = Expression('x[0]*.01')

mesh = RectangleMesh(0,0,1,1, 100, 100) #, 'crossed')
V1 = FunctionSpace(mesh, "Lagrange", 1)
V_vec = VectorFunctionSpace(mesh, "Lagrange", 1)

c = Function(V1)

u = TrialFunction(V_vec)
test_u = TestFunction(V_vec)

left = CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && on_boundary")

bc = DirichletBC(V_vec, [0, 0] , left)
eps0 = as_vector((1,1,0))

c.interpolate(fcn)

au = inner(eps(u), eps(test_u))*dx 
Lu = c*inner(eps0, eps(test_u))*dx
u = Function(V_vec) 

m_el = inner(eps0, eps(test_u))*dx
M_el = assemble(m_el)

A = assemble(au)
bc.apply(A)

solver_u = PETScKrylovSolver('cg')
prm = solver_u.parameters
prm["absolute_tolerance"] = 1e-4
prm["relative_tolerance"] = 1e-5
prm["nonzero_initial_guess"] = True

# begin time loop
for i in range(10):
    b = project(as_vector((c,c)),V_vec).vector()*M_el
    b = assemble(Lu)
    bc.apply(b)
    # Other time loop stuff
    solver_u.solve(A, u.vector(), b)
    print i

plot(u, interactive=True)

Comments

Can't you assemble A before entering the time loop? Where is the time loop in your example?

---

Yes, I do assemble A beforehand. I will update the demo code to show the time loop a bit more clearly.

---

Are you looking for a matrix-free solver?

---

Not necessarily no, just to save a bit of time on assembly.

---

At this point I'm lost. You want to solve a linear system of equations, right? If yes then you either assemble the matrix or use a matrix-free iterative solver. In the last case you'll likely have to assemble a preconditioner. What am I missing?

---

Assembling the matrix/ vector can be a costly procedure, and can be avoided by breaking the matrix/vector into operations of other vectors (for example), some of which don't change from iteration to iteration. Take a look here:
http://fenicsproject.org/documentation/tutorial/timedep.html#avoiding-assembly

---

In

b = project(as_vector((c,c)),V_vec).vector()*M_el

you are performing an elementwise product between two vectors. At least, bot M_el and project.vector() are vectors. To get what you want (the right answer for a fcn that varies in space) shouldn't it be a matrix-vector product? Right now you get the right result with a constant fcn, but I'm not sure that it should give you the correct answer for an fcn that varies in space.