.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/PhaseField/LShape.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_PhaseField_LShape.py: LShape ====== Damage simulation for a L-part. .. GENERATED FROM PYTHON SOURCE LINES 11-287 .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_001.gif :alt: LShape :srcset: /examples/PhaseField/images/sphx_glr_LShape_001.gif :class: sphx-glr-single-img .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_002.png :alt: $\phi$ :srcset: /examples/PhaseField/images/sphx_glr_LShape_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_003.png :alt: TRI3: Ne = 4234, Nn = 2178 :srcset: /examples/PhaseField/images/sphx_glr_LShape_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_004.png :alt: Boundary conditions :srcset: /examples/PhaseField/images/sphx_glr_LShape_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_005.png :alt: LShape :srcset: /examples/PhaseField/images/sphx_glr_LShape_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_006.png :alt: LShape :srcset: /examples/PhaseField/images/sphx_glr_LShape_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_LShape_007.png :alt: Summary :srcset: /examples/PhaseField/images/sphx_glr_LShape_007.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /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 | .. code-block:: Python :lineno-start: 12 from EasyFEA import ( Display, Folder, Models, plt, np, Tic, ElemType, Simulations, Paraview, PyVista, ) from EasyFEA.Geoms import Point, Points, Domain, Circle if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- # simu options doSimu = True meshTest = True optimMesh = True # outputs pltIter = False pltLoad = False makeMovie = True makeParaview = False # geom dim = 2 L = 250 # mm nL = 50 ep = 100 # model E = 2e4 # MPa v = 0.18 split = "Miehe" regu = "AT1" Gc = 130 # J/m2 Gc *= 1000 / 1e6 # mJ/mm2 # convergence tolConv = 1e-0 convOption = 2 # load uMax = 1 # mm inc0 = uMax / 25 inc1 = inc0 / 2 config = f""" uMax = {uMax} inc0 = {inc0} inc1 = {inc1} while ud <= uMax: ud += inc0 if simu.damage.max() < 0.6 else inc1 simu.add_dirichlet(nodes_circle, [0], ['d'], "damage") simu.add_dirichlet(nodes_y0, [0]*dim, simu.Get_dofs()) simu.add_dirichlet(nodes_load, [ud], ['y']) """ folder = Folder.PhaseField_Folder( f"LShape_{dim}D", "Isot", split, regu, "", tolConv, "", meshTest, optimMesh, nL=nL, ) # ---------------------------------------------- # Mesh # ---------------------------------------------- l0 = L / nL if meshTest: hC = l0 else: hC = 0.5 # hC = 0.25 p1 = Point() p2 = Point(L, 0) p3 = Point(L, L) p4 = Point(2 * L - 30, L) p5 = Point(2 * L, L) p6 = Point(2 * L, 2 * L) p7 = Point(0, 2 * L) if optimMesh: # hauteur zone rafinée h = 100 refineDomain = Domain(Point(0, L - h / 3), Point(L + h / 3, L + h), hC) hD = hC * 5 else: refineDomain = None hD = hC contour = Points([p1, p2, p3, p4, p5, p6, p7], hD) circle = Circle(p5, 100) if dim == 2: mesh = contour.Mesh_2D([], ElemType.TRI3, refineGeoms=[refineDomain]) else: mesh = contour.Mesh_Extrude( [], [0, 0, -ep], [3], ElemType.HEXA8, refineGeoms=[refineDomain] ) # Display.Plot_Mesh(mesh) # Display.Plot_Tags(mesh) # from EasyFEA import PyVista # PyVista.Plot_Mesh(mesh).show() nodes_y0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0) nodes_load = mesh.Nodes_Conditions(lambda x, y, z: (y == L) & (x >= 2 * L - 30)) node3 = mesh.Nodes_Point(p3) node4 = mesh.Nodes_Point(p4) nodes_circle = mesh.Nodes_Cylinder(circle, [0, 0, ep]) nodes_edges = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) | (y == 0)) # ---------------------------------------------- # Simulation # ---------------------------------------------- material = Models.ElasIsot(dim, E, v, True, ep) pfm = Models.PhaseField(material, split, regu, Gc, l0) folder_save = Folder.PhaseField_Folder( folder, "", pfm.split, pfm.regularization, "CP", tolConv, "", meshTest, optimMesh, nL=nL, ) Display.MyPrint(folder_save, "green") if doSimu: simu = Simulations.PhaseFieldSimu(mesh, pfm) simu.Results_Set_Bc_Summary(config) dofsY_load = simu.Bc_dofs_nodes(nodes_load, ["y"]) if pltIter: axIter = Display.Plot_Result(simu, "damage") axLoad = Display.Init_Axes() axLoad.set_xlabel("displacement [mm]") axLoad.set_ylabel("load [kN]") displacement = [] force = [] ud = -inc0 iter = -1 while ud <= uMax: # update displacement iter += 1 ud += inc0 if simu.damage.max() < 0.6 else inc1 # update boundary conditions simu.Bc_Init() simu.add_dirichlet(nodes_circle, [0], ["d"], "damage") simu.add_dirichlet(nodes_y0, [0] * dim, simu.Get_unknowns()) simu.add_dirichlet(nodes_load, [ud], ["y"]) # solve u, d, Kglob, convergence = simu.Solve(tolConv, 500, convOption) # calc load fr = np.sum(Kglob[dofsY_load, :] @ u) # saves load and displacement displacement.append(ud) force.append(fr) # print iter simu.Results_Set_Iteration_Summary(iter, ud, "mm", ud / uMax, True) # saves iteration simu.Save_Iter() if pltIter: plt.figure(axIter.figure) Display.Plot_Result(simu, "damage", ax=axIter) plt.pause(1e-12) plt.figure(axLoad.figure) axLoad.scatter(ud, fr / 1000, c="black") plt.pause(1e-12) if not convergence or np.max(d[nodes_edges]) >= 1: # stops simulation if damage occurs on edges or convergence has not been reached break # saves load and displacement displacement = np.asarray(displacement) force = np.asarray(force) Simulations.Save_pickle( (force, displacement), folder_save, "force-displacement" ) # saves the simulation simu.Save(folder_save) else: simu = Simulations.Load_Simu(folder_save) mesh = simu.mesh force, displacement = Simulations.Load_pickle(folder_save, "force-displacement") # ---------------------------------------------- # Results # ---------------------------------------------- Display.Plot_Result(simu, "damage", folder=folder_save) Display.Plot_Mesh(simu) Display.Plot_BoundaryConditions(simu) axLoad = Display.Init_Axes() axLoad.set_xlabel("displacement [mm]") axLoad.set_ylabel("load [kN]") axLoad.plot(displacement, force / 1000, c="blue") Display.Save_fig(folder_save, "forcedep") Display.Plot_Iter_Summary(simu, folder_save) if makeMovie: simu.Set_Iter(-1) depMax = simu.Result("displacement_norm").max() deformFactor = L * 0.05 / depMax iterations = np.arange(0, simu.Niter, simu.Niter // 20) def Func(plotter, iter): simu.Set_Iter(iterations[iter]) grid = PyVista._pvMesh(simu, "damage", deformFactor) tresh = grid.threshold((0, 0.8)) PyVista.Plot( tresh, "damage", deformFactor, plotMesh=True, plotter=plotter, clim=(0, 1), ) PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif") if makeParaview: Paraview.Save_simu(simu, folder_save) if doSimu: Tic.Plot_History(folder_save, False) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 18.367 seconds) .. _sphx_glr_download_examples_PhaseField_LShape.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: LShape.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: LShape.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: LShape.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_