.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/PhaseField/PlateWithHole.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_PlateWithHole.py: PlateWithHole ============= Damage simulation for a plate with a hole subjected to compression. .. GENERATED FROM PYTHON SOURCE LINES 11-369 .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_001.gif :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_001.gif :class: sphx-glr-single-img .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_002.png :alt: $\phi$ :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_003.png :alt: TRI3: Ne = 12578, Nn = 6444 :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_004.png :alt: Boundary conditions :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_005.png :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_006.png :alt: PlateWithHole :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_PlateWithHole_007.png :alt: Summary :srcset: /examples/PhaseField/images/sphx_glr_PlateWithHole_007.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh 2 : 0.000 μm, [0.00e+00; 0.00e+00], 1:0.867 s, tol=0.00e+00 3 : 0.800 μm, [0.00e+00; 0.00e+00], 1:0.834 s, tol=1.00e+00 3.20 % -> 50.45 s 4 : 1.600 μm, [1.03e-07; 3.70e-02], 1:0.840 s, tol=8.04e-01 6.40 % -> 36.86 s 5 : 2.400 μm, [3.74e-05; 3.70e-02], 1:0.837 s, tol=5.54e-01 9.60 % -> 31.52 s 6 : 3.200 μm, [8.41e-05; 3.70e-02], 1:0.845 s, tol=4.37e-01 12.80 % -> 28.77 s 7 : 4.000 μm, [1.50e-04; 3.70e-02], 1:0.830 s, tol=3.60e-01 16.00 % -> 26.14 s 8 : 4.800 μm, [2.34e-04; 3.70e-02], 1:0.836 s, tol=3.05e-01 19.20 % -> 24.62 s 9 : 5.600 μm, [3.37e-04; 3.70e-02], 1:0.839 s, tol=2.65e-01 22.40 % -> 23.26 s 10 : 6.400 μm, [4.60e-04; 3.71e-02], 1:0.836 s, tol=2.34e-01 25.60 % -> 21.86 s 11 : 7.200 μm, [6.01e-04; 4.57e-02], 1:0.831 s, tol=2.10e-01 28.80 % -> 20.56 s 12 : 8.000 μm, [7.63e-04; 5.81e-02], 1:0.838 s, tol=1.90e-01 32.00 % -> 19.59 s 13 : 8.800 μm, [9.46e-04; 7.21e-02], 1:0.843 s, tol=1.73e-01 35.20 % -> 18.61 s 14 : 9.600 μm, [1.15e-03; 8.78e-02], 1:0.834 s, tol=1.59e-01 38.40 % -> 17.40 s 15 : 10.400 μm, [1.38e-03; 1.05e-01], 1:0.837 s, tol=1.47e-01 41.60 % -> 16.45 s 16 : 11.200 μm, [1.62e-03; 1.25e-01], 1:0.837 s, tol=1.37e-01 44.80 % -> 15.47 s 17 : 12.000 μm, [1.89e-03; 1.46e-01], 1:0.833 s, tol=1.28e-01 48.00 % -> 14.43 s 18 : 12.800 μm, [2.19e-03; 1.70e-01], 1:0.842 s, tol=1.21e-01 51.20 % -> 13.65 s 19 : 13.600 μm, [2.50e-03; 1.97e-01], 1:0.832 s, tol=1.14e-01 54.40 % -> 12.56 s 20 : 14.400 μm, [2.85e-03; 2.29e-01], 1:0.836 s, tol=1.07e-01 57.60 % -> 11.69 s 21 : 15.200 μm, [3.21e-03; 2.65e-01], 1:0.840 s, tol=1.02e-01 60.80 % -> 10.83 s 22 : 16.000 μm, [3.61e-03; 3.04e-01], 1:0.837 s, tol=9.68e-02 64.00 % -> 9.88 s 23 : 16.800 μm, [4.03e-03; 3.48e-01], 1:0.830 s, tol=9.23e-02 67.20 % -> 8.91 s 24 : 17.600 μm, [4.47e-03; 3.96e-01], 1:0.833 s, tol=8.82e-02 70.40 % -> 8.06 s 25 : 18.400 μm, [4.95e-03; 4.49e-01], 1:0.842 s, tol=8.44e-02 73.60 % -> 7.25 s 26 : 19.200 μm, [5.45e-03; 5.05e-01], 1:0.833 s, tol=8.10e-02 76.80 % -> 6.29 s 27 : 20.000 μm, [5.98e-03; 5.61e-01], 1:0.843 s, tol=7.80e-02 80.00 % -> 5.48 s 28 : 20.800 μm, [6.53e-03; 6.38e-01], 1:0.836 s, tol=7.53e-02 83.20 % -> 4.56 s 29 : 21.000 μm, [7.11e-03; 7.36e-01], 1:0.840 s, tol=2.04e-02 84.00 % -> 4.48 s 30 : 21.200 μm, [7.30e-03; 8.35e-01], 1:0.843 s, tol=2.10e-02 84.80 % -> 4.38 s 31 : 21.400 μm, [7.46e-03; 9.20e-01], 1:0.839 s, tol=2.34e-02 85.60 % -> 4.23 s 32 : 21.600 μm, [7.64e-03; 9.75e-01], 1:0.839 s, tol=2.99e-02 86.40 % -> 4.09 s 33 : 21.800 μm, [7.87e-03; 9.98e-01], 1:0.844 s, tol=4.11e-02 87.20 % -> 3.96 s 34 : 22.000 μm, [8.18e-03; 1.00e+00], 1:0.832 s, tol=5.57e-02 88.00 % -> 3.75 s 35 : 22.200 μm, [8.57e-03; 1.00e+00], 1:0.837 s, tol=6.35e-02 88.80 % -> 3.59 s 36 : 22.400 μm, [8.97e-03; 1.00e+00], 1:0.834 s, tol=8.48e-02 89.60 % -> 3.39 s 37 : 22.600 μm, [9.40e-03; 1.00e+00], 1:0.841 s, tol=8.07e-02 90.40 % -> 3.21 s 38 : 22.800 μm, [9.82e-03; 1.00e+00], 1:0.833 s, tol=7.54e-02 91.20 % -> 2.98 s 39 : 23.000 μm, [1.01e-02; 1.00e+00], 1:0.842 s, tol=7.34e-02 92.00 % -> 2.78 s 40 : 23.200 μm, [1.04e-02; 1.00e+00], 1:0.839 s, tol=7.26e-02 92.80 % -> 2.54 s 41 : 23.400 μm, [1.06e-02; 1.00e+00], 1:0.834 s, tol=7.46e-02 93.60 % -> 2.28 s 42 : 23.600 μm, [1.08e-02; 1.00e+00], 1:0.841 s, tol=7.51e-02 94.40 % -> 2.05 s 43 : 23.800 μm, [1.10e-02; 1.00e+00], 1:0.837 s, tol=7.56e-02 95.20 % -> 1.77 s 44 : 24.000 μm, [1.11e-02; 1.00e+00], 1:0.836 s, tol=7.53e-02 96.00 % -> 1.50 s 45 : 24.200 μm, [1.12e-02; 1.00e+00], 1:0.838 s, tol=7.52e-02 96.80 % -> 1.22 s 46 : 24.400 μm, [1.14e-02; 1.00e+00], 1:0.835 s, tol=7.46e-02 97.60 % -> 924.04 ms 47 : 24.600 μm, [1.15e-02; 1.00e+00], 1:0.834 s, tol=7.42e-02 98.40 % -> 623.52 ms 48 : 24.800 μm, [1.16e-02; 1.00e+00], 1:0.837 s, tol=7.27e-02 99.20 % -> 317.14 ms 49 : 25.000 μm, [1.17e-02; 1.00e+00], 1:0.836 s, tol=7.13e-02 100.00 % -> -0.00 µs Saved: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/force-displacement.pickle Saved: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/simulation.pickle Saved: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/PlateWithHole/Test/Elas_Isot_Miehe_AT2_DP_optimMesh/summary.txt Generate movie 01/24 (4.17 %) 2.64 s Generate movie 02/24 (8.33 %) 2.34 s Generate movie 03/24 (12.50 %) 2.24 s Generate movie 04/24 (16.67 %) 2.06 s Generate movie 05/24 (20.83 %) 1.94 s Generate movie 06/24 (25.00 %) 1.85 s Generate movie 07/24 (29.17 %) 1.80 s Generate movie 08/24 (33.33 %) 1.64 s Generate movie 09/24 (37.50 %) 1.52 s Generate movie 10/24 (41.67 %) 1.44 s Generate movie 11/24 (45.83 %) 1.32 s Generate movie 12/24 (50.00 %) 1.24 s Generate movie 13/24 (54.17 %) 1.12 s Generate movie 14/24 (58.33 %) 1.02 s Generate movie 15/24 (62.50 %) 929.11 ms Generate movie 16/24 (66.67 %) 815.38 ms Generate movie 17/24 (70.83 %) 725.39 ms Generate movie 18/24 (75.00 %) 606.98 ms Generate movie 19/24 (79.17 %) 514.75 ms Generate movie 20/24 (83.33 %) 407.71 ms Generate movie 21/24 (87.50 %) 313.89 ms Generate movie 22/24 (91.67 %) 206.79 ms Generate movie 23/24 (95.83 %) 104.60 ms Generate movie 24/24 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 12 from EasyFEA import ( Display, Folder, Models, plt, np, Tic, ElemType, Mesh, Simulations, PyVista, Paraview, ) from EasyFEA.Geoms import Domain, Circle import multiprocessing # Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- # simu options doSimu = True meshTest = True optimMesh = True useParallel = False nProcs = 4 # number of processes in parallel # outputs folder = Folder.Results_Dir() plotMesh = False plotIter = False plotResult = True plotEnergy = False showFig = True saveParaview = False makeMovie = True # models # splits = ["Bourdin","Amor","Miehe","Stress"] # Splits Isotropes # splits = ["He","AnisotStrain","AnisotStress","Zhang"] # Splits Anisotropes # splits = ["Bourdin","Amor","Miehe","Stress","He","AnisotStrain","AnisotStress","Zhang"] # splits = ["Zhang"] # splits = ["AnisotStrain","AnisotStress","Zhang"] splits = ["Miehe"] regus = ["AT2"] # ["AT1", "AT2"] # regus = ["AT1", "AT2"] l0 = 0.12e-3 # convergence solver = ( Models.PhaseField.SolverType.History ) # ["History", "HistoryDamage", "BoundConstrain"] maxIter = 1000 tolConv = 1e-0 # ---------------------------------------------- # Mesh # ---------------------------------------------- def DoMesh( L: float, h: float, diam: float, thickness: float, l0: float, split: str ) -> Mesh: clC = l0 * 2 if meshTest else l0 / 2 if optimMesh: clD = l0 * 4 refineZone = diam * 1.5 / 2 if split in ["Bourdin", "Amor"]: refineGeom = Domain((0, h / 2 - refineZone), (L, h / 2 + refineZone), clC) else: refineGeom = Domain((L / 2 - refineZone, 0), (L / 2 + refineZone, h), clC) else: clD = l0 if meshTest else l0 / 2 refineGeom = None domain = Domain((0, 0), (L, h), clD) circle = Circle((L / 2, h / 2), diam, clD, isHollow=True) # ax = Display.Init_Axes() # domain.Plot(ax, color="k", plotPoints=False) # circle.Plot(ax, color="k", plotPoints=False) # # if refineGeom != None: # # refineGeom.Plot(ax, color='k', plotPoints=False) # # ax.scatter(((L+diam)/2, L/2), (h/2, (h+diam)/2), c='k') # ax.axis("off") # Display.Save_fig(folder, "sample", True) mesh = domain.Mesh_2D([circle], ElemType.TRI3, refineGeoms=[refineGeom]) # ax = Display.Plot_Mesh(mesh, lw=0.3, facecolors="white") # ax.axis("off") # ax.set_title("") # Display.Save_fig(folder, "mesh", transparent=True) return mesh # ---------------------------------------------- # Do Simu # ---------------------------------------------- def DoSimu(split: str, regu: str): folder_save = Simulations.PhaseField.Folder( folder, "Elas_Isot", split, regu, "DP", tolConv, solver, meshTest, optimMesh ) Display.MyPrint(folder_save, "green") # ---------------------------------------------- # Geom # ---------------------------------------------- L = 15e-3 h = 30e-3 thickness = 1 diam = 6e-3 # load units unitU = "μm" unitF = "kN/mm" unit = 1e6 # ---------------------------------------------- # Mesh # ---------------------------------------------- if doSimu: mesh = DoMesh(L, h, diam, thickness, l0, split) # Get Nodes nodes_lower = mesh.Nodes_Conditions(lambda x, y, z: y == 0) nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h) nodes_x0y0 = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (y == 0)) nodes_edges = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"]) nodes_upper = mesh.Nodes_Conditions(lambda x, y, z: y == h) # ---------------------------------------------- # Boundary conditions # ---------------------------------------------- threshold = 0.6 u_max = 25e-6 # u_max = 35e-6 uinc0 = 8e-7 if meshTest else 8e-8 uinc1 = 2e-7 if meshTest else 2e-8 config = f""" while ud <= u_max: ud += uinc0 if simu.damage.max() < threshold else uinc1 u_max = {u_max} uinc0 = {uinc0:.1e} (simu.damage.max() < {threshold}) uinc1 = {uinc1:.1e} simu.add_dirichlet(nodes_lower, [0], ["y"]) simu.add_dirichlet(nodes_x0y0, [0], ["x"]) simu.add_dirichlet(nodes_upper, [-ud], ["y"]) if dim == 3: simu.add_dirichlet(nodes_y0z0, [0], ["z"]) """ # ---------------------------------------------- # Material # ---------------------------------------------- E = 12e9 v = 0.3 planeStress = False material = Models.Elastic.Isotropic(2, E, v, planeStress, thickness) gc = 1.4 pfm = Models.PhaseField(material, split, regu, gc, l0, solver=solver) # ---------------------------------------------- # Simulation # ---------------------------------------------- simu = Simulations.PhaseField(mesh, pfm, verbosity=False) simu.Results_Set_Bc_Summary(config) dofsY_upper = simu.Bc_dofs_nodes(nodes_upper, ["y"]) def Apply_BC(ud: float): simu.Bc_Init() simu.add_dirichlet(nodes_lower, [0], ["y"]) simu.add_dirichlet(nodes_x0y0, [0], ["x"]) simu.add_dirichlet(nodes_upper, [-ud], ["y"]) # INIT displacement = [] force = [] ud = -uinc0 iter = 0 nDetect = 0 if plotIter: axIter = Display.Plot_Result(simu, "damage", nodeValues=True) force = np.asarray(force) displacement = np.asarray(displacement) _, axLoad = Display.Plot_Force_Displacement( force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]" ) while ud <= u_max: iter += 1 ud += uinc0 if simu.damage.max() < threshold else uinc1 Apply_BC(ud) u, d, Ku, convergence = simu.Solve(tolConv, maxIter) simu.Save_Iter() # stop if the simulation does not converge if not convergence: break f = np.sum(Ku[dofsY_upper, :] @ u) simu.Results_Set_Iteration_Summary(iter, ud * unit, unitU, ud / u_max, True) # Detect damaged edges if np.any(d[nodes_edges] >= 1): nDetect += 1 if nDetect == 10: break displacement = np.concatenate((displacement, [ud])) force = np.concatenate((force, [f])) if plotIter: Display.Plot_Result(simu, "damage", nodeValues=True, ax=axIter) plt.figure(axIter.figure) plt.pause(1e-12) force = np.asarray(force) displacement = np.asarray(displacement) Display.Plot_Force_Displacement( force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]", ax=axLoad, )[1] plt.figure(axLoad.figure) plt.pause(1e-12) # ---------------------------------------------- # Saving # ---------------------------------------------- force = np.asarray(force) displacement = np.asarray(displacement) print() Simulations.Save_pickle( (force, displacement), folder_save, "force-displacement" ) simu.Save(folder_save) else: # ---------------------------------------------- # Load # ---------------------------------------------- simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save) force, displacement = Simulations.Load_pickle(folder_save, "force-displacement") # ---------------------------------------------- # Results # --------------------------------------------- if plotEnergy: Display.Plot_Energy(simu, force, displacement, N=400, folder=folder_save) if plotResult: Display.Plot_Result( simu, "damage", nodeValues=True, colorbarIsClose=True, folder=folder_save, filename="damage", ) Display.Plot_Mesh(simu) Display.Plot_BoundaryConditions(simu) Display.Plot_Iter_Summary(simu, folder_save, None, None) Display.Plot_Force_Displacement( force / unit, displacement * unit, f"ud [{unitU}]", f"f [{unitF}]", folder_save, ) if plotMesh: Display.Plot_Mesh(mesh) if saveParaview: Paraview.Save_simu(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 doSimu: Tic.Plot_History(folder_save, details=False) if showFig: plt.show() # plt.close("all") Tic.Clear() if __name__ == "__main__": # generates configs Splits = [] Regus = [] for split in splits.copy(): for regu in regus.copy(): Splits.append(split) Regus.append(regu) if useParallel: items = [(split, regu) for split, regu in zip(Splits, Regus)] with multiprocessing.Pool(nProcs) as pool: for result in pool.starmap(DoSimu, items): pass else: [DoSimu(split, regu) for split, regu in zip(Splits, Regus)] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 45.798 seconds) .. _sphx_glr_download_examples_PhaseField_PlateWithHole.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: PlateWithHole.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: PlateWithHole.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: PlateWithHole.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_