### Simply supported beam with uniform load

65

views

0

Hello All,

I am trying to create a simply supported beam as shown in the following figure:

Code so far:

I am trying to create a simply supported beam as shown in the following figure:

I want to apply load to the top facet of the beam. For this I have so far written the following code but it is giving JIT failure. I am unable to find my mistake and thus will be very grateful if I could please get some help.

Code so far:

```
from __future__ import print_function
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 26 13:49:52 2018
@author: Aby
"""
from fenics import *
#select top facet of the beam
class Top(SubDomain):
def inside(self,x,on_boundary):
return near(x[1],0.2)
top = Top()
L = 1
W = 0.2
mu = 1
lambda_ = 1.25
#creating mesh and function space
mesh = BoxMesh(Point(0,0,0),Point(L,W,W),10,2,2)
V = VectorFunctionSpace(mesh,"P",1)
boundries = FacetFunction("size_t", mesh, mesh.topology().dim() - 1)
boundries.set_all(0)
top.mark(boundries,1)
ds = Measure('ds', domain=mesh, subdomain_data=boundries)
#defining boundary conditions:
tol = 1E-14
def fixed_boundary_left(x,on_boundary):
return on_boundary and x[0]<tol and x[1]<tol
def roller_boundary_bottom(x,on_boundary):
return on_boundary and x[0]==1 and x[1]<tol
bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left)
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom)
bc = [bc_fixed,bc_roller]
def epsilon(u):
strain_u=0.5*(nabla_grad(u)+nabla_grad(u).T)
return strain_u
def sigma(u):
stress_u=lambda_*nabla_div(u)*Identity(d)+2*mu*epsilon(u)
return stress_u
#defining the variational formulation
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()
T = Constant((0,0,0))
f = Constant((0.0,-0.2,0.0))
a = inner(sigma(u),epsilon(v))*dx
L = dot(T,v)*dx+dot(f,v)*ds
#computing the solution
u = Function(V)
solve(a==L,u,bc)
plot(u,title='Displacement',mode='displacement')
vtkfile_u = File('Displacement.pvd')
vtkfile_u << u
```

Community: FEniCS Project

### 1 Answer

3

Running your code, I don't get a JIT error.

I only get the warning that no domains for the BCs

I only get the warning that no domains for the BCs

```
bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left)
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom)
```

are found.

This is because your domain definitions only include one node. Add the pointwise kwarg to the definition

```
bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left, method='pointwise')
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom, method='pointwise')
```

*and* drop the `on_boundary`

since it defaults to `false`

for pointwise BCs:

```
def fixed_boundary_left(x,on_boundary):
return x[0]<tol and x[1]<tol
def roller_boundary_bottom(x,on_boundary):
return x[0]==1. and x[1]<tol
```

1

Yeah I didn't think that through enough. It's obviously all nodes on the edge. But the reasoning still applies (from the docs):

"Available methods are: topological approach (default), geometric approach, and pointwise approach. The topological approach is faster, but will only identify degrees of freedom that are located on a facet that is entirely on the boundary."

"Available methods are: topological approach (default), geometric approach, and pointwise approach. The topological approach is faster, but will only identify degrees of freedom that are located on a facet that is entirely on the boundary."

written
19 days ago by
klunkean

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

Thank you very much for your guidance. It did solve my problem.

I have one more quick question:

I did not understand your comment: "This is because your domain definitions only include one node. "

Will it not include all nodes on the left bottom edge? can you please explain a bit more about this?

Thank you very much.