.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/PhaseField/CT.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_CT.py: CT == Performs damage simulation on a CT specimen. .. GENERATED FROM PYTHON SOURCE LINES 11-256 .. image-sg:: /examples/PhaseField/images/sphx_glr_CT_001.gif :alt: CT :srcset: /examples/PhaseField/images/sphx_glr_CT_001.gif :class: sphx-glr-single-img .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/PhaseField/images/sphx_glr_CT_002.png :alt: $\phi$ :srcset: /examples/PhaseField/images/sphx_glr_CT_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_CT_003.png :alt: $uy$ :srcset: /examples/PhaseField/images/sphx_glr_CT_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_CT_004.png :alt: TRI3: Ne = 8982, Nn = 4661 :srcset: /examples/PhaseField/images/sphx_glr_CT_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/PhaseField/images/sphx_glr_CT_005.png :alt: CT :srcset: /examples/PhaseField/images/sphx_glr_CT_005.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/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100 1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.276 s, tol=0.00e+00 2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.273 s, tol=1.00e+00 2.00 % -> 13.38 s 3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.279 s, tol=1.00e+00 4.00 % -> 13.40 s 4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.243 s, tol=5.56e-01 6.00 % -> 11.42 s 5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.294 s, tol=4.37e-01 8.00 % -> 13.53 s 6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.283 s, tol=3.60e-01 10.00 % -> 12.74 s 7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.278 s, tol=3.06e-01 12.00 % -> 12.24 s 8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.274 s, tol=2.65e-01 14.00 % -> 11.79 s 9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.313 s, tol=2.32e-01 16.00 % -> 13.16 s 10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.304 s, tol=2.04e-01 18.00 % -> 12.47 s 11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.240 s, tol=1.81e-01 20.00 % -> 9.60 s 12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.305 s, tol=1.65e-01 22.00 % -> 11.90 s 13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.326 s, tol=1.52e-01 24.00 % -> 12.39 s 14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.301 s, tol=1.59e-01 26.00 % -> 11.15 s 15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.304 s, tol=2.32e-01 28.00 % -> 10.94 s 16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.316 s, tol=3.72e-01 30.00 % -> 11.06 s 17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.282 s, tol=4.13e-01 32.00 % -> 9.60 s 18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.324 s, tol=3.98e-01 34.00 % -> 10.69 s 19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.299 s, tol=2.98e-01 36.00 % -> 9.57 s 20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.321 s, tol=2.73e-01 38.00 % -> 9.95 s 21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.317 s, tol=2.61e-01 40.00 % -> 9.52 s 22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.311 s, tol=2.47e-01 42.00 % -> 9.02 s 23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.312 s, tol=2.42e-01 44.00 % -> 8.74 s 24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.285 s, tol=2.07e-01 46.00 % -> 7.70 s 25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.318 s, tol=2.11e-01 48.00 % -> 8.28 s 26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.327 s, tol=1.77e-01 50.00 % -> 8.17 s 27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.316 s, tol=1.67e-01 52.00 % -> 7.59 s 28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.309 s, tol=1.74e-01 54.00 % -> 7.11 s 29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.286 s, tol=1.49e-01 56.00 % -> 6.29 s 30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.269 s, tol=1.39e-01 58.00 % -> 5.64 s 31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.204 s, tol=1.29e-01 60.00 % -> 4.08 s 32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.298 s, tol=1.31e-01 62.00 % -> 5.66 s 33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:0.302 s, tol=1.26e-01 64.00 % -> 5.44 s 34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.302 s, tol=1.12e-01 66.00 % -> 5.14 s 35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:0.317 s, tol=1.11e-01 68.00 % -> 5.08 s 36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:0.309 s, tol=1.09e-01 70.00 % -> 4.63 s 37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:0.268 s, tol=1.08e-01 72.00 % -> 3.76 s 38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:0.323 s, tol=1.09e-01 74.00 % -> 4.20 s 39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:0.318 s, tol=1.13e-01 76.00 % -> 3.82 s 40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:0.310 s, tol=1.24e-01 78.00 % -> 3.41 s 41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:0.320 s, tol=1.55e-01 80.00 % -> 3.20 s 42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:0.308 s, tol=4.20e-01 82.00 % -> 2.77 s 43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:0.311 s, tol=7.43e-01 84.00 % -> 2.49 s 44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:0.259 s, tol=1.25e-01 86.00 % -> 1.81 s 45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:0.320 s, tol=1.33e-01 88.00 % -> 1.92 s 46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:0.311 s, tol=1.92e-02 90.00 % -> 1.56 s Saved: /results/PhaseField/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100/simulation.pickle Saved: /results/PhaseField/CT_2D/Test/Isot_Miehe_AT1_optimMesh nL=100/summary.txt Generate movie 0/22 Generate movie 1/22 (4.55 %) 3.41 s Generate movie 2/22 (9.09 %) 3.27 s Generate movie 3/22 (13.64 %) 3.10 s Generate movie 4/22 (18.18 %) 2.94 s Generate movie 5/22 (22.73 %) 2.76 s Generate movie 6/22 (27.27 %) 2.62 s Generate movie 7/22 (31.82 %) 2.46 s Generate movie 8/22 (36.36 %) 2.31 s Generate movie 9/22 (40.91 %) 2.13 s Generate movie 10/22 (45.45 %) 1.97 s Generate movie 11/22 (50.00 %) 1.79 s Generate movie 12/22 (54.55 %) 1.64 s Generate movie 13/22 (59.09 %) 1.48 s Generate movie 14/22 (63.64 %) 1.31 s Generate movie 15/22 (68.18 %) 1.16 s Generate movie 16/22 (72.73 %) 959.82 ms Generate movie 17/22 (77.27 %) 814.57 ms Generate movie 18/22 (81.82 %) 647.96 ms Generate movie 19/22 (86.36 %) 490.64 ms Generate movie 20/22 (90.91 %) 327.02 ms Generate movie 21/22 (95.45 %) 166.23 ms Generate movie 22/22 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 12 from EasyFEA import ( Display, Folder, Models, plt, np, ElemType, Simulations, PyVista, Paraview, ) from EasyFEA.Geoms import Point, Points, Circle, Line, Contour, Domain if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Config # ---------------------------------------------- dim = 2 # simu options doSimu = True meshTest = True optimMesh = True # outputs plotGeom = False plotIter = False makeParaview = False makeMovie = True # geom L = 60 # mm e = 4 t = 30 r = 2 t2 = 15 diam = 8 thickness = 8 # model Gc = 100 nL = 100 l0 = L / nL split = "Miehe" regu = "AT1" folder_save = Folder.PhaseField_Folder( f"CT_{dim}D", "Isot", split, regu, "", 1, "", meshTest, optimMesh, nL=nL ) Display.MyPrint(folder_save, "green") # ---------------------------------------------- # Geom # ---------------------------------------------- clC = l0 if meshTest else l0 / 2 # meshSize on the crack clD = l0 * 2 if optimMesh else clC mS = l0 / 2 pt1 = Point(0, -L / 2) pt2 = Point(L, -L / 2) pt3 = Point(L, L / 2) pt4 = Point(0, L / 2) pt5 = Point(0, e / 2) pt6 = Point(t + r, e / 2, r=r) pt7 = Point(t + r, -e / 2, r=r) pt8 = Point(0, -e / 2) points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], clD) contour = points.Get_Contour() circle1 = Circle(Point(t2, -L / 2 + t2), diam, clD) circle2 = Circle(Point(t2, L / 2 - t2), diam, clD) if plotGeom: ax = contour.Plot() contour.Plot_Geoms([circle1, circle2], ax=ax) # ---------------------------------------------- # Mesh # ---------------------------------------------- refineGeom = ( Domain(Point(t, -e * 1.5), Point(L, e * 1.5, thickness), clC) if optimMesh else None ) if plotGeom: refineGeom.Plot(ax=ax) if dim == 2: crack = Line( Point(t + r, 0, isOpen=True), Point(t + r + 6, 0), clC, isOpen=True ) mesh = contour.Mesh_2D( [circle1, circle2], ElemType.TRI3, cracks=[crack], refineGeoms=[refineGeom], ) else: elemType = ElemType.PRISM6 pc1 = Point(t + r, 0, 0, True) pc2 = Point(t + r + 6, 0, 0) pc3 = pc2 + [0, 0, thickness] pc4 = Point(t + r, 0, thickness, True) # line1 = Line(pc1, pc2, clC, True) line2 = Line(pc2, pc3, clC, False) line3 = Line(pc3, pc4, clC, True) line4 = Line(pc4, pc1, clC, True) # crack = Contour([line1, line2, line3, line4], isOpen=True) cracks = [crack] if plotGeom: crack.Plot(ax=ax) mesh = contour.Mesh_Extrude( [circle1, circle2], [0, 0, thickness], [4], elemType, cracks=cracks, additionalLines=[line1], refineGeoms=[refineGeom], ) # PyVista.Plot_Mesh(mesh).show() # PyVista.Plot_Nodes(mesh, mesh.orphanNodes).show() nodes_1 = mesh.Nodes_Cylinder(circle1) nodes_2 = mesh.Nodes_Cylinder(circle2) nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L) # ---------------------------------------------- # Simu # ---------------------------------------------- mat = Models.ElasIsot(dim, thickness=thickness, planeStress=True) pfm = Models.PhaseField(mat, split, regu, Gc, l0) if doSimu: displacements = np.linspace(0, L / 40, 50) config = """ displacements = np.linspace(0, L/40, 50) for i, dep in enumerate(displacements): if dim == 2: simu.add_dirichlet(nodes_1, [0,-dep], ["x","y"]) simu.add_dirichlet(nodes_2, [0,dep], ["x","y"]) else: simu.add_dirichlet(nodes_1, [0,-dep, 0], ["x", "y", "z"]) simu.add_dirichlet(nodes_2, [0,dep, 0], ["x", "y", "z"]) """ simu = Simulations.PhaseFieldSimu(mesh, pfm) simu.Results_Set_Bc_Summary(config) if plotIter: ax = Display.Plot_Result(simu, "damage", 1, plotMesh=True) for i, dep in enumerate(displacements): simu.Bc_Init() if dim == 2: simu.add_dirichlet(nodes_1, [0, -dep], ["x", "y"]) simu.add_dirichlet(nodes_2, [0, dep], ["x", "y"]) else: simu.add_dirichlet(nodes_1, [0, -dep, 0], ["x", "y", "z"]) simu.add_dirichlet(nodes_2, [0, dep, 0], ["x", "y", "z"]) # PyVista.Plot_BoundaryConditions(simu).show() u, d, K, converg = simu.Solve() simu.Results_Set_Iteration_Summary( i, dep, "mm", i / displacements.size, remove=True ) assert converg simu.Save_Iter() if np.any(d[nodes_xL] >= 1): break if plotIter: Display.Plot_Result(simu, "damage", 1, plotMesh=True, ax=ax) plt.pause(1e-12) simu.Save(folder_save) else: simu: Simulations.PhaseFieldSimu = Simulations.Load_Simu(folder_save) # ---------------------------------------------- # Results # ---------------------------------------------- if makeParaview: 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") Display.Plot_Result(simu, "damage", folder=folder_save) Display.Plot_Result(simu, "uy", deformFactor=1) Display.Plot_Mesh(mesh) Display.Plot_Tags(mesh, folder=folder_save) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 21.263 seconds) .. _sphx_glr_download_examples_PhaseField_CT.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: CT.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: CT.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: CT.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_