Note
Go to the end to download the full example code.
LShape#
Damage simulation for a L-part.
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/envs/v1.5.4/lib/python3.11/site-packages/results/PhaseField/LShape_2D/Test/Isot_Miehe_AT1_optimMesh nL=50/Test/Miehe_AT1_CP_optimMesh nL=50 1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.046 s, tol=0.00e+00
2 : 0.040 mm, [0.00e+00; 0.00e+00], 1:0.035 s, tol=1.00e+00 4.00 % -> 844.63 ms
3 : 0.080 mm, [0.00e+00; 0.00e+00], 1:0.036 s, tol=9.95e-01 8.00 % -> 820.99 ms
4 : 0.120 mm, [0.00e+00; 0.00e+00], 1:0.036 s, tol=5.56e-01 12.00 % -> 784.07 ms
5 : 0.160 mm, [0.00e+00; 0.00e+00], 1:0.036 s, tol=4.37e-01 16.00 % -> 756.20 ms
6 : 0.200 mm, [0.00e+00; 0.00e+00], 1:0.036 s, tol=3.60e-01 20.00 % -> 713.62 ms
7 : 0.240 mm, [0.00e+00; 0.00e+00], 1:0.036 s, tol=3.06e-01 24.00 % -> 681.41 ms
8 : 0.280 mm, [0.00e+00; 8.98e-03], 1:0.036 s, tol=2.65e-01 28.00 % -> 648.51 ms
9 : 0.320 mm, [0.00e+00; 4.57e-02], 1:0.036 s, tol=2.32e-01 32.00 % -> 609.99 ms
10 : 0.360 mm, [0.00e+00; 1.10e-01], 1:0.037 s, tol=2.07e-01 36.00 % -> 595.05 ms
11 : 0.400 mm, [0.00e+00; 2.07e-01], 1:0.037 s, tol=1.86e-01 40.00 % -> 560.52 ms
12 : 0.440 mm, [0.00e+00; 3.48e-01], 1:0.037 s, tol=1.69e-01 44.00 % -> 524.24 ms
13 : 0.480 mm, [0.00e+00; 5.46e-01], 1:0.036 s, tol=1.55e-01 48.00 % -> 472.71 ms
14 : 0.520 mm, [0.00e+00; 7.78e-01], 1:0.036 s, tol=1.62e-01 52.00 % -> 433.57 ms
15 : 0.540 mm, [0.00e+00; 9.48e-01], 1:0.036 s, tol=2.22e-01 54.00 % -> 433.07 ms
16 : 0.560 mm, [0.00e+00; 9.93e-01], 1:0.036 s, tol=3.80e-01 56.00 % -> 424.88 ms
17 : 0.580 mm, [0.00e+00; 1.00e+00], 1:0.037 s, tol=4.44e-01 58.00 % -> 424.75 ms
18 : 0.600 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=3.82e-01 60.00 % -> 409.46 ms
19 : 0.620 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=3.28e-01 62.00 % -> 399.91 ms
20 : 0.640 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=2.91e-01 64.00 % -> 384.29 ms
21 : 0.660 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=2.78e-01 66.00 % -> 370.11 ms
22 : 0.680 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=2.54e-01 68.00 % -> 356.99 ms
23 : 0.700 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=2.53e-01 70.00 % -> 338.77 ms
24 : 0.720 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=2.02e-01 72.00 % -> 324.39 ms
25 : 0.740 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.76e-01 74.00 % -> 303.96 ms
26 : 0.760 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.87e-01 76.00 % -> 284.47 ms
27 : 0.780 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.72e-01 78.00 % -> 263.49 ms
28 : 0.800 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.65e-01 80.00 % -> 242.71 ms
29 : 0.820 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.90e-01 82.00 % -> 222.86 ms
30 : 0.840 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.65e-01 84.00 % -> 199.27 ms
31 : 0.860 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.51e-01 86.00 % -> 175.43 ms
32 : 0.880 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.15e-01 88.00 % -> 153.61 ms
33 : 0.900 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.10e-01 90.00 % -> 129.35 ms
34 : 0.920 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=9.93e-02 92.00 % -> 104.23 ms
35 : 0.940 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=9.94e-02 94.00 % -> 78.96 ms
36 : 0.960 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=9.12e-02 96.00 % -> 53.01 ms
37 : 0.980 mm, [0.00e+00; 1.00e+00], 1:0.036 s, tol=1.02e-01 98.00 % -> 26.55 ms
38 : 1.000 mm, [0.00e+00; 1.00e+00], 1:0.037 s, tol=9.35e-02 100.00 % -> -0.00 µs
Saved:
/results/PhaseField/LShape_2D/Test/Isot_Miehe_AT1_optimMesh nL=50/Test/Miehe_AT1_CP_optimMesh nL=50/force-displacement.pickle
Saved:
/results/PhaseField/LShape_2D/Test/Isot_Miehe_AT1_optimMesh nL=50/Test/Miehe_AT1_CP_optimMesh nL=50/simulation.pickle
Saved:
/results/PhaseField/LShape_2D/Test/Isot_Miehe_AT1_optimMesh nL=50/Test/Miehe_AT1_CP_optimMesh nL=50/summary.txt
Loaded:
/results/PhaseField/LShape_2D/Test/Isot_Miehe_AT1_optimMesh nL=50/Test/Miehe_AT1_CP_optimMesh nL=50/force-displacement.pickle
Movie_func 0/37
Movie_func 1/37 (2.70 %) 4.40 s
Movie_func 2/37 (5.41 %) 4.03 s
Movie_func 3/37 (8.11 %) 4.00 s
Movie_func 4/37 (10.81 %) 3.86 s
Movie_func 5/37 (13.51 %) 3.70 s
Movie_func 6/37 (16.22 %) 3.68 s
Movie_func 7/37 (18.92 %) 3.75 s
Movie_func 8/37 (21.62 %) 3.35 s
Movie_func 9/37 (24.32 %) 3.27 s
Movie_func 10/37 (27.03 %) 3.18 s
Movie_func 11/37 (29.73 %) 3.13 s
Movie_func 12/37 (32.43 %) 2.89 s
Movie_func 13/37 (35.14 %) 2.78 s
Movie_func 14/37 (37.84 %) 2.72 s
Movie_func 15/37 (40.54 %) 2.55 s
Movie_func 16/37 (43.24 %) 2.43 s
Movie_func 17/37 (45.95 %) 2.34 s
Movie_func 18/37 (48.65 %) 2.19 s
Movie_func 19/37 (51.35 %) 2.09 s
Movie_func 20/37 (54.05 %) 2.01 s
Movie_func 21/37 (56.76 %) 1.87 s
Movie_func 22/37 (59.46 %) 1.73 s
Movie_func 23/37 (62.16 %) 1.61 s
Movie_func 24/37 (64.86 %) 1.50 s
Movie_func 25/37 (67.57 %) 1.38 s
Movie_func 26/37 (70.27 %) 1.27 s
Movie_func 27/37 (72.97 %) 1.17 s
Movie_func 28/37 (75.68 %) 1.06 s
Movie_func 29/37 (78.38 %) 914.72 ms
Movie_func 30/37 (81.08 %) 806.77 ms
Movie_func 31/37 (83.78 %) 686.85 ms
Movie_func 32/37 (86.49 %) 579.38 ms
Movie_func 33/37 (89.19 %) 458.91 ms
Movie_func 34/37 (91.89 %) 344.61 ms
Movie_func 35/37 (94.59 %) 234.47 ms
Movie_func 36/37 (97.30 %) 114.55 ms
Movie_func 37/37 (100.00 %) 0.00 µs
12 from EasyFEA import (
13 Display,
14 Folder,
15 Models,
16 plt,
17 np,
18 Tic,
19 ElemType,
20 Simulations,
21 Paraview,
22 PyVista,
23 )
24 from EasyFEA.Geoms import Point, Points, Domain, Circle
25
26 if __name__ == "__main__":
27 Display.Clear()
28
29 # ----------------------------------------------
30 # Configuration
31 # ----------------------------------------------
32
33 # simu options
34 doSimu = True
35 meshTest = True
36 optimMesh = True
37
38 # outputs
39 pltIter = False
40 pltLoad = False
41 makeMovie = True
42 makeParaview = False
43
44 # geom
45 dim = 2
46 L = 250 # mm
47 nL = 50
48 ep = 100
49
50 # model
51 E = 2e4 # MPa
52 v = 0.18
53
54 split = "Miehe"
55 regu = "AT1"
56 Gc = 130 # J/m2
57 Gc *= 1000 / 1e6 # mJ/mm2
58
59 # convergence
60 tolConv = 1e-0
61 convOption = 2
62
63 # load
64 uMax = 1 # mm
65 inc0 = uMax / 25
66 inc1 = inc0 / 2
67
68 config = f"""
69 uMax = {uMax}
70
71 inc0 = {inc0}
72 inc1 = {inc1}
73
74 while ud <= uMax:
75
76 ud += inc0 if simu.damage.max() < 0.6 else inc1
77
78 simu.add_dirichlet(nodes_circle, [0], ['d'], "damage")
79 simu.add_dirichlet(nodes_y0, [0]*dim, simu.Get_dofs())
80 simu.add_dirichlet(nodes_load, [ud], ['y'])
81 """
82
83 folder = Folder.PhaseField_Folder(
84 f"LShape_{dim}D",
85 "Isot",
86 split,
87 regu,
88 "",
89 tolConv,
90 "",
91 meshTest,
92 optimMesh,
93 nL=nL,
94 )
95
96 # ----------------------------------------------
97 # Mesh
98 # ----------------------------------------------
99 l0 = L / nL
100
101 if meshTest:
102 hC = l0
103 else:
104 hC = 0.5
105 # hC = 0.25
106
107 p1 = Point()
108 p2 = Point(L, 0)
109 p3 = Point(L, L)
110 p4 = Point(2 * L - 30, L)
111 p5 = Point(2 * L, L)
112 p6 = Point(2 * L, 2 * L)
113 p7 = Point(0, 2 * L)
114 if optimMesh:
115 # hauteur zone rafinée
116 h = 100
117 refineDomain = Domain(Point(0, L - h / 3), Point(L + h / 3, L + h), hC)
118 hD = hC * 5
119 else:
120 refineDomain = None
121 hD = hC
122
123 contour = Points([p1, p2, p3, p4, p5, p6, p7], hD)
124
125 circle = Circle(p5, 100)
126
127 if dim == 2:
128 mesh = contour.Mesh_2D([], ElemType.TRI3, refineGeoms=[refineDomain])
129 else:
130 mesh = contour.Mesh_Extrude(
131 [], [0, 0, -ep], [3], ElemType.HEXA8, refineGeoms=[refineDomain]
132 )
133
134 # Display.Plot_Mesh(mesh)
135 # Display.Plot_Tags(mesh)
136 # from EasyFEA import PyVista
137 # PyVista.Plot_Mesh(mesh).show()
138
139 nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
140 nodes_load = mesh.Nodes_Conditions(lambda x, y, z: (y == L) & (x >= 2 * L - 30))
141 node3 = mesh.Nodes_Point(p3)
142 node4 = mesh.Nodes_Point(p4)
143 nodes_circle = mesh.Nodes_Cylinder(circle, [0, 0, ep])
144 nodes_edges = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) | (y == 0))
145
146 # ----------------------------------------------
147 # Simulation
148 # ----------------------------------------------
149 material = Models.ElasIsot(dim, E, v, True, ep)
150 pfm = Models.PhaseField(material, split, regu, Gc, l0)
151
152 folder_save = Folder.PhaseField_Folder(
153 folder,
154 "",
155 pfm.split,
156 pfm.regularization,
157 "CP",
158 tolConv,
159 "",
160 meshTest,
161 optimMesh,
162 nL=nL,
163 )
164
165 Display.MyPrint(folder_save, "green")
166
167 if doSimu:
168 simu = Simulations.PhaseFieldSimu(mesh, pfm)
169 simu.Results_Set_Bc_Summary(config)
170
171 dofsY_load = simu.Bc_dofs_nodes(nodes_load, ["y"])
172
173 if pltIter:
174 axIter = Display.Plot_Result(simu, "damage")
175
176 axLoad = Display.Init_Axes()
177 axLoad.set_xlabel("displacement [mm]")
178 axLoad.set_ylabel("load [kN]")
179
180 displacement = []
181 force = []
182 ud = -inc0
183 iter = -1
184
185 while ud <= uMax:
186 # update displacement
187 iter += 1
188 ud += inc0 if simu.damage.max() < 0.6 else inc1
189
190 # update boundary conditions
191 simu.Bc_Init()
192 simu.add_dirichlet(nodes_circle, [0], ["d"], "damage")
193 simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
194 simu.add_dirichlet(nodes_load, [ud], ["y"])
195
196 # solve
197 u, d, Kglob, convergence = simu.Solve(tolConv, 500, convOption)
198 # calc load
199 fr = np.sum(Kglob[dofsY_load, :] @ u)
200
201 # saves load and displacement
202 displacement.append(ud)
203 force.append(fr)
204
205 # print iter
206 simu.Results_Set_Iteration_Summary(iter, ud, "mm", ud / uMax, True)
207
208 # saves iteration
209 simu.Save_Iter()
210
211 if pltIter:
212 plt.figure(axIter.figure)
213 Display.Plot_Result(simu, "damage", ax=axIter)
214 plt.pause(1e-12)
215
216 plt.figure(axLoad.figure)
217 axLoad.scatter(ud, fr / 1000, c="black")
218 plt.pause(1e-12)
219
220 if not convergence or np.max(d[nodes_edges]) >= 1:
221 # stops simulation if damage occurs on edges or convergence has not been reached
222 break
223
224 # saves load and displacement
225 displacement = np.asarray(displacement)
226 force = np.asarray(force)
227 Simulations.Save_pickle(
228 (force, displacement), folder_save, "force-displacement"
229 )
230
231 # saves the simulation
232 simu.Save(folder_save)
233
234 else:
235 simu = Simulations.Load_Simu(folder_save)
236 mesh = simu.mesh
237
238 force, displacement = Simulations.Load_pickle(folder_save, "force-displacement")
239
240 # ----------------------------------------------
241 # Results
242 # ----------------------------------------------
243 Display.Plot_Result(simu, "damage", folder=folder_save)
244 Display.Plot_Mesh(simu)
245 Display.Plot_BoundaryConditions(simu)
246
247 axLoad = Display.Init_Axes()
248 axLoad.set_xlabel("displacement [mm]")
249 axLoad.set_ylabel("load [kN]")
250 axLoad.plot(displacement, force / 1000, c="blue")
251 Display.Save_fig(folder_save, "forcedep")
252
253 Display.Plot_Iter_Summary(simu, folder_save)
254
255 if makeMovie:
256 simu.Set_Iter(-1)
257 depMax = simu.Result("displacement_norm").max()
258 deformFactor = L * 0.05 / depMax
259
260 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
261
262 def Func(plotter, iter):
263 simu.Set_Iter(iterations[iter])
264
265 grid = PyVista._pyVistaMesh(simu, "damage", deformFactor)
266
267 tresh = grid.threshold((0, 0.8))
268
269 PyVista.Plot(
270 tresh,
271 "damage",
272 deformFactor,
273 plotMesh=True,
274 plotter=plotter,
275 clim=(0, 1),
276 )
277
278 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
279
280 if makeParaview:
281 Paraview.Save_simu(simu, folder_save)
282
283 if doSimu:
284 Tic.Plot_History(folder_save, False)
285
286 plt.show()
Total running time of the script: (0 minutes 9.644 seconds)