Note
Go to the end to download the full example code.
PlateWithHole#
Damage simulation for a plate with a hole subjected to compression.
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh 2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:0.892 s, tol=0.00e+00
3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.865 s, tol=1.00e+00 3.20 % -> 52.31 s
4 : 1.600 μm, [1.03e-07; 3.70e-02], 1:0.869 s, tol=8.04e-01 6.40 % -> 38.13 s
5 : 2.400 μm, [3.74e-05; 3.70e-02], 1:0.858 s, tol=5.54e-01 9.60 % -> 32.30 s
6 : 3.200 μm, [8.41e-05; 3.70e-02], 1:0.881 s, tol=4.37e-01 12.80 % -> 30.01 s
7 : 4.000 μm, [1.50e-04; 3.70e-02], 1:0.891 s, tol=3.60e-01 16.00 % -> 28.08 s
8 : 4.800 μm, [2.34e-04; 3.70e-02], 1:0.898 s, tol=3.05e-01 19.20 % -> 26.45 s
9 : 5.600 μm, [3.37e-04; 3.70e-02], 1:0.885 s, tol=2.65e-01 22.40 % -> 24.52 s
10 : 6.400 μm, [4.60e-04; 3.71e-02], 1:0.891 s, tol=2.34e-01 25.60 % -> 23.30 s
11 : 7.200 μm, [6.01e-04; 4.57e-02], 1:0.870 s, tol=2.10e-01 28.80 % -> 21.52 s
12 : 8.000 μm, [7.63e-04; 5.81e-02], 1:0.896 s, tol=1.90e-01 32.00 % -> 20.93 s
13 : 8.800 μm, [9.46e-04; 7.21e-02], 1:0.860 s, tol=1.73e-01 35.20 % -> 19.01 s
14 : 9.600 μm, [1.15e-03; 8.78e-02], 1:0.886 s, tol=1.59e-01 38.40 % -> 18.48 s
15 : 10.400 μm, [1.38e-03; 1.05e-01], 1:0.860 s, tol=1.47e-01 41.60 % -> 16.90 s
16 : 11.200 μm, [1.62e-03; 1.25e-01], 1:0.883 s, tol=1.37e-01 44.80 % -> 16.32 s
17 : 12.000 μm, [1.89e-03; 1.46e-01], 1:0.904 s, tol=1.28e-01 48.00 % -> 15.67 s
18 : 12.800 μm, [2.19e-03; 1.70e-01], 1:0.891 s, tol=1.21e-01 51.20 % -> 14.43 s
19 : 13.600 μm, [2.50e-03; 1.97e-01], 1:0.864 s, tol=1.14e-01 54.40 % -> 13.04 s
20 : 14.400 μm, [2.85e-03; 2.29e-01], 1:0.888 s, tol=1.07e-01 57.60 % -> 12.43 s
21 : 15.200 μm, [3.21e-03; 2.65e-01], 1:0.874 s, tol=1.02e-01 60.80 % -> 11.27 s
22 : 16.000 μm, [3.61e-03; 3.04e-01], 1:0.892 s, tol=9.68e-02 64.00 % -> 10.53 s
23 : 16.800 μm, [4.03e-03; 3.48e-01], 1:0.887 s, tol=9.23e-02 67.20 % -> 9.53 s
24 : 17.600 μm, [4.47e-03; 3.96e-01], 1:0.868 s, tol=8.82e-02 70.40 % -> 8.40 s
25 : 18.400 μm, [4.95e-03; 4.49e-01], 1:0.855 s, tol=8.44e-02 73.60 % -> 7.36 s
26 : 19.200 μm, [5.45e-03; 5.05e-01], 1:0.882 s, tol=8.10e-02 76.80 % -> 6.66 s
27 : 20.000 μm, [5.98e-03; 5.61e-01], 1:0.872 s, tol=7.80e-02 80.00 % -> 5.67 s
28 : 20.800 μm, [6.53e-03; 6.38e-01], 1:0.905 s, tol=7.53e-02 83.20 % -> 4.94 s
29 : 21.000 μm, [7.11e-03; 7.36e-01], 1:0.872 s, tol=2.04e-02 84.00 % -> 4.65 s
30 : 21.200 μm, [7.30e-03; 8.35e-01], 1:0.904 s, tol=2.10e-02 84.80 % -> 4.70 s
31 : 21.400 μm, [7.46e-03; 9.20e-01], 1:0.879 s, tol=2.34e-02 85.60 % -> 4.43 s
32 : 21.600 μm, [7.64e-03; 9.75e-01], 1:0.894 s, tol=2.99e-02 86.40 % -> 4.36 s
33 : 21.800 μm, [7.87e-03; 9.98e-01], 1:0.875 s, tol=4.11e-02 87.20 % -> 4.11 s
34 : 22.000 μm, [8.18e-03; 1.00e+00], 1:0.877 s, tol=5.57e-02 88.00 % -> 3.95 s
35 : 22.200 μm, [8.57e-03; 1.00e+00], 1:0.870 s, tol=6.35e-02 88.80 % -> 3.73 s
36 : 22.400 μm, [8.97e-03; 1.00e+00], 1:0.882 s, tol=8.48e-02 89.60 % -> 3.58 s
37 : 22.600 μm, [9.40e-03; 1.00e+00], 1:0.864 s, tol=8.07e-02 90.40 % -> 3.30 s
38 : 22.800 μm, [9.82e-03; 1.00e+00], 1:0.877 s, tol=7.54e-02 91.20 % -> 3.13 s
39 : 23.000 μm, [1.01e-02; 1.00e+00], 1:0.883 s, tol=7.34e-02 92.00 % -> 2.92 s
40 : 23.200 μm, [1.04e-02; 1.00e+00], 1:0.889 s, tol=7.26e-02 92.80 % -> 2.69 s
41 : 23.400 μm, [1.06e-02; 1.00e+00], 1:0.887 s, tol=7.46e-02 93.60 % -> 2.43 s
42 : 23.600 μm, [1.08e-02; 1.00e+00], 1:0.872 s, tol=7.51e-02 94.40 % -> 2.12 s
43 : 23.800 μm, [1.10e-02; 1.00e+00], 1:0.863 s, tol=7.56e-02 95.20 % -> 1.83 s
44 : 24.000 μm, [1.11e-02; 1.00e+00], 1:0.877 s, tol=7.53e-02 96.00 % -> 1.57 s
45 : 24.200 μm, [1.12e-02; 1.00e+00], 1:0.858 s, tol=7.52e-02 96.80 % -> 1.25 s
46 : 24.400 μm, [1.14e-02; 1.00e+00], 1:0.895 s, tol=7.46e-02 97.60 % -> 990.34 ms
47 : 24.600 μm, [1.15e-02; 1.00e+00], 1:0.894 s, tol=7.42e-02 98.40 % -> 668.48 ms
48 : 24.800 μm, [1.16e-02; 1.00e+00], 1:0.897 s, tol=7.27e-02 99.20 % -> 340.00 ms
49 : 25.000 μm, [1.17e-02; 1.00e+00], 1:0.856 s, tol=7.13e-02 100.00 % -> -0.00 µs
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/force-displacement.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/simulation.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.3/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/summary.txt
Generate movie 01/24 (4.17 %) 2.99 s
Generate movie 02/24 (8.33 %) 2.32 s
Generate movie 03/24 (12.50 %) 2.20 s
Generate movie 04/24 (16.67 %) 2.08 s
Generate movie 05/24 (20.83 %) 2.07 s
Generate movie 06/24 (25.00 %) 1.92 s
Generate movie 07/24 (29.17 %) 1.78 s
Generate movie 08/24 (33.33 %) 1.74 s
Generate movie 09/24 (37.50 %) 1.55 s
Generate movie 10/24 (41.67 %) 1.45 s
Generate movie 11/24 (45.83 %) 1.33 s
Generate movie 12/24 (50.00 %) 1.32 s
Generate movie 13/24 (54.17 %) 1.13 s
Generate movie 14/24 (58.33 %) 1.07 s
Generate movie 15/24 (62.50 %) 952.01 ms
Generate movie 16/24 (66.67 %) 837.46 ms
Generate movie 17/24 (70.83 %) 741.29 ms
Generate movie 18/24 (75.00 %) 647.76 ms
Generate movie 19/24 (79.17 %) 537.03 ms
Generate movie 20/24 (83.33 %) 425.95 ms
Generate movie 21/24 (87.50 %) 322.93 ms
Generate movie 22/24 (91.67 %) 219.85 ms
Generate movie 23/24 (95.83 %) 108.74 ms
Generate movie 24/24 (100.00 %) 0.00 µs
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import (
17 Display,
18 Folder,
19 Models,
20 Tic,
21 ElemType,
22 Mesh,
23 Simulations,
24 PyVista,
25 Paraview,
26 )
27 from EasyFEA.Geoms import Domain, Circle
28
29 import multiprocessing
30
31 # Display.Clear()
32
33
34 # ----------------------------------------------
35 # Configuration
36 # ----------------------------------------------
37
38 # simu options
39 doSimu = True
40 meshTest = True
41 optimMesh = True
42 useParallel = False
43 nProcs = 4 # number of processes in parallel
44
45 # outputs
46 folder = Folder.Results_Dir()
47 plotMesh = False
48 plotIter = False
49 plotResult = True
50 plotEnergy = False
51 showFig = True
52
53 saveParaview = False
54 makeMovie = True
55
56 # models
57
58 # splits = ["Bourdin","Amor","Miehe","Stress"] # Splits Isotropes
59 # splits = ["He","AnisotStrain","AnisotStress","Zhang"] # Splits Anisotropes
60 # splits = ["Bourdin","Amor","Miehe","Stress","He","AnisotStrain","AnisotStress","Zhang"]
61 # splits = ["Zhang"]
62 # splits = ["AnisotStrain","AnisotStress","Zhang"]
63 splits = ["Miehe"]
64
65 regus = ["AT2"] # ["AT1", "AT2"]
66 # regus = ["AT1", "AT2"]
67
68 l0 = 0.12e-3
69
70 # convergence
71 solver = (
72 Models.PhaseField.SolverType.History
73 ) # ["History", "HistoryDamage", "BoundConstrain"]
74 maxIter = 1000
75 tolConv = 1e-0
76
77 # ----------------------------------------------
78 # Mesh
79 # ----------------------------------------------
80
81
82 def DoMesh(
83 L: float, h: float, diam: float, thickness: float, l0: float, split: str
84 ) -> Mesh:
85 clC = l0 * 2 if meshTest else l0 / 2
86 if optimMesh:
87 clD = l0 * 4
88 refineZone = diam * 1.5 / 2
89 if split in ["Bourdin", "Amor"]:
90 refineGeom = Domain((0, h / 2 - refineZone), (L, h / 2 + refineZone), clC)
91 else:
92 refineGeom = Domain((L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC)
93 else:
94 clD = l0 if meshTest else l0 / 2
95 refineGeom = None
96
97 domain = Domain((0, 0), (L, h), clD)
98 circle = Circle((L / 2, h / 2), diam, clD, isHollow=True)
99
100 # ax = Display.Init_Axes()
101 # domain.Plot(ax, color="k", plotPoints=False)
102 # circle.Plot(ax, color="k", plotPoints=False)
103 # # if refineGeom != None:
104 # # refineGeom.Plot(ax, color='k', plotPoints=False)
105 # # ax.scatter(((L+diam)/2, L/2), (h/2, (h+diam)/2), c='k')
106 # ax.axis("off")
107 # Display.Save_fig(folder, "sample", True)
108
109 mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom])
110
111 # ax = Display.Plot_Mesh(mesh, lw=0.3, facecolors="white")
112 # ax.axis("off")
113 # ax.set_title("")
114 # Display.Save_fig(folder, "mesh", transparent=True)
115
116 return mesh
117
118
119 # ----------------------------------------------
120 # Do Simu
121 # ----------------------------------------------
122 def DoSimu(split: str, regu: str) -> str:
123
124 folder_save = Simulations.PhaseField.Folder(
125 folder, "Elas_Isot", split, regu, "DP", tolConv, solver, meshTest, optimMesh
126 )
127
128 Display.MyPrint(folder_save, "green")
129
130 # ----------------------------------------------
131 # Geom
132 # ----------------------------------------------
133
134 L = 15e-3
135 h = 30e-3
136 thickness = 1
137 diam = 6e-3
138
139 # load units
140 unitU = "μm"
141 unitF = "kN/mm"
142 unit = 1e6
143
144 # ----------------------------------------------
145 # Mesh
146 # ----------------------------------------------
147
148 if doSimu:
149 mesh = DoMesh(L, h, diam, thickness, l0, split)
150
151 # Get Nodes
152 nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
153 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
154 nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0))
155 nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"])
156 nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h)
157
158 # ----------------------------------------------
159 # Boundary conditions
160 # ----------------------------------------------
161
162 threshold = 0.6
163 u_max = 25e-6
164 # u_max = 35e-6
165
166 uinc0 = 8e-7 if meshTest else 8e-8
167 uinc1 = 2e-7 if meshTest else 2e-8
168
169 config = f"""
170 while ud <= u_max:
171
172 ud += uinc0 if simu.damage.max() < threshold else uinc1
173
174 u_max = {u_max}
175 uinc0 = {uinc0:.1e} (simu.damage.max() < {threshold})
176 uinc1 = {uinc1:.1e}
177
178 simu.add_dirichlet(nodes_lower, [0], ["y"])
179 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
180 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
181 if dim == 3:
182 simu.add_dirichlet(nodes_y0z0, [0], ["z"])
183 """
184
185 # ----------------------------------------------
186 # Material
187 # ----------------------------------------------
188 E = 12e9
189 v = 0.3
190 planeStress = False
191 material = Models.Elastic.Isotropic(2, E, v, planeStress, thickness)
192
193 gc = 1.4
194 pfm = Models.PhaseField(material, split, regu, gc, l0, solver=solver)
195
196 # ----------------------------------------------
197 # Simulation
198 # ----------------------------------------------
199 simu = Simulations.PhaseField(mesh, pfm, verbosity=False)
200
201 simu.Results_Set_Bc_Summary(config)
202
203 dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"])
204
205 def Apply_BC(ud: float):
206 simu.Bc_Init()
207 simu.add_dirichlet(nodes_lower, [0], ["y"])
208 simu.add_dirichlet(nodes_x0y0, [0], ["x"])
209 simu.add_dirichlet(nodes_upper, [-ud], ["y"])
210
211 # INIT
212 displacement = []
213 force = []
214 ud = -uinc0
215 iter = 0
216 nDetect = 0
217
218 if plotIter:
219 axIter = Display.Plot_Result(simu, "damage", nodeValues=True)
220
221 force = np.asarray(force)
222 displacement = np.asarray(displacement)
223 _, axLoad = Display.Plot_Force_Displacement(
224 force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]"
225 )
226
227 while ud <= u_max:
228 iter += 1
229 ud += uinc0 if simu.damage.max() < threshold else uinc1
230
231 Apply_BC(ud)
232
233 u, d, Ku, convergence = simu.Solve(tolConv, maxIter)
234 simu.Save_Iter()
235
236 # stop if the simulation does not converge
237 if not convergence:
238 break
239
240 f = np.sum(Ku[dofsY_upper, :] @ u)
241
242 simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True)
243
244 # Detect damaged edges
245 if np.any(d[nodes_edges] >= 1):
246 nDetect += 1
247 if nDetect == 10:
248 break
249
250 displacement = np.concatenate((displacement, [ud]))
251 force = np.concatenate((force, [f]))
252
253 if plotIter:
254 Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter)
255 plt.figure(axIter.figure)
256 plt.pause(1e-12)
257
258 force = np.asarray(force)
259 displacement = np.asarray(displacement)
260 Display.Plot_Force_Displacement(
261 force / unit,
262 displacement * unit,
263 f"ud [{unitU}]",
264 f"f [{unitF}]",
265 ax=axLoad,
266 )[1]
267 plt.figure(axLoad.figure)
268 plt.pause(1e-12)
269
270 # ----------------------------------------------
271 # Saving
272 # ----------------------------------------------
273 force = np.asarray(force)
274 displacement = np.asarray(displacement)
275 print()
276 Simulations.Save_pickle(
277 (force, displacement), folder_save, "force-displacement"
278 )
279 simu.Save(folder_save)
280
281 else:
282 # ----------------------------------------------
283 # Load
284 # ----------------------------------------------
285 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
286 force, displacement = Simulations.Load_pickle(folder_save, "force-displacement")
287
288 # ----------------------------------------------
289 # Results
290 # ---------------------------------------------
291 if plotEnergy:
292 Display.Plot_Energy(simu, force, displacement, N=400, folder=folder_save)
293
294 if plotResult:
295 Display.Plot_Result(
296 simu,
297 "damage",
298 nodeValues=True,
299 colorbarIsClose=True,
300 folder=folder_save,
301 filename="damage",
302 )
303 Display.Plot_Mesh(simu)
304 Display.Plot_BoundaryConditions(simu)
305 Display.Plot_Iter_Summary(simu, folder_save, None, None)
306 Display.Plot_Force_Displacement(
307 force / unit,
308 displacement * unit,
309 f"ud [{unitU}]",
310 f"f [{unitF}]",
311 folder_save,
312 )
313
314 if plotMesh:
315 Display.Plot_Mesh(mesh)
316
317 if saveParaview:
318 Paraview.Save_simu(simu, folder_save)
319
320 if makeMovie:
321 simu.Set_Iter(-1)
322 depMax = simu.Result("displacement_norm").max()
323 deformFactor = L * 0.05 / depMax
324
325 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
326
327 def Func(plotter, iter):
328 simu.Set_Iter(iterations[iter])
329
330 grid = PyVista._pvMesh(simu, "damage", deformFactor)
331
332 tresh = grid.threshold((0, 0.8))
333
334 PyVista.Plot(
335 tresh,
336 "damage",
337 deformFactor,
338 plotMesh=True,
339 plotter=plotter,
340 clim=(0, 1),
341 )
342
343 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
344
345 if doSimu:
346 Tic.Plot_History(folder_save, details=False)
347
348 if showFig:
349 plt.show()
350
351 # Display.plt.close("all")
352 Tic.Clear()
353
354 return folder_save
355
356
357 if __name__ == "__main__":
358 # generates configs
359 Splits = []
360 Regus = []
361 for split in splits.copy():
362 for regu in regus.copy():
363 Splits.append(split)
364 Regus.append(regu)
365
366 list_folder: list[str] = []
367 if useParallel:
368 items = [(split, regu) for split, regu in zip(Splits, Regus)]
369 with multiprocessing.Pool(nProcs) as pool:
370 for result in pool.starmap(DoSimu, items):
371 list_folder.append(result)
372 else:
373 for split, regu in zip(Splits, Regus):
374 list_folder.append(DoSimu(split, regu))
Total running time of the script: (0 minutes 47.845 seconds)