Skip to content
This repository was archived by the owner on Jun 16, 2022. It is now read-only.

Commit f3a858f

Browse files
diehlpkdiehlpk
committed
Add hard loading case
1 parent 341bc09 commit f3a858f

1 file changed

Lines changed: 335 additions & 0 deletions

File tree

‎bondbased2D-plate-hard.py‎

Lines changed: 335 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,335 @@
1+
#Protoype for the 2D bond-based analytic stiffness matrix including fracture
2+
# using a pre-crack square plate
3+
#@author Patrick Diehl (patrickdiehl@lsu.edu)
4+
#@date September 2020
5+
import numpy as np
6+
from scipy import linalg
7+
from shapely.geometry import LineString
8+
import matplotlib.pyplot as plt
9+
import matplotlib as mpl
10+
import sys
11+
import pickle
12+
from matplotlib import rc
13+
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
14+
rc('text', usetex=True)
15+
rc('font', size=14)
16+
17+
np.set_printoptions(precision=4)
18+
#mpl.style.use('seaborn')
19+
20+
class Compute:
21+
22+
nodes = []
23+
neighbors = []
24+
25+
# Config
26+
E=4000
27+
G=E/(2*(1+1/3))
28+
beta=10
29+
C=300000000
30+
rbar=np.sqrt(0.5/beta)
31+
32+
def __init__(self,h,delta_factor,iter):
33+
# Generate the mesh
34+
self.h = h
35+
self.delta_factor = delta_factor
36+
self.delta = self.delta_factor*h
37+
n = int(15/self.h) + 1
38+
self.fix = []
39+
self.loadT = []
40+
self.loadB = []
41+
self.z = []
42+
self.iter = iter
43+
44+
#Generate grid
45+
index = 0
46+
47+
for i in range(0,n):
48+
for j in range(0,n):
49+
self.nodes.append([i*h,j*h])
50+
51+
if j*h < self.delta and i * h < 13 * self.delta :
52+
self.loadB.append(index)
53+
#plt.scatter(i*h,j*h)
54+
55+
if j * h > 15 + h - self.delta and i * h < 13 * self.delta:
56+
self.loadT.append(index)
57+
#plt.scatter(i*h,j*h)
58+
59+
if j*h < self.delta and i * h > 15 + h - self.delta:
60+
self.fix.append(index)
61+
#plt.scatter(i*h,j*h)
62+
63+
if j * h > 15 + h - self.delta and i * h > 15 + h - self.delta:
64+
self.fix.append(index)
65+
#plt.scatter(i*h,j*h)
66+
67+
index += 1
68+
69+
70+
#plt.show()
71+
#sys.exit(1)
72+
73+
self.fix = np.sort(self.fix)
74+
75+
self.nodes = np.array(self.nodes)
76+
77+
self.V = np.empty(len(self.nodes))
78+
self.V.fill(h*h)
79+
80+
self.d = np.zeros(len(self.nodes))
81+
82+
self.VB = np.pi * self.delta * self.delta
83+
# Search neighbors
84+
self.searchNeighbors()
85+
86+
self.f = np.zeros(2*len(self.nodes))
87+
88+
if self.iter == 1:
89+
90+
# Initialize
91+
self.uCurrent = np.zeros(2*len(self.nodes)).reshape((len(self.nodes),2))
92+
93+
else:
94+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(self.iter-1)+"-displacement.npy", "rb")
95+
#self.nodes += np.load(filehandler)
96+
self.uCurrent = np.load(filehandler)
97+
98+
99+
100+
#Apply the load to the body force vector
101+
self.b = np.zeros(2*len(self.nodes))
102+
self.b2 = np.zeros(2*len(self.nodes))
103+
self.b3 = np.zeros(2*len(self.nodes))
104+
for i in range(0,len(self.nodes)):
105+
if i in self.loadT:
106+
self.b[2*i+1] = 4e3 / (2*self.delta*self.delta)
107+
self.b2[2*i+1] = 4e5 / (2*self.delta*self.delta)
108+
self.b3[2*i+1] = 4e6 / (2*self.delta*self.delta)
109+
if i in self.loadB:
110+
self.b[2*i+1] = -(4e3) / (2*self.delta*self.delta)
111+
self.b2[2*i+1] = -(4e5) / (2*self.delta*self.delta)
112+
self.b3[2*i+1] = -(4e6) / (2*self.delta*self.delta)
113+
114+
print("Matrix size "+str(2*len(self.nodes))+"x"+str(2*len(self.nodes)))
115+
116+
117+
def searchNeighbors(self):
118+
left = np.array([0,7.5])
119+
right = np.array([7.5,7.5])
120+
neighbor = []
121+
122+
123+
#fig = plt.gcf()
124+
#fig.set_size_inches(4,50)
125+
126+
for i in range(0,len(self.nodes)):
127+
self.neighbors.append([])
128+
for j in range(0,len(self.nodes)):
129+
if i != j and self.length(j,i) <= self.delta:
130+
if not self.intersect(left,right,self.nodes[i],self.nodes[j]):
131+
self.neighbors[i].append(j)
132+
#plt.plot([self.nodes[i][0],self.nodes[j][0]],[self.nodes[i][1],self.nodes[j][1]],c="#91A3B0",alpha=0.15)
133+
neighbor.append(len(self.neighbors[i]))
134+
135+
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=neighbor)
136+
#plt.xlabel("Position $x$",fontsize = 30)
137+
#plt.ylabel("Position $y$",fontsize = 30)
138+
#clb = plt.colorbar()
139+
#clb.set_label(r'$B_\delta(x)$',labelpad=5)
140+
#plt.show()
141+
#sys.exit()
142+
143+
def L(self,i,j):
144+
r = np.sqrt(self.length(i,j)) * self.S(i,j)
145+
return (2./self.VB) * self.w(self.length(i,j))/self.delta * self.f1(r) / self.length(i,j) * self.e(i,j) * self.V[j]
146+
147+
148+
def residual(self,iter):
149+
self.f.fill(0)
150+
self.f += self.b * (iter)+ self.b * 7 + 13 * self.b3 + 1 * self.b2
151+
152+
for i in range(0,len(self.nodes)):
153+
for j in self.neighbors[i]:
154+
if not i in self.fix:
155+
tmp = self.L(i,j)
156+
self.f[2*i] += tmp[0]
157+
self.f[2*i+1] += tmp[1]
158+
159+
self.damage(i,j)
160+
161+
def damage(self,i,j):
162+
z = self.d[i]
163+
sr = 0
164+
sij = self.S(i,j)
165+
rij = self.length(i,j)
166+
167+
if rij > 1e-12:
168+
sr = np.abs(sij) / (self.rbar / np.sqrt(rij) )
169+
if z < sr :
170+
z = sr
171+
172+
self.d[i] = z
173+
174+
175+
def length(self,i,j):
176+
return np.sqrt((self.nodes[j][0]-self.nodes[i][0]) * (self.nodes[j][0]-self.nodes[i][0]) + (self.nodes[j][1]-self.nodes[i][1]) * (self.nodes[j][1]-self.nodes[i][1]) )
177+
178+
def S(self,i,j):
179+
return np.dot((self.uCurrent[j] - self.uCurrent[i]) / self.length(i,j), self.e(i,j))
180+
181+
def w(self,r):
182+
return 1
183+
184+
def e(self,i,j):
185+
return (self.nodes[j]-self.nodes[i]) / self.length(i,j)
186+
187+
def f1(self,r):
188+
return 2*r*self.C*self.beta* np.exp(-self.beta * r * r )
189+
190+
def f2(self,r):
191+
return 2*self.C*self.beta*np.exp(-r*r*self.beta)-4*self.C*r*r*self.beta*self.beta*np.exp(-r*r*self.beta)
192+
193+
def A(self,i,j):
194+
r = np.sqrt(self.length(i,j)) * self.S(i,j)
195+
return (2./self.VB) * self.w(self.length(i,j))/self.delta * self.f2(r) * ( 1. / self.length(i,j)) * self.E(i,j)
196+
197+
198+
def assemblymatrix(self):
199+
self.matrix = np.zeros((2*len(self.nodes),2*len(self.nodes)),dtype=np.double)
200+
for i in range(0,len(self.nodes)):
201+
if not i in self.fix:
202+
for j in self.neighbors[i]:
203+
tmp = self.A(i,j)
204+
205+
#Set the matrix entries for the neighbors
206+
self.matrix[i*2][j*2] += tmp[0,0]
207+
self.matrix[i*2][j*2+1] += tmp[0,1]
208+
self.matrix[i*2+1][j*2] += tmp[1,0]
209+
self.matrix[i*2+1][j*2+1] += tmp[1,1]
210+
#set the matrix entries for the node it self
211+
self.matrix[i*2][i*2] -= tmp[0,0]
212+
self.matrix[i*2][i*2+1] -= tmp[0,1]
213+
self.matrix[i*2+1][i*2] -= tmp[1,0]
214+
self.matrix[i*2+1][i*2+1] -= tmp[1,1]
215+
216+
def E(self,i,j):
217+
return np.tensordot(self.e(i,j),self.e(j,i),axes=0)
218+
219+
def intersect(self,A,B,C,D):
220+
line = LineString([(A[0], A[1]), (B[0], B[1])])
221+
other = LineString([(C[0], C[1]), (D[0], D[1])])
222+
return line.intersects(other)
223+
224+
def solve(self,maxIt,epsilion):
225+
226+
for iter in range(1,maxIt+1):
227+
228+
print(" ##### Load step: " + str(iter+self.iter) + " #####")
229+
self.residual(iter)
230+
print("Residual with intial guess",np.linalg.norm(self.f))
231+
232+
233+
residual = np.finfo(np.float).max
234+
235+
it = 1
236+
while(residual > epsilion):
237+
238+
self.assemblymatrix()
239+
240+
b = np.copy(self.f)
241+
242+
for i in range(0,len(self.fix)):
243+
244+
index = 2* self.fix[len(self.fix)-1-i]
245+
b = np.delete(b,index+1)
246+
b = np.delete(b,index)
247+
self.matrix = np.delete(self.matrix,index+1,1)
248+
self.matrix = np.delete(self.matrix,index+1,0)
249+
self.matrix = np.delete(self.matrix,index,1)
250+
self.matrix = np.delete(self.matrix,index,0)
251+
252+
253+
#print("Con:" + str(np.linalg.cond(self.matrix))+" Det: "+str(np.linalg.det(self.matrix)))
254+
255+
#res = linalg.solve(self.matrix,b)
256+
257+
inv = linalg.inv(self.matrix)
258+
259+
res = inv.dot(b)
260+
261+
262+
unew = np.zeros(2*len(self.nodes)).reshape((len(self.nodes),2))
263+
j = 0
264+
for i in range(0,len(self.uCurrent)):
265+
if not i in self.fix:
266+
unew[i] = np.array([res[2*j],res[2*j+1]])
267+
j += 1
268+
269+
270+
271+
self.uCurrent += unew
272+
273+
274+
275+
self.residual(iter)
276+
residual = np.linalg.norm(self.f)
277+
print("Iteration ",it," Residual: ",residual)
278+
it += 1
279+
280+
self.plot(iter)
281+
self.dump(iter)
282+
283+
def plot(self,iter):
284+
step = iter + self.iter
285+
# Plot u_x
286+
plt.scatter(self.nodes[:,0]+self.uCurrent[:,0],self.nodes[:,1]+self.uCurrent[:,1],c=abs(self.uCurrent[:,0]))
287+
ax = plt.gca()
288+
ax.set_facecolor('#F0F8FF')
289+
v = np.linspace(min(abs(self.uCurrent[:,0])), max(abs(self.uCurrent[:,0])), 10, endpoint=True)
290+
clb = plt.colorbar(ticks=v)
291+
clb.set_label(r'Displacement $ u_x $',labelpad=5)
292+
plt.xlabel("Position $x$")
293+
plt.ylabel("Position $y$")
294+
plt.savefig("bond-based-2d-plate-u-x-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+".pdf",bbox_inches='tight')
295+
plt.clf()
296+
# Plot u_y
297+
plt.scatter(self.nodes[:,0]+self.uCurrent[:,0],self.nodes[:,1]+self.uCurrent[:,1],c=self.uCurrent[:,1])
298+
ax = plt.gca()
299+
ax.set_facecolor('#F0F8FF')
300+
v = np.linspace(min(self.uCurrent[:,1]), max(self.uCurrent[:,1]), 10, endpoint=True)
301+
clb = plt.colorbar(ticks=v)
302+
clb.set_label(r'Displacement $ u_y $',labelpad=5)
303+
plt.xlabel("Position $x$")
304+
plt.ylabel("Position $y$")
305+
plt.savefig("bond-based-2d-plate-u-y-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+".pdf",bbox_inches='tight')
306+
plt.clf()
307+
# Plot damage
308+
plt.scatter(self.nodes[:,0]+self.uCurrent[:,0],self.nodes[:,1]+self.uCurrent[:,1],c=self.d)
309+
ax = plt.gca()
310+
ax.set_facecolor('#F0F8FF')
311+
v = np.linspace(min(self.d), max(self.d), 10, endpoint=True)
312+
clb = plt.colorbar(ticks=v)
313+
clb.set_label(r'Damage',labelpad=5)
314+
plt.xlabel("Position $x$")
315+
plt.ylabel("Position $y$")
316+
plt.savefig("bond-based-2d-plate-d-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+".pdf",bbox_inches='tight')
317+
plt.clf()
318+
319+
def dump(self,iter):
320+
step = iter + self.iter
321+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-displacement.npy", "wb")
322+
np.save(filehandler, self.uCurrent)
323+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-damage.npy", "wb")
324+
np.save(filehandler,self.d)
325+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-b.npy", "wb")
326+
np.save(filehandler,self.b)
327+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-f.npy", "wb")
328+
np.save(filehandler,self.f)
329+
330+
331+
if __name__=="__main__":
332+
333+
c = Compute(float(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3]))
334+
c.solve(100,1e-3)
335+
#c.plot()

0 commit comments

Comments
 (0)