Note

Go to the end to download the full example code.

LShape#

Damage simulation for a L-part.

LShape
  • $\phi$
  • TRI3: Ne = 4234, Nn = 2178
  • Boundary conditions
  • LShape
  • LShape
  • Summary
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/envs/v1.5.5/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.213 s, tol=0.00e+00
   2 : 0.040 mm, [0.00e+00; 0.00e+00], 1:0.222 s, tol=1.00e+00  4.00 % -> 5.32 s
   3 : 0.080 mm, [0.00e+00; 0.00e+00], 1:0.222 s, tol=9.95e-01  8.00 % -> 5.11 s
   4 : 0.120 mm, [0.00e+00; 0.00e+00], 1:0.222 s, tol=5.56e-01  12.00 % -> 4.89 s
   5 : 0.160 mm, [0.00e+00; 0.00e+00], 1:0.221 s, tol=4.38e-01  16.00 % -> 4.64 s
   6 : 0.200 mm, [0.00e+00; 0.00e+00], 1:0.237 s, tol=3.60e-01  20.00 % -> 4.74 s
   7 : 0.240 mm, [0.00e+00; 0.00e+00], 1:0.223 s, tol=3.06e-01  24.00 % -> 4.24 s
   8 : 0.280 mm, [0.00e+00; 8.98e-03], 1:0.255 s, tol=2.65e-01  28.00 % -> 4.59 s
   9 : 0.320 mm, [0.00e+00; 4.57e-02], 1:0.248 s, tol=2.32e-01  32.00 % -> 4.22 s
  10 : 0.360 mm, [0.00e+00; 1.10e-01], 1:0.247 s, tol=2.07e-01  36.00 % -> 3.96 s
  11 : 0.400 mm, [0.00e+00; 2.07e-01], 1:0.248 s, tol=1.86e-01  40.00 % -> 3.72 s
  12 : 0.440 mm, [0.00e+00; 3.48e-01], 1:0.245 s, tol=1.69e-01  44.00 % -> 3.43 s
  13 : 0.480 mm, [0.00e+00; 5.46e-01], 1:0.255 s, tol=1.55e-01  48.00 % -> 3.32 s
  14 : 0.520 mm, [0.00e+00; 7.78e-01], 1:0.261 s, tol=1.62e-01  52.00 % -> 3.13 s
  15 : 0.540 mm, [0.00e+00; 9.48e-01], 1:0.172 s, tol=2.22e-01  54.00 % -> 2.05 s
  16 : 0.560 mm, [0.00e+00; 9.93e-01], 1:0.248 s, tol=3.80e-01  56.00 % -> 2.93 s
  17 : 0.580 mm, [0.00e+00; 1.00e+00], 1:0.244 s, tol=4.44e-01  58.00 % -> 2.83 s
  18 : 0.600 mm, [0.00e+00; 1.00e+00], 1:0.252 s, tol=3.82e-01  60.00 % -> 2.86 s
  19 : 0.620 mm, [0.00e+00; 1.00e+00], 1:0.247 s, tol=3.28e-01  62.00 % -> 2.73 s
  20 : 0.640 mm, [0.00e+00; 1.00e+00], 1:0.261 s, tol=2.91e-01  64.00 % -> 2.79 s
  21 : 0.660 mm, [0.00e+00; 1.00e+00], 1:0.253 s, tol=2.78e-01  66.00 % -> 2.61 s
  22 : 0.680 mm, [0.00e+00; 1.00e+00], 1:0.261 s, tol=2.54e-01  68.00 % -> 2.58 s
  23 : 0.700 mm, [0.00e+00; 1.00e+00], 1:0.227 s, tol=2.53e-01  70.00 % -> 2.14 s
  24 : 0.720 mm, [0.00e+00; 1.00e+00], 1:0.252 s, tol=2.02e-01  72.00 % -> 2.26 s
  25 : 0.740 mm, [0.00e+00; 1.00e+00], 1:0.263 s, tol=1.76e-01  74.00 % -> 2.22 s
  26 : 0.760 mm, [0.00e+00; 1.00e+00], 1:0.251 s, tol=1.87e-01  76.00 % -> 1.98 s
  27 : 0.780 mm, [0.00e+00; 1.00e+00], 1:0.260 s, tol=1.72e-01  78.00 % -> 1.91 s
  28 : 0.800 mm, [0.00e+00; 1.00e+00], 1:0.260 s, tol=1.65e-01  80.00 % -> 1.76 s
  29 : 0.820 mm, [0.00e+00; 1.00e+00], 1:0.245 s, tol=1.90e-01  82.00 % -> 1.51 s
  30 : 0.840 mm, [0.00e+00; 1.00e+00], 1:0.252 s, tol=1.65e-01  84.00 % -> 1.39 s
  31 : 0.860 mm, [0.00e+00; 1.00e+00], 1:0.217 s, tol=1.51e-01  86.00 % -> 1.06 s
  32 : 0.880 mm, [0.00e+00; 1.00e+00], 1:0.252 s, tol=1.15e-01  88.00 % -> 1.07 s
  33 : 0.900 mm, [0.00e+00; 1.00e+00], 1:0.250 s, tol=1.10e-01  90.00 % -> 889.59 ms
  34 : 0.920 mm, [0.00e+00; 1.00e+00], 1:0.250 s, tol=9.93e-02  92.00 % -> 717.41 ms
  35 : 0.940 mm, [0.00e+00; 1.00e+00], 1:0.247 s, tol=9.94e-02  94.00 % -> 536.20 ms
  36 : 0.960 mm, [0.00e+00; 1.00e+00], 1:0.254 s, tol=9.12e-02  96.00 % -> 370.73 ms
  37 : 0.980 mm, [0.00e+00; 1.00e+00], 1:0.257 s, tol=1.02e-01  98.00 % -> 188.95 ms
  38 : 1.000 mm, [0.00e+00; 1.00e+00], 1:0.251 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

Generate movie 0/37
Generate movie 1/37 (2.70 %) 4.96 s
Generate movie 2/37 (5.41 %) 4.83 s
Generate movie 3/37 (8.11 %) 4.71 s
Generate movie 4/37 (10.81 %) 4.60 s
Generate movie 5/37 (13.51 %) 4.45 s
Generate movie 6/37 (16.22 %) 4.29 s
Generate movie 7/37 (18.92 %) 4.14 s
Generate movie 8/37 (21.62 %) 4.03 s
Generate movie 9/37 (24.32 %) 3.89 s
Generate movie 10/37 (27.03 %) 3.72 s
Generate movie 11/37 (29.73 %) 3.62 s
Generate movie 12/37 (32.43 %) 3.47 s
Generate movie 13/37 (35.14 %) 3.31 s
Generate movie 14/37 (37.84 %) 3.18 s
Generate movie 15/37 (40.54 %) 3.06 s
Generate movie 16/37 (43.24 %) 2.91 s
Generate movie 17/37 (45.95 %) 2.78 s
Generate movie 18/37 (48.65 %) 2.63 s
Generate movie 19/37 (51.35 %) 2.51 s
Generate movie 20/37 (54.05 %) 2.37 s
Generate movie 21/37 (56.76 %) 2.21 s
Generate movie 22/37 (59.46 %) 2.10 s
Generate movie 23/37 (62.16 %) 1.94 s
Generate movie 24/37 (64.86 %) 1.80 s
Generate movie 25/37 (67.57 %) 1.65 s
Generate movie 26/37 (70.27 %) 1.52 s
Generate movie 27/37 (72.97 %) 1.37 s
Generate movie 28/37 (75.68 %) 1.25 s
Generate movie 29/37 (78.38 %) 1.10 s
Generate movie 30/37 (81.08 %) 971.09 ms
Generate movie 31/37 (83.78 %) 826.97 ms
Generate movie 32/37 (86.49 %) 696.05 ms
Generate movie 33/37 (89.19 %) 560.84 ms
Generate movie 34/37 (91.89 %) 426.70 ms
Generate movie 35/37 (94.59 %) 275.22 ms
Generate movie 36/37 (97.30 %) 138.62 ms
Generate movie 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._pvMesh(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 18.367 seconds)

Gallery generated by Sphinx-Gallery