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

Commit 3da4d85

Browse files
diehlpkdiehlpk
committed
New version
1 parent 3b1b7f7 commit 3da4d85

2 files changed

Lines changed: 90 additions & 81 deletions

File tree

‎bondbased2D-plate-hard.py‎

Lines changed: 86 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -37,35 +37,18 @@ def __init__(self,h,delta_factor,iter):
3737
self.delta = self.delta_factor*h
3838
n = int(15/self.h) + 1
3939
self.fix = []
40-
#self.loadT = []
41-
#self.loadB = []
40+
self.loadT = []
41+
self.loadB = []
4242
self.z = []
4343
self.iter = iter
4444

4545
#Generate grid
46-
index = 0
46+
#index = 0
4747
for i in range(0,n):
4848
for j in range(0,n+2*self.delta_factor):
49-
self.nodes.append([i*h,(j-self.delta_factor)*h])
49+
self.nodes.append([i*h,(j-self.delta_factor)*h])
50+
#index += 1
5051

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 < 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-
68-
index += 1
6952
self.nodes = np.array(self.nodes)
7053

7154
#plt.show()
@@ -86,7 +69,7 @@ def __init__(self,h,delta_factor,iter):
8669

8770
self.f = np.zeros(2*len(self.nodes))
8871

89-
#self.b = np.full(2*len(self.nodes),1)
72+
self.b = np.zeros(2*len(self.nodes))
9073

9174
if self.iter == 0:
9275

@@ -101,22 +84,48 @@ def __init__(self,h,delta_factor,iter):
10184

10285

10386
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
87+
if self.nodes[i][1] <= 0 : # and self.nodes[i][0] < 2.4 :
88+
self.wCurrent[i][1] = 0.125
89+
self.loadB.append(i)
90+
91+
elif self.nodes[i][1] >= 15 : # and self.nodes[i][0] < 2.4 :
92+
self.wCurrent[i][1] = -0.125
93+
self.loadT.append(i)
10694

107-
if self.nodes[i][1] > 15 and self.nodes[i][0] < 13 * self.delta :
108-
self.wCurrent[i][1] = .25
95+
elif self.nodes[i][1] > 0 and self.nodes[i][1] <= 7.5 :
96+
self.wCurrent[i][1] = (.125 / 7.5) * (7.5 -self.nodes[i][1])
10997

110-
if self.nodes[i][1] > 15 and self.nodes[i][0] >= 15 * h - self.delta :
111-
self.fix.append(i)
98+
elif self.nodes[i][1] < 15 and self.nodes[i][1] >= 7.5 :
99+
self.wCurrent[i][1] = (-.125 / 7.5) * (self.nodes[i][1]-7.5)
100+
101+
102+
if self.nodes[i][0] > 15 - h and self.nodes[i][1] >= 7.5 and self.nodes[i][1] < 7.5 + 1*h :
103+
self.fix.append(i)
104+
self.wCurrent[i][1] = 0
105+
print("d")
112106

113-
if self.nodes[i][1] < 0 and self.nodes[i][0] >= 15* h - self.delta :
107+
elif self.nodes[i][0] > 15-h and self.nodes[i][1] < 7.5 and self.nodes[i][1] >= 7.5 - 1*h :
114108
self.fix.append(i)
109+
self.wCurrent[i][1] = 0
110+
print("d2")
111+
112+
#elif self.nodes[i][1] >= 15 and self.nodes[i][0] > 2.4 and self.nodes[i][0] <= 15 :
113+
# self.wCurrent[i][1] = -.125 / 5
115114

115+
#elif self.nodes[i][1] <= 0 and self.nodes[i][0] > 2.4 and self.nodes[i][0] <= 15 :
116+
# self.wCurrent[i][1] = .125 / 5
117+
118+
#else:
119+
#plt.scatter(self.nodes[i][0],self.nodes[i][1],color="red")
120+
116121
self.fix = np.sort(self.fix)
122+
print(self.fix)
123+
124+
#self.remove = np.concatenate([self.fix,self.loadT,self.loadB])
125+
#self.remove = np.sort(self.remove)
117126

118127

119-
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.wCurrent[:,1])
128+
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.wCurrent[:,1],cmap="seismic")
120129
#plt.colorbar()
121130
#plt.show()
122131
#sys.exit(1)
@@ -127,39 +136,6 @@ def __init__(self,h,delta_factor,iter):
127136
#self.nodes += np.load(filehandler)
128137
#self.uCurrent = np.load(filehandler)
129138

130-
131-
132-
#Apply the load to the body force vector
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)
162-
163139
print("Matrix size "+str(2*len(self.nodes))+"x"+str(2*len(self.nodes)))
164140

165141

@@ -193,20 +169,36 @@ def L(self,i,j):
193169
r = np.sqrt(self.length(i,j)) * self.S(i,j)
194170
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]
195171

172+
def computeLoad(self,iter):
173+
self.b.fill(0)
174+
175+
#uCurrentCopy = np.copy(self.uCurrent)
176+
177+
self.uCurrent += iter * self.wCurrent
178+
179+
for i in range(0,len(self.nodes)):
180+
for j in self.neighbors[i]:
181+
182+
tmp = self.L(i,j)
183+
self.b[2*i] += tmp[0]
184+
self.b[2*i+1] += tmp[1]
185+
186+
self.uCurrent -= iter * self.wCurrent #np.copy(uCurrentCopy)
187+
196188

197189
def residual(self,iter):
198190
self.f.fill(0)
199191

200-
self.uCurrent += iter * self.wCurrent
201-
192+
#self.uCurrent =+ iter * self.wCurrent
193+
202194
for i in range(0,len(self.nodes)):
203195
for j in self.neighbors[i]:
204196

205197
tmp = self.L(i,j)
206198
self.f[2*i] += tmp[0]
207199
self.f[2*i+1] += tmp[1]
208200

209-
self.uCurrent -= iter * self.wCurrent
201+
#self.uCurrent =- iter * self.wCurrent
210202

211203

212204
def computeDamage(self):
@@ -220,8 +212,6 @@ def computeDamage(self):
220212

221213
self.damage(i,j)
222214

223-
224-
225215
def damage(self,i,j):
226216
z = self.d[i]
227217
sr = 0
@@ -292,22 +282,25 @@ def solve(self,maxIt,epsilion):
292282

293283
print(" ##### Load step: " + str(iter+self.iter) + " #####")
294284
self.residual(iter)
295-
print("Residual with intial guess",np.linalg.norm(self.wCurrent)-np.linalg.norm(self.uCurrent))
296285

286+
print("Residual with intial guess",np.linalg.norm(self.f))
287+
288+
289+
residual = np.finfo(float).max
297290

298-
residual = np.finfo(np.float).max
291+
self.computeLoad(iter)
299292

300293
it = 1
301294
while(residual > epsilion):
302295

303296
self.assemblymatrix()
304297

305298

306-
b = np.copy(self.f)
299+
b = np.copy(self.f + self.b)
307300

308301
for i in range(0,len(self.fix)):
309-
310-
index = 2* self.fix[len(self.fix)-1-i]
302+
303+
index = int(2* self.fix[len(self.fix)-1-i])
311304
b = np.delete(b,index+1)
312305
b = np.delete(b,index)
313306
self.matrix = np.delete(self.matrix,index+1,1)
@@ -325,6 +318,10 @@ def solve(self,maxIt,epsilion):
325318
if not i in self.fix:
326319
unew[i] = np.array([res[2*j],res[2*j+1]])
327320
j += 1
321+
#elif i in self.loadT:
322+
# unew[i] = iter* np.array([0,0.125])
323+
#elif i in self.loadB:
324+
# unew[i] = iter * np.array([0,-.125])
328325

329326
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.uCurrent[:,1])
330327
#plt.colorbar()
@@ -336,9 +333,21 @@ def solve(self,maxIt,epsilion):
336333

337334
self.uCurrent += unew
338335

339-
self.residual(iter)
336+
#self.residual(iter)
337+
self.computeLoad(iter)
338+
339+
extension = []
340+
341+
342+
for i in range(0,len(self.nodes)):
343+
if i in self.loadT or i in self.loadB :
344+
#extension.append(abs(iter*self.wCurrent[i])-abs(self.uCurrent[i]))
345+
extension.append(self.b[i])
340346

341-
residual = np.linalg.norm(iter * self.wCurrent) - np.linalg.norm(self.uCurrent)
347+
#residual = np.linalg.norm(iter * self.wCurrent) - np.linalg.norm(extension)
348+
#residual = np.linalg.norm(iter * self.wCurrent-extension)
349+
residual = np.linalg.norm(extension)
350+
#residual = np.linalg.norm(b)
342351
print("Iteration ",it," Residual: ",residual)
343352
it += 1
344353

@@ -398,4 +407,4 @@ def dump(self,iter):
398407
if __name__=="__main__":
399408

400409
c = Compute(float(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3]))
401-
c.solve(1,1e-5)
410+
c.solve(500,1e-5)

‎plot-plate-hard.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ def grayscale_cmap(cmap):
5858

5959
for i in range(0,len(nodes)):
6060
if nodes[i][1] >=0 and nodes[i][1] <= 15:
61-
#if uCurrent[i][1] > 1e-4 :
61+
6262
nodes_small.append(nodes[i])
6363
u_small.append(uCurrent[i])
6464
d_small.append(damage[i])
@@ -70,7 +70,7 @@ def grayscale_cmap(cmap):
7070

7171

7272
ax = plt.gca()
73-
plt.scatter(nodes_small[:,0]+u_small[:,0],nodes_small[:,1]+u_small[:,1],c=u_small[:,1],cmap=grayscale_cmap(plt.get_cmap("viridis")),marker="s",s=np.sqrt(30))
73+
plt.scatter(nodes_small[:,0]+u_small[:,0],nodes_small[:,1]+u_small[:,1],c=u_small[:,1],cmap="seismic",marker="s",s=np.sqrt(30))
7474
#ax.xaxis.tick_top()
7575
#ax.xaxis.set_label_position('top')
7676
fig = plt.gcf()
@@ -86,8 +86,8 @@ def grayscale_cmap(cmap):
8686
plt.savefig("bond-based-2d-plate-"+str(h)+"-"+str(delta_factor)+"-"+str(iter)+"-u-y-rotated-hard.pdf",bbox_inches='tight')
8787
plt.clf()
8888

89-
#plt.figure(figsize=(20,1))
90-
plt.scatter(nodes_small[:,0],nodes_small[:,1],c=d_small,cmap=grayscale_cmap(plt.get_cmap("viridis")),marker="s",s=np.sqrt(30))
89+
#plt.figure(figsize=(20,1)) grayscale_cmap(plt.get_cmap("viridis"))
90+
plt.scatter(nodes_small[:,0],nodes_small[:,1],c=d_small,cmap="seismic",marker="s",s=np.sqrt(30))
9191
ax = plt.gca()
9292
#ax.xaxis.tick_top()
9393
#ax.xaxis.set_label_position('top')

0 commit comments

Comments
 (0)