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/checkouts/v1.7.1/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_CP_optimMesh nL=50   1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.212 s, tol=0.00e+00
   2 : 0.040 mm, [0.00e+00; 0.00e+00], 1:0.251 s, tol=1.00e+00  4.00 % -> 6.02 s
   3 : 0.080 mm, [0.00e+00; 0.00e+00], 1:0.174 s, tol=9.95e-01  8.00 % -> 4.00 s
   4 : 0.120 mm, [0.00e+00; 0.00e+00], 1:0.162 s, tol=5.56e-01  12.00 % -> 3.56 s
   5 : 0.160 mm, [0.00e+00; 0.00e+00], 1:0.242 s, tol=4.38e-01  16.00 % -> 5.08 s
   6 : 0.200 mm, [0.00e+00; 0.00e+00], 1:0.164 s, tol=3.60e-01  20.00 % -> 3.29 s
   7 : 0.240 mm, [0.00e+00; 0.00e+00], 1:0.265 s, tol=3.06e-01  24.00 % -> 5.04 s
   8 : 0.280 mm, [0.00e+00; 8.98e-03], 1:0.157 s, tol=2.65e-01  28.00 % -> 2.83 s
   9 : 0.320 mm, [0.00e+00; 4.57e-02], 1:0.252 s, tol=2.32e-01  32.00 % -> 4.28 s
  10 : 0.360 mm, [0.00e+00; 1.10e-01], 1:0.247 s, tol=2.07e-01  36.00 % -> 3.95 s
  11 : 0.400 mm, [0.00e+00; 2.07e-01], 1:0.248 s, tol=1.86e-01  40.00 % -> 3.73 s
  12 : 0.440 mm, [0.00e+00; 3.48e-01], 1:0.250 s, tol=1.69e-01  44.00 % -> 3.50 s
  13 : 0.480 mm, [0.00e+00; 5.46e-01], 1:0.180 s, tol=1.55e-01  48.00 % -> 2.34 s
  14 : 0.520 mm, [0.00e+00; 7.78e-01], 1:0.251 s, tol=1.62e-01  52.00 % -> 3.01 s
  15 : 0.540 mm, [0.00e+00; 9.48e-01], 1:0.250 s, tol=2.22e-01  54.00 % -> 2.98 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.254 s, tol=4.44e-01  58.00 % -> 2.94 s
  18 : 0.600 mm, [0.00e+00; 1.00e+00], 1:0.251 s, tol=3.82e-01  60.00 % -> 2.85 s
  19 : 0.620 mm, [0.00e+00; 1.00e+00], 1:0.248 s, tol=3.28e-01  62.00 % -> 2.74 s
  20 : 0.640 mm, [0.00e+00; 1.00e+00], 1:0.250 s, tol=2.91e-01  64.00 % -> 2.67 s
  21 : 0.660 mm, [0.00e+00; 1.00e+00], 1:0.251 s, tol=2.78e-01  66.00 % -> 2.59 s
  22 : 0.680 mm, [0.00e+00; 1.00e+00], 1:0.258 s, tol=2.54e-01  68.00 % -> 2.55 s
  23 : 0.700 mm, [0.00e+00; 1.00e+00], 1:0.253 s, tol=2.53e-01  70.00 % -> 2.38 s
  24 : 0.720 mm, [0.00e+00; 1.00e+00], 1:0.252 s, tol=2.02e-01  72.00 % -> 2.25 s
  25 : 0.740 mm, [0.00e+00; 1.00e+00], 1:0.256 s, tol=1.76e-01  74.00 % -> 2.16 s
  26 : 0.760 mm, [0.00e+00; 1.00e+00], 1:0.256 s, tol=1.87e-01  76.00 % -> 2.02 s
  27 : 0.780 mm, [0.00e+00; 1.00e+00], 1:0.269 s, tol=1.72e-01  78.00 % -> 1.97 s
  28 : 0.800 mm, [0.00e+00; 1.00e+00], 1:0.258 s, tol=1.65e-01  80.00 % -> 1.74 s
  29 : 0.820 mm, [0.00e+00; 1.00e+00], 1:0.260 s, tol=1.90e-01  82.00 % -> 1.60 s
  30 : 0.840 mm, [0.00e+00; 1.00e+00], 1:0.262 s, tol=1.65e-01  84.00 % -> 1.45 s
  31 : 0.860 mm, [0.00e+00; 1.00e+00], 1:0.274 s, tol=1.51e-01  86.00 % -> 1.34 s
  32 : 0.880 mm, [0.00e+00; 1.00e+00], 1:0.264 s, tol=1.15e-01  88.00 % -> 1.11 s
  33 : 0.900 mm, [0.00e+00; 1.00e+00], 1:0.257 s, tol=1.10e-01  90.00 % -> 912.68 ms
  34 : 0.920 mm, [0.00e+00; 1.00e+00], 1:0.266 s, tol=9.93e-02  92.00 % -> 762.97 ms
  35 : 0.940 mm, [0.00e+00; 1.00e+00], 1:0.261 s, tol=9.94e-02  94.00 % -> 567.12 ms
  36 : 0.960 mm, [0.00e+00; 1.00e+00], 1:0.180 s, tol=9.12e-02  96.00 % -> 262.38 ms
  37 : 0.980 mm, [0.00e+00; 1.00e+00], 1:0.264 s, tol=1.02e-01  98.00 % -> 193.92 ms
  38 : 1.000 mm, [0.00e+00; 1.00e+00], 1:0.274 s, tol=9.35e-02  100.00 % -> -0.00 µs
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_CP_optimMesh nL=50/force-displacement.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_CP_optimMesh nL=50/simulation.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_CP_optimMesh nL=50/summary.txt
Loaded:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.1/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_CP_optimMesh nL=50/force-displacement.pickle

Generate movie 01/38 (2.63 %) 8.38 s
Generate movie 02/38 (5.26 %) 7.03 s
Generate movie 03/38 (7.89 %) 6.96 s
Generate movie 04/38 (10.53 %) 6.74 s
Generate movie 05/38 (13.16 %) 6.51 s
Generate movie 06/38 (15.79 %) 4.89 s
Generate movie 07/38 (18.42 %) 6.06 s
Generate movie 08/38 (21.05 %) 5.95 s
Generate movie 09/38 (23.68 %) 5.72 s
Generate movie 10/38 (26.32 %) 5.55 s
Generate movie 11/38 (28.95 %) 5.32 s
Generate movie 12/38 (31.58 %) 5.12 s
Generate movie 13/38 (34.21 %) 4.92 s
Generate movie 14/38 (36.84 %) 4.74 s
Generate movie 15/38 (39.47 %) 4.56 s
Generate movie 16/38 (42.11 %) 4.36 s
Generate movie 17/38 (44.74 %) 4.15 s
Generate movie 18/38 (47.37 %) 3.97 s
Generate movie 19/38 (50.00 %) 3.77 s
Generate movie 20/38 (52.63 %) 3.57 s
Generate movie 21/38 (55.26 %) 3.35 s
Generate movie 22/38 (57.89 %) 3.13 s
Generate movie 23/38 (60.53 %) 2.98 s
Generate movie 24/38 (63.16 %) 2.78 s
Generate movie 25/38 (65.79 %) 2.56 s
Generate movie 26/38 (68.42 %) 2.37 s
Generate movie 27/38 (71.05 %) 2.17 s
Generate movie 28/38 (73.68 %) 1.99 s
Generate movie 29/38 (76.32 %) 1.78 s
Generate movie 30/38 (78.95 %) 1.58 s
Generate movie 31/38 (81.58 %) 1.38 s
Generate movie 32/38 (84.21 %) 1.20 s
Generate movie 33/38 (86.84 %) 982.79 ms
Generate movie 34/38 (89.47 %) 789.65 ms
Generate movie 35/38 (92.11 %) 591.13 ms
Generate movie 36/38 (94.74 %) 394.39 ms
Generate movie 37/38 (97.37 %) 195.25 ms
Generate movie 38/38 (100.00 %) 0.00 µs


 12 import matplotlib.pyplot as plt
 13 import numpy as np
 14
 15 from EasyFEA import (
 16     Display,
 17     Folder,
 18     Models,
 19     Tic,
 20     ElemType,
 21     Simulations,
 22     Paraview,
 23     PyVista,
 24 )
 25 from EasyFEA.Geoms import Point, Points, Domain, Circle
 26
 27 if __name__ == "__main__":
 28     Display.Clear()
 29
 30     # ----------------------------------------------
 31     # Configuration
 32     # ----------------------------------------------
 33
 34     # simu options
 35     doSimu = True
 36     meshTest = True
 37     optimMesh = True
 38
 39     # outputs
 40     folder = Folder.Results_Dir()
 41     pltIter = False
 42     pltLoad = False
 43     makeMovie = True
 44     makeParaview = False
 45
 46     # geom
 47     dim = 2
 48     L = 250  # mm
 49     nL = 50
 50     ep = 100
 51
 52     # model
 53     E = 2e4  # MPa
 54     v = 0.18
 55
 56     split = "Miehe"
 57     regu = "AT1"
 58     Gc = 130  # J/m2
 59     Gc *= 1000 / 1e6  # mJ/mm2
 60
 61     # convergence
 62     tolConv = 1e-0
 63     convOption = 2
 64
 65     # load
 66     uMax = 1  # mm
 67     inc0 = uMax / 25
 68     inc1 = inc0 / 2
 69
 70     config = f"""
 71     uMax = {uMax}
 72
 73     inc0 = {inc0}
 74     inc1 = {inc1}
 75
 76     while ud <= uMax:
 77
 78     ud += inc0 if simu.damage.max() < 0.6 else inc1
 79
 80     simu.add_dirichlet(nodes_circle, [0], ['d'], "damage")
 81     simu.add_dirichlet(nodes_y0, [0]*dim, simu.Get_dofs())
 82     simu.add_dirichlet(nodes_load, [ud], ['y'])
 83     """
 84
 85     # ----------------------------------------------
 86     # Mesh
 87     # ----------------------------------------------
 88     l0 = L / nL
 89
 90     if meshTest:
 91         hC = l0
 92     else:
 93         hC = 0.5
 94         # hC = 0.25
 95
 96     p1 = Point()
 97     p2 = Point(L, 0)
 98     p3 = Point(L, L)
 99     p4 = Point(2 * L - 30, L)
100     p5 = Point(2 * L, L)
101     p6 = Point(2 * L, 2 * L)
102     p7 = Point(0, 2 * L)
103     if optimMesh:
104         # hauteur zone rafinée
105         h = 100
106         refineDomain = Domain(Point(0, L - h / 3), Point(L + h / 3, L + h), hC)
107         hD = hC * 5
108     else:
109         refineDomain = None
110         hD = hC
111
112     contour = Points([p1, p2, p3, p4, p5, p6, p7], hD)
113
114     circle = Circle(p5, 100)
115
116     if dim == 2:
117         mesh = contour.Mesh_2D([], ElemType.TRI3, refineGeoms=[refineDomain])
118     else:
119         mesh = contour.Mesh_Extrude(
120             [], [0, 0, -ep], [3], ElemType.HEXA8, refineGeoms=[refineDomain]
121         )
122
123     # Display.Plot_Mesh(mesh)
124     # Display.Plot_Tags(mesh)
125     # from EasyFEA import PyVista
126     # PyVista.Plot_Mesh(mesh).show()
127
128     nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
129     nodes_load = mesh.Nodes_Conditions(lambda x, y, z: (y == L) & (x >= 2 * L - 30))
130     node3 = mesh.Nodes_Point(p3)
131     node4 = mesh.Nodes_Point(p4)
132     nodes_circle = mesh.Nodes_Cylinder(circle, direction=[0, 0, ep])
133     nodes_edges = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) | (y == 0))
134
135     # ----------------------------------------------
136     # Simulation
137     # ----------------------------------------------
138     material = Models.Elastic.Isotropic(dim, E, v, True, ep)
139     pfm = Models.PhaseField(material, split, regu, Gc, l0)
140
141     folder_save = Simulations.PhaseField.Folder(
142         f"{folder}{dim}D",
143         "",
144         pfm.split,
145         pfm.regularization,
146         "CP",
147         tolConv,
148         "",
149         meshTest,
150         optimMesh,
151         nL=nL,
152     )
153
154     Display.MyPrint(folder_save, "green")
155
156     if doSimu:
157         simu = Simulations.PhaseField(mesh, pfm)
158         simu.Results_Set_Bc_Summary(config)
159
160         dofsY_load = simu.Bc_dofs_nodes(nodes_load, ["y"])
161
162         if pltIter:
163             axIter = Display.Plot_Result(simu, "damage")
164
165             axLoad = Display.Init_Axes()
166             axLoad.set_xlabel("displacement [mm]")
167             axLoad.set_ylabel("load [kN]")
168
169         displacement = []
170         force = []
171         ud = -inc0
172         iter = -1
173
174         while ud <= uMax:
175             # update displacement
176             iter += 1
177             ud += inc0 if simu.damage.max() < 0.6 else inc1
178
179             # update boundary conditions
180             simu.Bc_Init()
181             simu.add_dirichlet(nodes_circle, [0], ["d"], "damage")
182             simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns())
183             simu.add_dirichlet(nodes_load, [ud], ["y"])
184
185             # solve
186             u, d, Ku, convergence = simu.Solve(tolConv, 500, convOption)
187             # calc load
188             fr = np.sum(Ku[dofsY_load, :] @ u)
189
190             # saves load and displacement
191             displacement.append(ud)
192             force.append(fr)
193
194             # print iter
195             simu.Results_Set_Iteration_Summary(iter, ud, "mm", ud / uMax, True)
196
197             # saves iteration
198             simu.Save_Iter()
199
200             if pltIter:
201                 plt.figure(axIter.figure)
202                 Display.Plot_Result(simu, "damage", ax=axIter)
203                 plt.pause(1e-12)
204
205                 plt.figure(axLoad.figure)
206                 axLoad.scatter(ud, fr / 1000, c="black")
207                 plt.pause(1e-12)
208
209             if not convergence or np.max(d[nodes_edges]) >= 1:
210                 # stops simulation if damage occurs on edges or convergence has not been reached
211                 break
212
213         # saves load and displacement
214         displacement = np.asarray(displacement)
215         force = np.asarray(force)
216         Simulations.Save_pickle(
217             (force, displacement), folder_save, "force-displacement"
218         )
219
220         # saves the simulation
221         simu.Save(folder_save)
222
223     else:
224         simu = Simulations.Load_Simu(folder_save)
225         mesh = simu.mesh
226
227     force, displacement = Simulations.Load_pickle(folder_save, "force-displacement")
228
229     # ----------------------------------------------
230     # Results
231     # ----------------------------------------------
232     Display.Plot_Result(simu, "damage", folder=folder_save)
233     Display.Plot_Mesh(simu)
234     Display.Plot_BoundaryConditions(simu)
235
236     axLoad = Display.Init_Axes()
237     axLoad.set_xlabel("displacement [mm]")
238     axLoad.set_ylabel("load [kN]")
239     axLoad.plot(displacement, force / 1000, c="blue")
240     Display.Save_fig(folder_save, "forcedep")
241
242     Display.Plot_Iter_Summary(simu, folder_save)
243
244     if makeMovie:
245         simu.Set_Iter(-1)
246         depMax = simu.Result("displacement_norm").max()
247         deformFactor = L * 0.05 / depMax
248
249         iterations = np.arange(0, simu.Niter, simu.Niter // 20)
250
251         def Func(plotter, iter):
252             simu.Set_Iter(iterations[iter])
253
254             grid = PyVista._pvMesh(simu, "damage", deformFactor)
255
256             tresh = grid.threshold((0, 0.8))
257
258             PyVista.Plot(
259                 tresh,
260                 "damage",
261                 deformFactor,
262                 plotMesh=True,
263                 plotter=plotter,
264                 clim=(0, 1),
265             )
266
267         PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
268
269     if makeParaview:
270         Paraview.Save_simu(simu, folder_save)
271
272     if doSimu:
273         Tic.Plot_History(folder_save, False)
274
275     plt.show()

Total running time of the script: (0 minutes 22.436 seconds)

Gallery generated by Sphinx-Gallery