-
Notifications
You must be signed in to change notification settings - Fork 0
/
Elasticity.py
55 lines (45 loc) · 1.19 KB
/
Elasticity.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
from fenics import *
mesh = UnitSquareMesh(2,2)
ndim = mesh.geometry().dim()
print "number of unknown",mesh.num_vertices()
print "number of elements",mesh.num_cells()
plot(mesh, interactive = True)
V = VectorFunctionSpace(mesh,"CG",1)
E = 1.0
nu = 0.3
print E, nu
#E = (E*(1.0+2.0*nu) )/((1.0+nu)**2)
E = E/(1.0-nu**2)
print E
nu = nu/(1.0- nu)
print nu
lmbda = E*nu/((1.0 + nu )*(1.0-2.0*nu))
mu = E/(2.*(1.+nu))
u = TrialFunction(V)
v = TestFunction(V)
class Left(SubDomain):
def inside(self,x,on_boundary):
return (near (x[0],0.0) and on_boundary)
class Right(SubDomain):
def inside(self , x, on_boundary ):
return (near (x[0],1.0) and on_boundary)
right = Right()
left = Left()
bc_fixed = DirichletBC(V, Constant ((0.0,0.0)), left)
disp_right = DirichletBC(V.sub(0), Constant(0.1), right)
bcs = [bc_fixed, disp_right]
def epsilon(v):
return sym(grad(v))
def sigma(v):
return 2.0 * mu * epsilon(v) + lmbda * tr(epsilon(v)) * Identity (ndim)
a = inner(grad(v), sigma(u))*dx
u = Function(V)
solve(lhs(a)==rhs(a),u,bcs)
(u1, u2) = u.split(True)
#plot(u1)
#plot(u2)
#plot(u, mode = 'displacement')
#interactive()
point = (0.5,0.5)
print u.vector().array()
print u(point)