.. 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 12-239 .. 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/checkouts/v1.8.5/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100 1: 0.000mm, [0.00e+00;0.00e+00], 1:0.438s, tol=0.00e+00 2: 0.031mm, [0.00e+00;0.00e+00], 1:0.442s, tol=1.00e+00 2.00% -> 21.68s 3: 0.061mm, [0.00e+00;0.00e+00], 1:0.450s, tol=7.50e-01 4.00% -> 21.61s 4: 0.092mm, [0.00e+00;0.00e+00], 1:0.441s, tol=5.56e-01 6.00% -> 20.74s 5: 0.122mm, [0.00e+00;0.00e+00], 1:0.445s, tol=4.37e-01 8.00% -> 20.48s 6: 0.153mm, [0.00e+00;0.00e+00], 1:0.442s, tol=3.60e-01 10.00% -> 19.89s 7: 0.184mm, [0.00e+00;0.00e+00], 1:0.449s, tol=3.06e-01 12.00% -> 19.74s 8: 0.214mm, [0.00e+00;0.00e+00], 1:0.442s, tol=2.65e-01 14.00% -> 19.02s 9: 0.245mm, [1.63e-04;2.46e-02], 1:0.444s, tol=2.32e-01 16.00% -> 18.64s 10: 0.276mm, [2.76e-04;8.92e-02], 1:0.441s, tol=2.04e-01 18.00% -> 18.09s 11: 0.306mm, [3.39e-04;1.99e-01], 1:0.450s, tol=1.81e-01 20.00% -> 18.01s 12: 0.337mm, [3.21e-04;3.42e-01], 1:0.447s, tol=1.65e-01 22.00% -> 17.42s 13: 0.367mm, [2.96e-04;5.27e-01], 1:0.444s, tol=1.52e-01 24.00% -> 16.88s 14: 0.398mm, [2.68e-04;7.44e-01], 1:0.440s, tol=1.59e-01 26.00% -> 16.29s 15: 0.429mm, [2.44e-04;9.33e-01], 1:0.450s, tol=2.32e-01 28.00% -> 16.20s 16: 0.459mm, [2.15e-04;1.02e+00], 1:0.449s, tol=3.72e-01 30.00% -> 15.73s 17: 0.490mm, [1.82e-04;1.02e+00], 1:0.446s, tol=4.13e-01 32.00% -> 15.15s 18: 0.520mm, [1.52e-04;1.01e+00], 1:0.442s, tol=3.98e-01 34.00% -> 14.59s 19: 0.551mm, [1.24e-04;1.01e+00], 1:0.454s, tol=2.99e-01 36.00% -> 14.53s 20: 0.582mm, [1.03e-04;1.01e+00], 1:0.444s, tol=2.73e-01 38.00% -> 13.76s 21: 0.612mm, [8.73e-05;1.01e+00], 1:0.445s, tol=2.61e-01 40.00% -> 13.36s 22: 0.643mm, [7.52e-05;1.01e+00], 1:0.445s, tol=2.47e-01 42.00% -> 12.90s 23: 0.673mm, [6.58e-05;1.01e+00], 1:0.448s, tol=2.42e-01 44.00% -> 12.54s 24: 0.704mm, [5.74e-05;1.00e+00], 1:0.447s, tol=2.07e-01 46.00% -> 12.06s 25: 0.735mm, [5.02e-05;1.00e+00], 1:0.471s, tol=2.11e-01 48.00% -> 12.25s 26: 0.765mm, [4.43e-05;1.00e+00], 1:0.467s, tol=1.77e-01 50.00% -> 11.66s 27: 0.796mm, [3.96e-05;1.01e+00], 1:0.475s, tol=1.67e-01 52.00% -> 11.40s 28: 0.827mm, [3.54e-05;1.01e+00], 1:0.474s, tol=1.74e-01 54.00% -> 10.91s 29: 0.857mm, [3.16e-05;1.00e+00], 1:0.475s, tol=1.49e-01 56.00% -> 10.44s 30: 0.888mm, [2.87e-05;1.00e+00], 1:0.465s, tol=1.39e-01 58.00% -> 9.77s 31: 0.918mm, [2.61e-05;1.00e+00], 1:0.478s, tol=1.29e-01 60.00% -> 9.55s 32: 0.949mm, [2.38e-05;1.00e+00], 1:0.475s, tol=1.31e-01 62.00% -> 9.03s 33: 0.980mm, [2.20e-05;1.00e+00], 1:0.500s, tol=1.26e-01 64.00% -> 9.01s 34: 1.010mm, [2.02e-05;1.00e+00], 1:0.497s, tol=1.12e-01 66.00% -> 8.45s 35: 1.041mm, [1.86e-05;1.00e+00], 1:0.504s, tol=1.11e-01 68.00% -> 8.06s 36: 1.071mm, [1.75e-05;1.00e+00], 1:0.591s, tol=1.09e-01 70.00% -> 8.86s 37: 1.102mm, [1.64e-05;1.01e+00], 1:0.594s, tol=1.08e-01 72.00% -> 8.31s 38: 1.133mm, [1.54e-05;1.01e+00], 1:0.589s, tol=1.09e-01 74.00% -> 7.66s 39: 1.163mm, [1.48e-05;1.01e+00], 1:0.598s, tol=1.13e-01 76.00% -> 7.18s 40: 1.194mm, [1.45e-05;1.01e+00], 1:0.597s, tol=1.24e-01 78.00% -> 6.56s 41: 1.224mm, [1.46e-05;1.01e+00], 1:0.596s, tol=1.55e-01 80.00% -> 5.96s 42: 1.255mm, [1.48e-05;1.00e+00], 1:0.598s, tol=4.20e-01 82.00% -> 5.38s 43: 1.286mm, [1.54e-05;1.00e+00], 1:0.620s, tol=7.43e-01 84.00% -> 4.96s 44: 1.316mm, [1.59e-05;1.00e+00], 1:0.630s, tol=1.25e-01 86.00% -> 4.41s 45: 1.347mm, [1.64e-05;1.01e+00], 1:0.625s, tol=1.33e-01 88.00% -> 3.75s 46: 1.378mm, [1.71e-05;1.00e+00], 1:0.627s, tol=1.92e-02 90.00% -> 3.13s Saved mesh data in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100/Meshes Saved simulation and summary in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100 Generate movie 01/23 (4.35 %) 2.84 s Generate movie 02/23 (8.70 %) 2.21 s Generate movie 03/23 (13.04 %) 2.09 s Generate movie 04/23 (17.39 %) 1.95 s Generate movie 05/23 (21.74 %) 1.83 s Generate movie 06/23 (26.09 %) 1.74 s Generate movie 07/23 (30.43 %) 1.64 s Generate movie 08/23 (34.78 %) 1.56 s Generate movie 09/23 (39.13 %) 1.41 s Generate movie 10/23 (43.48 %) 1.32 s Generate movie 11/23 (47.83 %) 1.22 s Generate movie 12/23 (52.17 %) 1.19 s Generate movie 13/23 (56.52 %) 1.05 s Generate movie 14/23 (60.87 %) 964.45 ms Generate movie 15/23 (65.22 %) 839.84 ms Generate movie 16/23 (69.57 %) 735.61 ms Generate movie 17/23 (73.91 %) 646.79 ms Generate movie 18/23 (78.26 %) 529.15 ms Generate movie 19/23 (82.61 %) 440.93 ms Generate movie 20/23 (86.96 %) 319.21 ms Generate movie 21/23 (91.30 %) 211.63 ms Generate movie 22/23 (95.65 %) 107.92 ms Generate movie 23/23 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 13 import matplotlib.pyplot as plt import numpy as np from EasyFEA import Display, Folder, Models, 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 folder = Folder.Results_Dir() 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 = Models.PhaseField.SplitType.Miehe regu = Models.PhaseField.ReguType.AT1 folder_save = Simulations.PhaseField.Folder( f"{folder}{dim}D", "", split, regu, "", 1, "", meshTest, optimMesh, nL=nL ) Display.MyPrint(folder_save, "green", end="\n") # ---------------------------------------------- # 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.Elastic.Isotropic(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.PhaseField(mesh, pfm, folder=folder_save) 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, converg = simu.Solve() simu.Results_Set_Iteration_Summary( i, dep, "mm", i / displacements.size, remove=True ) assert converg simu.Save_Iter() if simu.Detect_Damage(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.PhaseField = Simulations.Load_Simu(folder_save) # ---------------------------------------------- # Results # ---------------------------------------------- if makeParaview: Paraview.Save_simu(simu, folder_save) if makeMovie: simu.Set_Iter(-1) deformFactor = L * 0.05 / simu.Result("displacement_norm").max() iterations = np.arange(0, simu.Niter, simu.Niter // 20) def Func(plotter, iter): simu.Set_Iter(iterations[iter]) thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8)) PyVista.Plot(thresh, "damage", 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 28.864 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 `_