forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
12.5.2.py
101 lines (101 loc) · 4.16 KB
/
12.5.2.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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
from sympy import Matrix,simplify,diff,integrate,lambdify
import sympy as sp
import matplotlib.pyplot as plt
from mpmath import quad,mp
xi,eta = sp.symbols("xi eta")
coordinates = Matrix([[0,0],[3,1],[35/10,32/10],[5/10,3]])
# coordinates = Matrix([[0,0],[1,0],[1,1],[0,1]])
# plot coordinates
a = Matrix([[0,0],[3,1],[35/10,32/10],[5/10,3],[0,0]])
# a = Matrix([[0,0],[1,0],[1,1],[0,1],[0,0]])
t = 1
Ee = 1
Nu = 3/10
Cc = Ee/(1+Nu)*Matrix([[1/(1 - Nu), Nu/(1 - Nu), 0],
[Nu/(1 - Nu), 1/(1 - Nu), 0],
[0, 0, 1/2]])
fig = plt.figure()
ax = fig.add_subplot(111)
cp = ax.plot(a[:,0], a[:,1])
ax.grid(True, which='both')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Shape")
ax.axhline(y = 0, color = 'k')
ax.axvline(x = 0, color = 'k')
Shapefun = Matrix([[(1-xi)*(1-eta)],[(1+xi)*(1-eta)],
[(1+xi)*(1+eta)],[(1-xi)*(1+eta)]])/4
print("Mapping Functions")
Nn = Matrix([[0 for x in range(8)] for y in range(2)])
for i in range(4):
Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
x = simplify(sum([Shapefun[i]*coordinates[i,0] for i in range(4)]))
print("x =",x)
y = simplify(sum([Shapefun[i]*coordinates[i,1] for i in range(4)]))
print("y =",y)
J = Matrix([[diff(x,xi),diff(x,eta)],[diff(y,xi),diff(y,eta)]])
print("Gradient of the Map")
print("J =",J)
JinvT = J.T.inv()
print("J^-1 =",JinvT)
mat1 = Matrix([[1,0,0,0],[0,0,0,1],[0,1,1,0]])
mat2 = Matrix([[JinvT[0,0],JinvT[0,1],0,0],
[JinvT[1,0],JinvT[1,1],0,0],
[0,0,JinvT[0,0],JinvT[0,1]],
[0,0,JinvT[1,0],JinvT[1,1]]])
mat3 = Matrix([[0 for x in range(8)] for y in range(4)])
for i in range(4):
mat3[0,2*i] = mat3[2,2*i+1] = diff(Shapefun[i],xi)
mat3[1,2*i] = mat3[3,2*i+1] = diff(Shapefun[i],eta)
B = simplify(mat1*mat2*mat3)
print("B =",B)
kbeforeintegration = simplify(B.T*Cc*B)
print("Stiffness Matrix")
K1 = Matrix([[0 for x in range(8)] for y in range(8)])
for i in range(8):
for j in range(8):
integrand = kbeforeintegration[i,j]*J.det()
intlam = lambdify((eta,xi),integrand)
mp.dps = 3
K1[i,j] = quad(intlam,[-1,1],[-1,1])
print("K =",K1)
xiset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
etaset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
kfull = t*Matrix([[simplify(sum([sum([(kbeforeintegration[i,j]*J.det()).subs({xi:k,eta:l}) for k in xiset]) for l in etaset]))for i in range(8)]for j in range(8)])
print("k_{Full Integration} =",kfull)
xiset = Matrix([[0]])
etaset = Matrix([[0]])
kreduced = 4*t*Matrix([[simplify(sum([sum([(kbeforeintegration[i,j]*J.det()).subs({xi:k,eta:l}) for k in xiset]) for l in etaset]))for i in range(8)]for j in range(8)])
print("k_{Reduced Integration}=",kreduced)
print("Body Forces")
rb = Matrix([0,-1])
print("\u03C1_b =",rb)
fbeforeintegration = t*Nn.T*rb*J.det()
f1 = Matrix([integrate(fbeforeintegration[i],(xi,-1,1),(eta,-1,1))for i in range(8)])
print("Body Force Vetor")
print("f_{Exact Integration} =",f1)
xiset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
etaset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
ffull = Matrix([simplify(sum([sum([(fbeforeintegration[i]).subs({xi:l,eta:k}) for k in etaset]) for l in xiset]))for i in range(8)])
print("f_{Full Integration} =",ffull)
xiset = Matrix([[0]])
etaset = Matrix([[0]])
freduced = 4*Matrix([simplify(sum([sum([(fbeforeintegration[i]).subs({xi:l,eta:k}) for k in etaset]) for l in xiset]))for i in range(8)])
print("f_{Reduced Integration}=",freduced)
print("Traction Vectors")
dldeta = sp.sqrt(J[0,1]**2 + J[1,1]**2).subs({xi:1})
tn = Matrix([-1,0])
print("t_n =",tn)
print("Spatial Coordinate System =",dldeta)
fbeforeintegration = t*Nn.T*tn*dldeta
f1 = Matrix([(integrate(fbeforeintegration[i],(eta,-1,1))).subs({xi:1})for i in range(8)])
print("Nodal Forces due to Traction Forces")
print("f_{Exact Integration}",f1)
xiset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
etaset = Matrix([[-1/sp.sqrt(3),1/sp.sqrt(3)]])
ffull = Matrix([simplify(sum([(fbeforeintegration[i]).subs({xi:1,eta:k}) for k in etaset]))for i in range(8)])
print("f_{Full Integration} =",ffull)
xiset = Matrix([[0]])
etaset = Matrix([[0]])
freduced = 2*Matrix([simplify(sum([(fbeforeintegration[i]).subs({xi:1,eta:k}) for k in etaset]))for i in range(8)])
print("f_{Reduced Integration}",freduced)