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

Commit 4c31c77

Browse files
diehlpkdiehlpk
committed
Working version of the hard load
1 parent 00d2935 commit 4c31c77

2 files changed

Lines changed: 252 additions & 80 deletions

File tree

‎bondbased2D-plate-hard.py‎

Lines changed: 146 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#Protoype for the 2D bond-based analytic stiffness matrix including fracture
2-
# using a pre-crack square plate
2+
# using a pre-crack square plate (hard loading)
33
#@author Patrick Diehl (patrickdiehl@lsu.edu)
44
#@date September 2020
55
import numpy as np
@@ -9,6 +9,7 @@
99
import matplotlib as mpl
1010
import sys
1111
import pickle
12+
import random
1213
from matplotlib import rc
1314
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
1415
rc('text', usetex=True)
@@ -36,43 +37,43 @@ def __init__(self,h,delta_factor,iter):
3637
self.delta = self.delta_factor*h
3738
n = int(15/self.h) + 1
3839
self.fix = []
39-
self.loadT = []
40-
self.loadB = []
40+
#self.loadT = []
41+
#self.loadB = []
4142
self.z = []
4243
self.iter = iter
4344

4445
#Generate grid
4546
index = 0
46-
4747
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)
48+
for j in range(0,n+2*self.delta_factor):
49+
self.nodes.append([i*h,(j-self.delta_factor)*h])
6250

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)
51+
#if j*h < self.delta and i * h < 13 * self.delta :
52+
# self.loadB.append(index)
53+
# #plt.scatter(i*h,j*h)
6654

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 < s 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+
#plt.scatter(i*h,(j-self.delta_factor)*h)
67+
6768
index += 1
68-
69-
69+
self.nodes = np.array(self.nodes)
70+
7071
#plt.show()
7172
#sys.exit(1)
7273

73-
self.fix = np.sort(self.fix)
74+
#self.fix = np.sort(self.fix)
7475

75-
self.nodes = np.array(self.nodes)
76+
#self.nodes = np.array(self.nodes)
7677

7778
self.V = np.empty(len(self.nodes))
7879
self.V.fill(h*h)
@@ -85,31 +86,79 @@ def __init__(self,h,delta_factor,iter):
8586

8687
self.f = np.zeros(2*len(self.nodes))
8788

89+
#self.b = np.full(2*len(self.nodes),1)
90+
8891
if self.iter == 1:
8992

9093
# Initialize
9194
self.uCurrent = np.zeros(2*len(self.nodes)).reshape((len(self.nodes),2))
95+
self.wCurrent = np.zeros(2*len(self.nodes)).reshape((len(self.nodes),2))
96+
97+
#Apply the boundary load to the extension
98+
99+
100+
#scale = 1/7
101+
102+
103+
for i in range(0,len(self.nodes)):
104+
if self.nodes[i][1] < 0 and self.nodes[i][0] < 13 * self.delta :
105+
self.wCurrent[i][1] = .25
106+
107+
if self.nodes[i][1] > 15 and self.nodes[i][0] < 13 * self.delta :
108+
self.wCurrent[i][1] = -.25
109+
110+
if self.nodes[i][1] > 15 and self.nodes[i][0] >= 15 * h - self.delta :
111+
self.fix.append(i)
112+
113+
if self.nodes[i][1] < 0 and self.nodes[i][0] >= 15* h - self.delta :
114+
self.fix.append(i)
115+
116+
self.fix = np.sort(self.fix)
117+
118+
119+
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.wCurrent[:,1])
120+
#plt.colorbar()
121+
#plt.show()
122+
#sys.exit(1)
92123

93124
else:
94-
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(self.iter-1)+"-displacement.npy", "rb")
125+
print("Restart not implemented yet")
126+
#filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(self.iter-1)+"-displacement.npy", "rb")
95127
#self.nodes += np.load(filehandler)
96-
self.uCurrent = np.load(filehandler)
128+
#self.uCurrent = np.load(filehandler)
97129

98130

99131

100132
#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)
133+
#self.b = np.zeros(2*len(self.nodes))
134+
#self.b2 = np.zeros(2*len(self.nodes))
135+
#self.b3 = np.zeros(2*len(self.nodes))
136+
137+
#self.b = np.zeros(2*len(self.nodes))
138+
139+
#self.uCurrent += self.wCurrent
140+
141+
#for i in range(0,len(self.nodes)):
142+
# for j in self.neighbors[i]:
143+
#
144+
# tmp = self.L(i,j)
145+
# self.b[2*i] += tmp[0]
146+
# self.b[2*i+1] += tmp[1]
147+
148+
#self.uCurrent -= self.wCurrent
149+
150+
#print(self.b)
151+
152+
#sys.exit(1)
153+
#for i in range(0,len(self.nodes)):
154+
# if i in self.loadT:
155+
# self.b[2*i+1] = 4e3 / (2*self.delta*self.delta)
156+
# self.b2[2*i+1] = 4e5 / (2*self.delta*self.delta)
157+
# self.b3[2*i+1] = 4e6 / (2*self.delta*self.delta)
158+
# if i in self.loadB:
159+
# self.b[2*i+1] = -(4e3) / (2*self.delta*self.delta)
160+
# self.b2[2*i+1] = -(4e5) / (2*self.delta*self.delta)
161+
# self.b3[2*i+1] = -(4e6) / (2*self.delta*self.delta)
113162

114163
print("Matrix size "+str(2*len(self.nodes))+"x"+str(2*len(self.nodes)))
115164

@@ -147,17 +196,32 @@ def L(self,i,j):
147196

148197
def residual(self,iter):
149198
self.f.fill(0)
150-
self.f += self.b * (iter)+ self.b * 7 + 13 * self.b3 + 1 * self.b2
199+
200+
self.uCurrent += iter * self.wCurrent
201+
202+
for i in range(0,len(self.nodes)):
203+
for j in self.neighbors[i]:
204+
205+
tmp = self.L(i,j)
206+
self.f[2*i] += tmp[0]
207+
self.f[2*i+1] += tmp[1]
208+
209+
self.uCurrent -= iter * self.wCurrent
210+
211+
212+
def computeDamage(self):
151213

152214
for i in range(0,len(self.nodes)):
153215
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]
216+
217+
tmp = self.L(i,j)
218+
self.f[2*i] += tmp[0]
219+
self.f[2*i+1] += tmp[1]
158220

159221
self.damage(i,j)
160222

223+
224+
161225
def damage(self,i,j):
162226
z = self.d[i]
163227
sr = 0
@@ -198,20 +262,20 @@ def A(self,i,j):
198262
def assemblymatrix(self):
199263
self.matrix = np.zeros((2*len(self.nodes),2*len(self.nodes)),dtype=np.double)
200264
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)
265+
266+
for j in self.neighbors[i]:
267+
tmp = self.A(i,j)
204268

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]
269+
#Set the matrix entries for the neighbors
270+
self.matrix[i*2][j*2] += tmp[0,0]
271+
self.matrix[i*2][j*2+1] += tmp[0,1]
272+
self.matrix[i*2+1][j*2] += tmp[1,0]
273+
self.matrix[i*2+1][j*2+1] += tmp[1,1]
274+
#set the matrix entries for the node it self
275+
self.matrix[i*2][i*2] -= tmp[0,0]
276+
self.matrix[i*2][i*2+1] -= tmp[0,1]
277+
self.matrix[i*2+1][i*2] -= tmp[1,0]
278+
self.matrix[i*2+1][i*2+1] -= tmp[1,1]
215279

216280
def E(self,i,j):
217281
return np.tensordot(self.e(i,j),self.e(j,i),axes=0)
@@ -225,9 +289,10 @@ def solve(self,maxIt,epsilion):
225289

226290
for iter in range(1,maxIt+1):
227291

292+
228293
print(" ##### Load step: " + str(iter+self.iter) + " #####")
229294
self.residual(iter)
230-
print("Residual with intial guess",np.linalg.norm(self.f))
295+
print("Residual with intial guess",np.linalg.norm(self.wCurrent)-np.linalg.norm(self.uCurrent))
231296

232297

233298
residual = np.finfo(np.float).max
@@ -237,10 +302,11 @@ def solve(self,maxIt,epsilion):
237302

238303
self.assemblymatrix()
239304

305+
240306
b = np.copy(self.f)
241-
307+
242308
for i in range(0,len(self.fix)):
243-
309+
244310
index = 2* self.fix[len(self.fix)-1-i]
245311
b = np.delete(b,index+1)
246312
b = np.delete(b,index)
@@ -249,15 +315,9 @@ def solve(self,maxIt,epsilion):
249315
self.matrix = np.delete(self.matrix,index,1)
250316
self.matrix = np.delete(self.matrix,index,0)
251317

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-
257318
inv = linalg.inv(self.matrix)
258319

259-
res = inv.dot(b)
260-
320+
res = inv.dot(-b)
261321

262322
unew = np.zeros(2*len(self.nodes)).reshape((len(self.nodes),2))
263323
j = 0
@@ -266,17 +326,26 @@ def solve(self,maxIt,epsilion):
266326
unew[i] = np.array([res[2*j],res[2*j+1]])
267327
j += 1
268328

329+
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.uCurrent[:,1])
330+
#plt.colorbar()
331+
#plt.xlabel("X")
332+
#plt.xlabel("Y")
333+
#plt.title("u current")
334+
#plt.show()
269335

270336

271337
self.uCurrent += unew
272338

273-
274-
275339
self.residual(iter)
276-
residual = np.linalg.norm(self.f)
340+
341+
residual = np.linalg.norm(iter * self.wCurrent) - np.linalg.norm(self.uCurrent)
277342
print("Iteration ",it," Residual: ",residual)
278343
it += 1
279344

345+
#if it == 10 :
346+
# break
347+
348+
self.computeDamage()
280349
self.plot(iter)
281350
self.dump(iter)
282351

@@ -291,7 +360,7 @@ def plot(self,iter):
291360
clb.set_label(r'Displacement $ u_x $',labelpad=5)
292361
plt.xlabel("Position $x$")
293362
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')
363+
plt.savefig("bond-based-2d-plate-u-x-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-hard.pdf",bbox_inches='tight')
295364
plt.clf()
296365
# Plot u_y
297366
plt.scatter(self.nodes[:,0]+self.uCurrent[:,0],self.nodes[:,1]+self.uCurrent[:,1],c=self.uCurrent[:,1])
@@ -302,7 +371,7 @@ def plot(self,iter):
302371
clb.set_label(r'Displacement $ u_y $',labelpad=5)
303372
plt.xlabel("Position $x$")
304373
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')
374+
plt.savefig("bond-based-2d-plate-u-y-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-hard.pdf",bbox_inches='tight')
306375
plt.clf()
307376
# Plot damage
308377
plt.scatter(self.nodes[:,0]+self.uCurrent[:,0],self.nodes[:,1]+self.uCurrent[:,1],c=self.d)
@@ -313,23 +382,20 @@ def plot(self,iter):
313382
clb.set_label(r'Damage',labelpad=5)
314383
plt.xlabel("Position $x$")
315384
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')
385+
plt.savefig("bond-based-2d-plate-d-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-hard.pdf",bbox_inches='tight')
317386
plt.clf()
318387

319388
def dump(self,iter):
320389
step = iter + self.iter
321-
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-displacement.npy", "wb")
390+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-displacement-hard.npy", "wb")
322391
np.save(filehandler, self.uCurrent)
323-
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-damage.npy", "wb")
392+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-damage-hard.npy", "wb")
324393
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")
394+
filehandler = open("bond-based-2d-plate-"+str(self.h)+"-"+str(self.delta_factor)+"-"+str(step)+"-f-hard.npy", "wb")
328395
np.save(filehandler,self.f)
329396

330397

331398
if __name__=="__main__":
332399

333400
c = Compute(float(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3]))
334-
c.solve(100,1e-3)
335-
#c.plot()
401+
c.solve(1,1e-5)

0 commit comments

Comments
 (0)