### Incorrect Boundary Reaction

117

views

-1

I am trying to solve a simple elasticity problem and trying to compute the reaction on the boundary. The code :

```
from dolfin import *
import numpy as np
spd = 2 # Spatial dimension
n = 16 # Size of mesh
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, 'CG', 1)
x = SpatialCoordinate(mesh)
# Define variational problem
du = TrialFunction(V)
delta_u = TestFunction(V)
# set up solution functions
uh = Function(V,name='displacement')
uh_array = uh.vector().array()
uh_array = np.random.rand(len(uh_array))
uh.vector()[:] = uh_array
force = Constant((0., -1.0))
# Define boundary condition
tol = 1E-14
class clamped_boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary())
# set boundary markers for clamped boundary
clamped = clamped_boundary()
domainBoundaries = FacetFunction("size_t", mesh)
domainBoundaries.set_all(0)
clamped.mark(domainBoundaries,1)
ds = Measure("ds")(subdomain_data=domainBoundaries)
clamped_surface_mark = ds(1)
def eps(u):
return sym(grad(u))
lmbda = 10.0
mu = 1.0
def sigma(u):
return lmbda*(tr(eps(u)))*Identity(spd) +2*mu*eps(u)
F_u = inner(sigma(uh),grad(delta_u))*dx + dot(force, delta_u)*dx
J_u = derivative(F_u, uh, du )
# Compute solution
problem = NonlinearVariationalProblem(F_u, uh, bc, J_u)
solver = NonlinearVariationalSolver(problem)
solver.solve()
V = FunctionSpace(mesh, 'CG', 1)
#Compute the boundary forces
n = FacetNormal(mesh)
stress = sigma(uh)
Traction = dot(stress,n)
fx=Traction[0]*clamped_surface_mark
fy=Traction[1]*clamped_surface_mark
print(assemble(fx))
print(assemble(fy))
interactive()
```

However, I am getting wrong values.I am getting:

fy : 1.01141 instead of 1

fx : 0.52087 instead of 0

Any thoughts on what I might be doing wrong?

Thanks

Community: FEniCS Project

### 1 Answer

2

Just raise the degree of the function space to 2 and you'll get the results you want.

Anyway, the code you posted is really a big mess with lot of unnecessary lines. This is much easier to follow (at least for me):

Anyway, the code you posted is really a big mess with lot of unnecessary lines. This is much easier to follow (at least for me):

```
from __future__ import print_function
from fenics import *
mesh = UnitSquareMesh(16, 16)
V = VectorFunctionSpace(mesh, 'P', 2)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
force = Constant((0.0, -1.0))
# Define boundary condition
clamped = CompiledSubDomain('on_boundary && near(x[0], 0.0)')
bc = DirichletBC(V, Constant((0.0, 0.0)), clamped)
# set boundary markers for clamped boundary
domainBoundaries = FacetFunction("size_t", mesh, 0)
clamped.mark(domainBoundaries, 1)
ds = Measure("ds", subdomain_data=domainBoundaries)
lmbda = 10.0
mu = 1.0
def eps(u): return sym(grad(u))
def sigma(u):
return lmbda*(tr(eps(u)))*Identity(2) + 2*mu*eps(u)
a = inner(sigma(u), grad(v))*dx
L = dot(force, v)*dx
# Compute the solution
u = Function(V)
solve(a == L, u, bc)
# Compute the boundary forces
n = FacetNormal(mesh)
Traction = dot(sigma(u), n)
fx = Traction[0]*ds(1)
fy = Traction[1]*ds(1)
print('fx:', assemble(fx))
print('fy:', assemble(fy))
```

Please login to add an answer/comment or follow this question.

I have updated the code. See if this helps in understanding it better. Thanks for the patience.