Note
Go to the end to download the full example code.
PostProcess#
Script to compare damage simulations.
Missing Simulations :
Amor_AT1_DP_optimMesh
Amor_AT1_DP_optimMesh_conv0.1
Amor_AT1_DP_optimMesh_conv0.01
Amor_AT1_DP_optimMesh_conv0.001
Amor_AT1_DP_optimMesh_conv0.0001
Amor_AT1_DP_optimMesh_conv1e-05
12 from EasyFEA import Display, Folder, Tic, plt, np, pd, Simulations
13
14 if __name__ == "__main__":
15 Display.Clear()
16
17 # ----------------------------------------------
18 # Configuration
19 # ----------------------------------------------
20 # "PlateWithHole_Benchmark", "Shear_Benchmark", "Tension_Benchmark" "L_Shape_Benchmark"
21 simulation = "Shear_Benchmark"
22 # simulation = "PlateWithHole_Benchmark"
23
24 meshTest = False
25 loadSimu = True
26 plotDamage = False
27 savefig = True
28
29 folder_results = Folder.Join(Folder.RESULTS_DIR, simulation, mkdir=True)
30
31 if savefig:
32 if meshTest:
33 folder_save = Folder.Join(folder_results, "Test", "_Post processing")
34 else:
35 folder_save = Folder.Join(folder_results, "_Post processing")
36 else:
37 folder_save = ""
38
39 list_mat = ["Elas_Isot"] # ["Elas_Isot", "Elas_IsotTrans", "ElasAnisot"]
40
41 # list_regu = ["AT2", "AT1"] # ["AT1", "AT2"]
42 list_regu = ["AT1"] # ["AT1", "AT2"]
43
44 list_simpli2D = ["DP"] # ["CP","DP"]
45 list_solver = ["History"]
46
47 # list_split = ["Bourdin","Amor","Miehe","He","Zhang"]
48 # list_split = ["Bourdin","Amor","Miehe","He","AnisotStrain","AnisotStress","Zhang"]
49 # list_split = ["Bourdin","He","AnisotStrain","AnisotStress","Zhang"]
50 # list_split = ["AnisotStrain","AnisotStress", "He", "Zhang"]
51 # list_split = ["AnisotStress", "He"]
52 # list_split = ["Bourdin","Amor","Miehe"]
53 list_split = ["Amor"]
54
55 # listOptimMesh=[False, True] # [True, False]
56 listOptimMesh = [True] # [True, False]
57
58 listTol = [1e-0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5] # [1e-0, 1e-1, 1e-2, 1e-3, 1e-4]
59 # listTol = [1e-0, 1e-2]
60 # listTol = [1e-0]
61
62 # listnL = [100] # [100] [100, 120, 140, 180, 200]
63 listnL = [0]
64
65 listTheta = [0]
66 # listTheta = [-0, -10, -20, -30, -45, -60, -70, -80, -90]
67
68 # snapshots = [18.5, 24.6, 25, 28, 35]
69 # snapshots = [9.05, 10.5, 14.5]
70 # snapshots = [18.6,24.6,24.8]
71 # snapshots = [24.6, 25]
72 snapshots = []
73
74 # depMax = 20e-5
75 # depMax = 2.5e-5
76 # depMax = 2.46e-5
77 depMax = 3.5e-5
78 # depMax = 2.5e-5
79
80 # Génération des configurations
81 listConfig = []
82
83 for theta in listTheta:
84 for comp in list_mat:
85 for split in list_split:
86 for regu in list_regu:
87 for simpli2D in list_simpli2D:
88 for solveur in list_solver:
89 for tol in listTol:
90 for optimMesh in listOptimMesh:
91 for nL in listnL:
92 listConfig.append(
93 [
94 comp,
95 regu,
96 simpli2D,
97 solveur,
98 split,
99 tol,
100 optimMesh,
101 nL,
102 theta,
103 ]
104 )
105
106 Nconfig = len(listConfig)
107
108 # ----------------------------------------------
109 # Loads all simulations
110 # ----------------------------------------------
111 ax_load = Display.Init_Axes() # superposition axis of force-displacement curves
112
113 missingSimulations = []
114
115 for config in listConfig:
116 comp = config[0]
117 regu = config[1]
118 simpli2D = config[2]
119 solveur = config[3]
120 split = config[4]
121 tolConv = config[5]
122 optimMesh = config[6]
123 nL = config[7]
124 theta = config[8]
125
126 tic = Tic()
127
128 foldername = Folder.PhaseField_Folder(
129 folder_results,
130 material=comp,
131 split=split,
132 regu=regu,
133 simpli2D=simpli2D,
134 tolConv=tolConv,
135 solver=solveur,
136 test=meshTest,
137 optimMesh=optimMesh,
138 closeCrack=False,
139 nL=nL,
140 theta=theta,
141 )
142
143 fileForceDep = Folder.Join(foldername, "force-displacement.pickle")
144 fileSimu = Folder.Join(foldername, "simulation.pickle")
145
146 nomSimu = foldername.split(comp + "_")[-1]
147
148 # text = nomSimu
149 # text = split
150 # text = f"{split}_{regu}_tol{tolConv:1.0e}"
151 text = f"{tolConv:1.0e}"
152 # text = f"{split}_{regu}"
153 # text = foldername.replace(Folder.Dir(foldername), "")[1:]
154
155 # text = f"{split} {regu}"
156 # if optimMesh:
157 # text += f" optim"
158 # else:
159 # text += f" unif"
160
161 if (loadSimu or plotDamage) and Folder.Exists(fileSimu):
162 # Load simulation
163 simu = Simulations.Load_Simu(foldername)
164 simu.mesh.groupElem.coord
165 results = pd.DataFrame(simu.results)
166 temps = results["timeIter"].values.sum()
167 temps_str, unite = Tic.Get_time_unity(temps)
168 print(len(results), f"-> {temps_str:.3} {unite}")
169
170 else:
171 if nomSimu not in missingSimulations:
172 missingSimulations.append(nomSimu)
173
174 # ----------------------------------------------
175 # Plot loads
176 # ----------------------------------------------
177
178 pltCrack = True
179
180 # Loads force and displacement
181 if Folder.Exists(fileForceDep):
182 force, displacement = Simulations.Load_Force_Displacement(foldername)
183
184 if pltCrack:
185 damage = np.asarray(
186 [
187 simu.Result("damage", iter=i).max()
188 for i in range(len(simu.results))
189 ]
190 )
191 i_crack = np.where(damage >= 1 - 1e-12)[0][0]
192
193 fc = np.abs(force[i_crack] * 1e-3)
194 print(f"fc = {fc:.2f} N/mm")
195 uc = displacement[i_crack] * 1e6
196 print(f"uc = {uc:.2f} µm")
197 # print(f"{displacement[-1]*1e6:.2f} µm")
198
199 if depMax == 0:
200 depMax = displacement[-1]
201
202 indexLim = np.where(displacement <= depMax)[0]
203
204 # text += f" ({temps_str:.3} {unite})"
205
206 ls = None
207 c = None
208
209 # ls = '--' if optimMesh else None
210 # c = ax_load.get_lines()[-1].get_color() if optimMesh else None
211
212 # ls = '--' if regu == "AT1" else None
213 # c = ax_load.get_lines()[-1].get_color() if regu == "AT1" else None
214
215 ax_load.plot(
216 displacement[indexLim] * 1e6,
217 np.abs(force[indexLim]) * 1e-6,
218 c=c,
219 label=text,
220 ls=ls,
221 )
222 if pltCrack:
223 ax_load.scatter(
224 displacement[i_crack] * 1e6,
225 np.abs(force[i_crack]) * 1e-6,
226 c=c,
227 marker="+",
228 s=80,
229 label=f"{fc:.2f} N/mm",
230 )
231
232 else:
233 if nomSimu not in missingSimulations:
234 missingSimulations.append(nomSimu)
235
236 # ----------------------------------------------
237 # Plot loads
238 # ----------------------------------------------
239 if plotDamage and Folder.Exists(fileSimu):
240 # Display.Plot_Mesh(simu.mesh)
241
242 if simulation == "PlateWithHole_Benchmark":
243 colorBarIsClose = True
244 else:
245 colorBarIsClose = False
246
247 # Displays last damage
248 filename = foldername.replace(Folder.Dir(foldername), "")[1:]
249
250 ax = Display.Plot_Result(simu, "damage", ncolors=21, colorbarIsClose=True)
251 ax.axis("off")
252 ax.set_title("")
253 Display.Save_fig(folder_save, f"{filename} last")
254
255 # Display.Plot_Result(simu, "damage", ncolors=21, nodeValues=True, colorbarIsClose=colorBarIsClose,
256 # folder=folder_save, filename=f"{split} tol{tolConv} last", plotMesh=False,
257 # title=f"{split}_{regu}_tol{tolConv}")
258
259 # Recover snapshot iterations
260
261 if len(snapshots) > 0:
262 # ax_load_2 = Display.Init_Axes()
263 # ax_load_2.plot(displacement[indexLim]*1e3, np.abs(force[indexLim])*1e-6, label=text)
264 # ax_load_2.set_xlabel("Déplacement [mm]")
265 # ax_load_2.set_ylabel("Force [kN/mm]")
266
267 for d, dep in enumerate(snapshots):
268 try:
269 i = np.where(displacement * 1e6 >= dep)[0][0]
270 except IndexError:
271 continue
272
273 simu.Set_Iter(i)
274 filenameDamage = (
275 f"{nomSimu}, ud = {np.round(displacement[i] * 1e6, 2)}"
276 )
277 # titleDamage = filenameDamage
278
279 # titleDamage = split
280 # filenameDamage = f"PlateBench {comp}_{split}_{regu}_{simpli2D}"
281
282 # Display.Plot_Result(simu, "damage", nodeValues=True, colorbarIsClose=colorBarIsClose, folder=folder_save, filename=filenameDamage,title=titleDamage)
283
284 ax = Display.Plot_Result(
285 simu, "damage", ncolors=21, colorbarIsClose=True
286 )
287 ax.axis("off")
288 ax.set_title("")
289 # Display._Remove_colorbar(ax)
290 filename = foldername.replace(Folder.Dir(foldername), "")[1:]
291 Display.Save_fig(
292 folder_save, f"{filename} snap={dep}", True
293 ) # , "png", 200
294
295 # ax_load_2.scatter(displacement[i]*1e3, force[i]*1e-6, c='k', marker='+', zorder=2)
296
297 # plt.figure(ax_load_2.figure)
298 # Display.Save_fig(folder_save, f"load displacement {filename}")
299
300 # text = text.replace("AnisotStrain","Spectral")
301 tic.Tac("PostProcessing", split, False)
302
303 ax_load.set_xlabel("Déplacement [µm]")
304 ax_load.set_ylabel("Force [kN/mm]")
305 # ax_load.set_xlabel("Déplacement [mm]")
306 # ax_load.set_ylabel("Force [kN/mm]")
307
308 # ax_load.set_xlabel("displacement [µm]")
309 # ax_load.set_ylabel("load [kN/mm]")
310 # ax_load.set_xlabel("displacement")
311 # ax_load.set_ylabel("load")
312 ax_load.grid()
313
314 if pltCrack:
315 handles, labels = ax_load.get_legend_handles_labels()
316 handles = np.concatenate((handles[::2], handles[1::2]), axis=0)
317 labels = np.concatenate((labels[::2], labels[1::2]), axis=0)
318 ax_load.legend(handles, labels, ncols=2)
319 else:
320 ax_load.legend()
321
322 plt.figure(ax_load.figure)
323 Display.Save_fig(folder_save, "load displacement")
324
325 print("\nMissing Simulations :")
326 [print(simul) for simul in missingSimulations]
327
328 plt.show()
Total running time of the script: (0 minutes 0.092 seconds)