Stiffness matrix and Reaction Force calculation

Question

Hi, I am solving elasticity problem directly from energy and calculating stiffness matrix and reaction force. But I have got wrong stiffness matrix and reaction force. any idea to fix this issue.

Thanks in advance.

from dolfin import *
mesh = UnitSquareMesh(2,2)
E, nu = Constant(10.0), Constant(0.3)
plot(mesh, title = 'mesh', interactive = True)
ndim = mesh.geometry().dim()
def eps(v):
      return sym(grad(v))
def eps_dev(du):
      return dev(eps(du))
zero_v=Constant((0.,)*ndim)   
f=zero_v
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
du = Function(V)

mu    = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu )*(1.0-2.0*nu))
EDef = ( mu*tr(eps(du)*eps(du)) + 0.5*lmbda*tr(eps(du))**2 )*dx - dot(f,du)*dx
EDer = derivative(EDef,du,v)
J = derivative(EDer,du,u)
class top(SubDomain):
       def inside(self,x,on_boundary):
               tol = 1e-10
                return abs(x[1]-1.0) < tol and on_boundary

class bottom(SubDomain):
       def inside(self,x,on_boundary):
               tol = 1e-10
               return abs(x[1]) < tol and on_boundary
Top = top()
Bottom = bottom()
bclx= DirichletBC(V.sub(0), Constant(0.0), Bottom)
bcly = DirichletBC(V.sub(1), Constant(0.0), Bottom)
bcty = DirichletBC(V.sub(1), Constant(1.0), Top)
bc_u = [bclx, bcly, bcty]

solve(EDer==0,du, J = J, bcs=bc_u)
# siffness matrix calculation
a = lhs(EDer)
A = assemble(EDer)
print A
point = (1.0,1.0)
print 'numerical u at the corner point:', du(point)
plot(du, interactive  = True)
u1, u2 = split(du)
plot(u1, interactive = True)
plot(u2, interactive = True)

Comments

What do you mean by stiffness matrix? Your problem is nonlinear, so there is no such matrix that Ax = b. The only stiffness matrices are for linearized problems.

You are assembling EDer after you compute the solution, du, so the "matrix" A that you compute is just a vector. The only unknown in your EDer is a test function v, therefore you have a vector.

---

How I will get reaction forces ?
K*du =reaction forces ?

---

Force energy is dot(f, du) and its derivative wrt. du in direction v is just dot(f, v). You can get the RHS force vector with print assemble(dot(f, v)*dx).array()

---

Thanks for the quick reply , I will try this.

---

After doing that we will get:
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
but i need reaction force.