### Incorrect Boundary Reaction

117
views
-1
3 months ago by
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):

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
Can you reproduce your problem in a simpler way?  I don't have time to figure out what you're doing from that jumble of code, sorry.
written 3 months ago by pf4d
Hi,
I have updated the code. See if this helps in understanding it better. Thanks for the patience.
written 3 months ago by Vivek Kumar

2
3 months ago by
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):
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))​