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

Commit 42ad1cd

Browse files
diehlpkdiehlpk
committed
Clean the code
1 parent 93ba3ef commit 42ad1cd

1 file changed

Lines changed: 15 additions & 125 deletions

File tree

‎hard/bondbased2D-plate-hard.py‎

Lines changed: 15 additions & 125 deletions
Original file line numberDiff line numberDiff line change
@@ -46,24 +46,17 @@ def __init__(self,h,delta_factor,iter):
4646
#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])
50-
#index += 1
49+
self.nodes.append([i*h,(j-self.delta_factor)*h])
5150

5251
self.nodes = np.array(self.nodes)
5352

54-
#plt.show()
55-
#sys.exit(1)
56-
57-
#self.fix = np.sort(self.fix)
58-
59-
#self.nodes = np.array(self.nodes)
60-
6153
self.V = np.empty(len(self.nodes))
6254
self.V.fill(h*h)
6355

6456
self.d = np.zeros(len(self.nodes))
6557

6658
self.VB = np.pi * self.delta * self.delta
59+
6760
# Search neighbors
6861
self.searchNeighbors()
6962

@@ -80,25 +73,8 @@ def __init__(self,h,delta_factor,iter):
8073

8174
#Apply the boundary load to the extension
8275

83-
#print(min(self.nodes[:,1]))
84-
#scale = 1/7
85-
8676
for i in range(0,len(self.nodes)):
87-
#if self.nodes[i][1] <= 0 : #and self.nodes[i][0] < 2.4 :
88-
# self.wCurrentExtension[i][1] = (.125 / 11.5) * (11.5 -self.nodes[i][1])
89-
# self.loadB.append(i)
90-
#
91-
#elif self.nodes[i][1] >= 15 : #and self.nodes[i][0] < 2.4 :
92-
# self.wCurrentExtension[i][1] = 0 #-0.125
93-
# self.loadT.append(i)
94-
#
95-
#elif self.nodes[i][1] > 0 and self.nodes[i][1] <= 2.5 :
96-
# self.wCurrentDomain[i][1] = (.125 / 11.5) * (11.5 -self.nodes[i][1])
97-
#
98-
#elif self.nodes[i][1] < 15 and self.nodes[i][1] >= 2.5 :
99-
# self.wCurrentDomain[i][1] = 0 # (-.125 / 7.5) * (self.nodes[i][1]-7.5)
100-
101-
77+
10278
if self.nodes[i][1] > 7.5 :
10379
self.wCurrentExtension[i][1] = (-1e-1 ) * (self.nodes[i][1] - 7.5 ) / 15.5 - 1e-2 # * (15 - self.nodes[i][0] ) / 15
10480
self.loadB.append(i)
@@ -107,23 +83,6 @@ def __init__(self,h,delta_factor,iter):
10783
self.wCurrentExtension[i][1] = (1e-1 ) * ( -7.5 + self.nodes[i][1] ) / -15.5 + 1e-2 # * (15 - self.nodes[i][0] ) / 15
10884
self.loadT.append(i)
10985

110-
#if self.nodes[i][1] >= 15 :
111-
# self.wCurrentExtension[i][1] = -1e-4
112-
113-
#if self.nodes[i][1] <= 0 :
114-
# self.wCurrentExtension[i][1] = 1e-4
115-
116-
117-
#if self.nodes[i][0] > 15 - h and self.nodes[i][1] >= 7.5 and self.nodes[i][1] < 7.5 + 1*h :
118-
# self.fix.append(i)
119-
# self.wCurrentExtension[i][1] = 0
120-
# print("d")
121-
122-
#elif self.nodes[i][0] > 15-h and self.nodes[i][1] < 7.5 and self.nodes[i][1] >= 7.5 - 1*h :
123-
# self.fix.append(i)
124-
# self.wCurrentExtension[i][1] = 0
125-
# print("d2")
126-
12786
if self.nodes[i][0] > 15 - h and self.nodes[i][1] > 15 + self.delta - h :
12887
self.fix.append(i)
12988
self.wCurrentExtension[i][1] = 0
@@ -133,29 +92,19 @@ def __init__(self,h,delta_factor,iter):
13392
self.fix.append(i)
13493
self.wCurrentExtension[i][1] = 0
13594
print("d2")
136-
137-
138-
#elif self.nodes[i][1] >= 15 and self.nodes[i][0] > 2.4 and self.nodes[i][0] <= 15 :
139-
# self.wCurrent[i][1] = -.125 / 5
140-
141-
#elif self.nodes[i][1] <= 0 and self.nodes[i][0] > 2.4 and self.nodes[i][0] <= 15 :
142-
# self.wCurrent[i][1] = .125 / 5
143-
144-
#else:
145-
#plt.scatter(self.nodes[i][0],self.nodes[i][1],color="red")
14695

14796
self.fix = np.sort(self.fix)
14897
print(self.fix)
14998

150-
#self.remove = np.concatenate([self.fix,self.loadT,self.loadB])
151-
#self.remove = np.sort(self.remove)
152-
153-
154-
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c= self.wCurrentExtension[:,1],cmap="seismic")
155-
#for i in self.fix :
156-
# plt.scatter(self.nodes[i,0],self.nodes[i,1],c="black")
157-
#plt.colorbar()
158-
#plt.show()
99+
#ax = plt.gca()
100+
#ax.set_facecolor('#F0F8FF')
101+
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=-self.wCurrentExtension[:,1],cmap="seismic",marker="s",s=np.sqrt(30))
102+
#v = np.linspace(min(self.wCurrentExtension[:,1]), max(self.wCurrentExtension[:,1]), 6, endpoint=True)
103+
#clb = plt.colorbar(ticks=v,format='%.1e')
104+
#clb.set_label(r'Displacement $w_y$')
105+
#plt.xlabel(r"Position $x$")
106+
#plt.ylabel(r"Poistion $y$")
107+
#plt.savefig("w.pdf",bbox_inches='tight')
159108
#sys.exit(1)
160109

161110
else:
@@ -172,26 +121,12 @@ def searchNeighbors(self):
172121
right = np.array([7.5,7.5])
173122
neighbor = []
174123

175-
176-
#fig = plt.gcf()
177-
#fig.set_size_inches(4,50)
178-
179124
for i in range(0,len(self.nodes)):
180125
self.neighbors.append([])
181126
for j in range(0,len(self.nodes)):
182127
if i != j and self.length(j,i) <= self.delta:
183128
if not self.intersect(left,right,self.nodes[i],self.nodes[j]):
184129
self.neighbors[i].append(j)
185-
#plt.plot([self.nodes[i][0],self.nodes[j][0]],[self.nodes[i][1],self.nodes[j][1]],c="#91A3B0",alpha=0.15)
186-
#neighbor.append(len(self.neighbors[i]))
187-
188-
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=neighbor)
189-
#plt.xlabel("Position $x$",fontsize = 30)
190-
#plt.ylabel("Position $y$",fontsize = 30)
191-
#clb = plt.colorbar()
192-
#clb.set_label(r'$B_\delta(x)$',labelpad=5)
193-
#plt.show()
194-
#sys.exit()
195130

196131
def L(self,i,j):
197132
r = np.sqrt(self.length(i,j)) * self.S(i,j)
@@ -204,9 +139,6 @@ def computeLoad(self,iter):
204139

205140
self.uCurrent += iter * self.wCurrentExtension
206141

207-
#if iter == 1:
208-
#self.uCurrent += iter * self.wCurrentDomain
209-
210142
for i in range(0,len(self.nodes)):
211143
for j in self.neighbors[i]:
212144

@@ -216,24 +148,16 @@ def computeLoad(self,iter):
216148

217149
self.uCurrent -= iter * self.wCurrentExtension
218150

219-
#if iter == 1:
220-
#self.uCurrent -= iter * self.wCurrentDomain
221-
222151
def residual(self,iter):
223152
self.f.fill(0)
224153

225-
#self.uCurrent =+ iter * self.wCurrent
226-
227154
for i in range(0,len(self.nodes)):
228155
for j in self.neighbors[i]:
229156

230157
tmp = self.L(i,j)
231158
self.f[2*i] += tmp[0]
232159
self.f[2*i+1] += tmp[1]
233160

234-
#self.uCurrent =- iter * self.wCurrent
235-
236-
237161
def computeDamage(self):
238162

239163
for i in range(0,len(self.nodes)):
@@ -316,11 +240,7 @@ def solve(self,maxIt,epsilion):
316240

317241

318242
print(" ##### Load step: " + str(iter+self.iter) + " #####")
319-
#self.residual(iter)
320-
321-
#print("Residual with intial guess",np.linalg.norm(self.f))
322-
323-
243+
324244
residual = np.finfo(float).max
325245
residual_old = 0
326246

@@ -337,6 +257,7 @@ def solve(self,maxIt,epsilion):
337257
self.matrix = np.delete(self.matrix,index,0)
338258

339259
it = 1
260+
340261
while(residual > epsilion):
341262

342263

@@ -374,47 +295,16 @@ def solve(self,maxIt,epsilion):
374295
if not i in self.fix:
375296
unew[i] = np.array([res[2*j],res[2*j+1]])
376297
j += 1
377-
#elif i in self.loadT:
378-
# unew[i] = iter* np.array([0,0.125])
379-
#elif i in self.loadB:
380-
# unew[i] = iter * np.array([0,-.125])
381-
382-
#plt.scatter(self.nodes[:,0],self.nodes[:,1],c=self.uCurrent[:,1])
383-
#plt.colorbar()
384-
#plt.xlabel("X")
385-
#plt.xlabel("Y")
386-
#plt.title("u current")
387-
#plt.show()
388-
389298

390299
self.uCurrent += unew
391300

392-
#self.residual(iter)
393301
self.computeLoad(iter)
394-
#self.residual(iter)
395-
396-
#extension1 = []
397-
#extension2 = []
398-
399-
#for i in range(0,len(self.nodes)):
400-
# if i in self.loadT or i in self.loadB :
401-
# #extension.append(abs(iter*self.wCurrentExtension[i])-abs(self.uCurrent[i]))
402-
# extension1.append(iter*abs(self.wCurrentExtension[i]))
403-
# extension2.append(self.uCurrent[i])
404-
405-
#print(np.linalg.norm(extension1),np.linalg.norm(extension2))
406-
#residual = np.linalg.norm(extension1) - np.linalg.norm(extension2)
407-
#residual = np.linalg.norm(iter * self.wCurrent-extension)
408-
#residual = np.linalg.norm(extension)
409-
#residual = np.linalg.norm(self.b+self.f)
302+
410303
residual = np.linalg.norm(self.uCurrent) - residual_old
411304
residual_old = np.linalg.norm(self.uCurrent)
412305
print("Iteration ",it," Residual: ",residual)
413306
it += 1
414307

415-
#if it == 10 :
416-
# break
417-
418308
self.computeDamage()
419309
self.plot(iter)
420310
self.dump(iter)

0 commit comments

Comments
 (0)