error while running a simple linear elasticity problem

Asked by Niraj Kumar Jha on 2013-04-26

hello I am new to Fenics and tried to run a simple example..it shows error

from dolfin import *

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create mesh and define function space

mesh = dolfin.UnitSquare(6, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary conditions

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

# Define Dirichlet boundary (x = 0)
c = Expression(("0.0", "0.0"))
bc = DirichletBC(V, c, left)

# Define functions

v = TestFunction(V) # Test function
U = TrialFunction(V)
f = Constant((0.0, -0.5)) # force on boundary

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

def epsilon(v):return 0.5*(grad(v)+transp(grad(v)))
def sigma(v):return 2*mu*epsilon(v)+lmbda*mult(trace(epsilon(v)),Identity(len(v)))

# Define variational problem

a= dot(grad(v),sigma(U))*dx
L= dot(v,f)*dx

# Compute solution

U = Function(V)
solve(a == L, U, bc)

# Plot solution and mesh

plot(U)

# Dump solution to file in VTK format
file = File('elasticity.pvd')
file << U

# Hold plot
interactive()

any suggestions???

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Blechta
Solved:
2013-04-26
Last query:
2013-04-26
Last reply:
2013-04-26
Best Jan Blechta (blechta) said : #1

On Fri, 26 Apr 2013 13:26:14 -0000
Niraj Kumar Jha <email address hidden> wrote:
> New question #227545 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227545
>
> hello I am new to Fenics and tried to run a simple example..it shows
> error
>
>
> from dolfin import *
>
> # Optimization options for the form compiler
> parameters["form_compiler"]["cpp_optimize"] = True
> ffc_options = {"optimize": True, \
> "eliminate_zeros": True, \
> "precompute_basis_const": True, \
> "precompute_ip_const": True}
>
> # Create mesh and define function space
>
> mesh = dolfin.UnitSquare(6, 4)
> V = VectorFunctionSpace(mesh, "Lagrange", 1)
>
> # Define boundary conditions
>
> left = ["(std::abs(x[0]) < DOLFIN_EPS) && on_boundary"]
left = "(std::abs(x[0]) < DOLFIN_EPS) && on_boundary"
>
> # Define Dirichlet boundary (x = 0)
> c = Expression(("0.0", "0.0"))
> bc = DirichletBC(V, c, left)
>
> # Define functions
>
> v = TestFunction(V) # Test function
> U = TrialFunction(V)
> f = Constant((0.0, -0.5)) # force on boundary
>
> # Elasticity parameters
> E, nu = 10.0, 0.3
> mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 -
> 2*nu)))
>
> def epsilon(v):return 0.5*(grad(v)+transp(grad(v)))
def epsilon(v):return 0.5*(grad(v)+transpose(grad(v)))
> def sigma(v):return
> 2*mu*epsilon(v)+lmbda*mult(trace(epsilon(v)),Identity(len(v)))
2*mu*epsilon(v)+lmbda*tr(epsilon(v))*Identity(V.cell().d)
>
> # Define variational problem
>
> a= dot(grad(v),sigma(U))*dx
a = inner(grad(v),sigma(U))*dx
> L= dot(v,f)*dx
>
>
> # Compute solution
>
> U = Function(V)
> solve(a == L, U, bc)
>
> # Plot solution and mesh
>
> plot(U)
>
> # Dump solution to file in VTK format
> file = File('elasticity.pvd')
> file << U
>
> # Hold plot
> interactive()
>
> any suggestions???
>
>

Niraj Kumar Jha (niraj-jha) said : #2

Thanks Jan Blechta, that solved my question.