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 (soft loading)
33#@author Patrick Diehl (patrickdiehl@lsu.edu)
44#@date September 2020
55import numpy as np
@@ -50,26 +50,22 @@ def __init__(self,h,delta_factor,iter):
5050
5151 if j * h < self .delta and i * h < 13 * self .delta :
5252 self .loadB .append (index )
53- #plt.scatter(i*h,j*h)
53+
5454
5555 if j * h > 15 + h - self .delta and i * h < 13 * self .delta :
5656 self .loadT .append (index )
57- #plt.scatter(i*h,j*h)
57+
5858
5959 if j * h < self .delta and i * h > 15 + h - self .delta :
6060 self .fix .append (index )
61- #plt.scatter(i*h,j*h)
61+
6262
6363 if j * h > 15 + h - self .delta and i * h > 15 + h - self .delta :
6464 self .fix .append (index )
65- #plt.scatter(i*h,j*h)
65+
6666
6767 index += 1
6868
69-
70- #plt.show()
71- #sys.exit(1)
72-
7369 self .fix = np .sort (self .fix )
7470
7571 self .nodes = np .array (self .nodes )
@@ -95,8 +91,6 @@ def __init__(self,h,delta_factor,iter):
9591 #self.nodes += np.load(filehandler)
9692 self .uCurrent = np .load (filehandler )
9793
98-
99-
10094 #Apply the load to the body force vector
10195 self .b = np .zeros (2 * len (self .nodes ))
10296 self .b2 = np .zeros (2 * len (self .nodes ))
@@ -119,26 +113,14 @@ def searchNeighbors(self):
119113 right = np .array ([7.5 ,7.5 ])
120114 neighbor = []
121115
122-
123- #fig = plt.gcf()
124- #fig.set_size_inches(4,50)
125-
126116 for i in range (0 ,len (self .nodes )):
127117 self .neighbors .append ([])
128118 for j in range (0 ,len (self .nodes )):
129119 if i != j and self .length (j ,i ) <= self .delta :
130120 if not self .intersect (left ,right ,self .nodes [i ],self .nodes [j ]):
131121 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)
133122 neighbor .append (len (self .neighbors [i ]))
134123
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()
142124
143125 def L (self ,i ,j ):
144126 r = np .sqrt (self .length (i ,j )) * self .S (i ,j )
@@ -249,16 +231,10 @@ def solve(self,maxIt,epsilion):
249231 self .matrix = np .delete (self .matrix ,index ,1 )
250232 self .matrix = np .delete (self .matrix ,index ,0 )
251233
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-
257234 inv = linalg .inv (self .matrix )
258235
259236 res = inv .dot (b )
260237
261-
262238 unew = np .zeros (2 * len (self .nodes )).reshape ((len (self .nodes ),2 ))
263239 j = 0
264240 for i in range (0 ,len (self .uCurrent )):
@@ -332,4 +308,3 @@ def dump(self,iter):
332308
333309 c = Compute (float (sys .argv [1 ]),int (sys .argv [2 ]),int (sys .argv [3 ]))
334310 c .solve (100 ,1e-3 )
335- #c.plot()
0 commit comments