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
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/latest/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50
   1: 0.000mm, [0.00e+00;0.00e+00], 1:0.062s, tol=0.00e+00
   2: 0.040mm, [0.00e+00;0.00e+00], 1:0.049s, tol=1.00e+00  4.00% -> 1.17s
   3: 0.080mm, [0.00e+00;0.00e+00], 1:0.049s, tol=7.50e-01  8.00% -> 1.13s
   4: 0.120mm, [0.00e+00;0.00e+00], 1:0.048s, tol=5.56e-01  12.00% -> 1.06s
   5: 0.160mm, [0.00e+00;0.00e+00], 1:0.049s, tol=4.38e-01  16.00% -> 1.04s
   6: 0.200mm, [0.00e+00;0.00e+00], 1:0.049s, tol=3.60e-01  20.00% -> 980.68ms
   7: 0.240mm, [0.00e+00;0.00e+00], 1:0.048s, tol=3.06e-01  24.00% -> 920.70ms
   8: 0.280mm, [0.00e+00;8.98e-03], 1:0.049s, tol=2.65e-01  28.00% -> 887.61ms
   9: 0.320mm, [0.00e+00;4.57e-02], 1:0.049s, tol=2.32e-01  32.00% -> 838.01ms
  10: 0.360mm, [0.00e+00;1.10e-01], 1:0.049s, tol=2.07e-01  36.00% -> 791.09ms
  11: 0.400mm, [0.00e+00;2.07e-01], 1:0.049s, tol=1.86e-01  40.00% -> 733.63ms
  12: 0.440mm, [0.00e+00;3.48e-01], 1:0.050s, tol=1.69e-01  44.00% -> 694.23ms
  13: 0.480mm, [0.00e+00;5.46e-01], 1:0.050s, tol=1.55e-01  48.00% -> 646.57ms
  14: 0.520mm, [0.00e+00;7.78e-01], 1:0.050s, tol=1.62e-01  52.00% -> 600.70ms
  15: 0.540mm, [0.00e+00;9.48e-01], 1:0.050s, tol=2.22e-01  54.00% -> 601.02ms
  16: 0.560mm, [0.00e+00;9.93e-01], 1:0.050s, tol=3.80e-01  56.00% -> 592.89ms
  17: 0.580mm, [0.00e+00;1.00e+00], 1:0.049s, tol=4.44e-01  58.00% -> 567.60ms
  18: 0.600mm, [0.00e+00;1.00e+00], 1:0.050s, tol=3.82e-01  60.00% -> 563.71ms
  19: 0.620mm, [0.00e+00;1.00e+00], 1:0.049s, tol=3.28e-01  62.00% -> 543.41ms
  20: 0.640mm, [0.00e+00;1.00e+00], 1:0.050s, tol=2.91e-01  64.00% -> 532.69ms
  21: 0.660mm, [0.00e+00;1.00e+00], 1:0.050s, tol=2.78e-01  66.00% -> 513.37ms
  22: 0.680mm, [0.00e+00;1.00e+00], 1:0.050s, tol=2.54e-01  68.00% -> 490.53ms
  23: 0.700mm, [0.00e+00;1.00e+00], 1:0.049s, tol=2.53e-01  70.00% -> 461.56ms
  24: 0.720mm, [0.00e+00;1.00e+00], 1:0.049s, tol=2.03e-01  72.00% -> 438.95ms
  25: 0.740mm, [0.00e+00;1.00e+00], 1:0.049s, tol=1.76e-01  74.00% -> 415.51ms
  26: 0.760mm, [0.00e+00;1.00e+00], 1:0.049s, tol=1.87e-01  76.00% -> 387.50ms
  27: 0.780mm, [0.00e+00;1.00e+00], 1:0.050s, tol=1.72e-01  78.00% -> 368.26ms
  28: 0.800mm, [0.00e+00;1.00e+00], 1:0.050s, tol=1.65e-01  80.00% -> 340.08ms
  29: 0.820mm, [0.00e+00;1.00e+00], 1:0.050s, tol=1.90e-01  82.00% -> 307.76ms
  30: 0.840mm, [0.00e+00;1.00e+00], 1:0.051s, tol=1.64e-01  84.00% -> 279.13ms
  31: 0.860mm, [0.00e+00;1.00e+00], 1:0.050s, tol=1.51e-01  86.00% -> 244.52ms
  32: 0.880mm, [0.00e+00;1.00e+00], 1:0.049s, tol=1.15e-01  88.00% -> 207.88ms
  33: 0.900mm, [0.00e+00;1.00e+00], 1:0.049s, tol=1.10e-01  90.00% -> 175.68ms
  34: 0.920mm, [0.00e+00;1.00e+00], 1:0.050s, tol=9.94e-02  92.00% -> 142.55ms
  35: 0.940mm, [0.00e+00;1.00e+00], 1:0.049s, tol=9.95e-02  94.00% -> 106.26ms
  36: 0.960mm, [0.00e+00;1.00e+00], 1:0.049s, tol=9.12e-02  96.00% -> 72.10ms
  37: 0.980mm, [0.00e+00;1.00e+00], 1:0.049s, tol=1.02e-01  98.00% -> 36.28ms
  38: 1.000mm, [0.00e+00;1.00e+00], 1:0.050s, tol=9.35e-02  100.00% -> -0.00µs
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/latest/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/latest/examples/PhaseField/results/LShape2D/Test/Miehe_AT1_optimMesh nL=50

Generate movie 01/38 (2.63 %) 5.32 s
Generate movie 02/38 (5.26 %) 4.87 s
Generate movie 03/38 (7.89 %) 4.72 s
Generate movie 04/38 (10.53 %) 4.56 s
Generate movie 05/38 (13.16 %) 4.48 s
Generate movie 06/38 (15.79 %) 4.31 s
Generate movie 07/38 (18.42 %) 4.16 s
Generate movie 08/38 (21.05 %) 4.02 s
Generate movie 09/38 (23.68 %) 3.91 s
Generate movie 10/38 (26.32 %) 3.76 s
Generate movie 11/38 (28.95 %) 3.62 s
Generate movie 12/38 (31.58 %) 3.53 s
Generate movie 13/38 (34.21 %) 3.42 s
Generate movie 14/38 (36.84 %) 3.23 s
Generate movie 15/38 (39.47 %) 3.09 s
Generate movie 16/38 (42.11 %) 3.01 s
Generate movie 17/38 (44.74 %) 2.85 s
Generate movie 18/38 (47.37 %) 2.68 s
Generate movie 19/38 (50.00 %) 2.57 s
Generate movie 20/38 (52.63 %) 2.41 s
Generate movie 21/38 (55.26 %) 2.32 s
Generate movie 22/38 (57.89 %) 2.17 s
Generate movie 23/38 (60.53 %) 2.03 s
Generate movie 24/38 (63.16 %) 1.90 s
Generate movie 25/38 (65.79 %) 1.76 s
Generate movie 26/38 (68.42 %) 1.63 s
Generate movie 27/38 (71.05 %) 1.52 s
Generate movie 28/38 (73.68 %) 1.36 s
Generate movie 29/38 (76.32 %) 1.24 s
Generate movie 30/38 (78.95 %) 1.08 s
Generate movie 31/38 (81.58 %) 946.68 ms
Generate movie 32/38 (84.21 %) 816.97 ms
Generate movie 33/38 (86.84 %) 675.03 ms
Generate movie 34/38 (89.47 %) 545.54 ms
Generate movie 35/38 (92.11 %) 403.74 ms
Generate movie 36/38 (94.74 %) 279.30 ms
Generate movie 37/38 (97.37 %) 135.92 ms
Generate movie 38/38 (100.00 %) 0.00 µs


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

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

Gallery generated by Sphinx-Gallery