Note
Go to the end to download the full example code.
CT#
Performs damage simulation on a CT specimen.
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.4/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
1: 0.000mm, [0.00e+00;0.00e+00], 1:0.537s, tol=0.00e+00
2: 0.031mm, [0.00e+00;0.00e+00], 1:0.520s, tol=1.00e+00 2.00% -> 25.48s
3: 0.061mm, [0.00e+00;0.00e+00], 1:0.514s, tol=7.50e-01 4.00% -> 24.67s
4: 0.092mm, [0.00e+00;0.00e+00], 1:0.520s, tol=5.56e-01 6.00% -> 24.43s
5: 0.122mm, [0.00e+00;0.00e+00], 1:0.515s, tol=4.37e-01 8.00% -> 23.69s
6: 0.153mm, [0.00e+00;0.00e+00], 1:0.514s, tol=3.60e-01 10.00% -> 23.15s
7: 0.184mm, [0.00e+00;0.00e+00], 1:0.514s, tol=3.06e-01 12.00% -> 22.61s
8: 0.214mm, [0.00e+00;0.00e+00], 1:0.513s, tol=2.65e-01 14.00% -> 22.05s
9: 0.245mm, [1.63e-04;2.46e-02], 1:0.513s, tol=2.32e-01 16.00% -> 21.56s
10: 0.276mm, [2.76e-04;8.92e-02], 1:0.511s, tol=2.04e-01 18.00% -> 20.95s
11: 0.306mm, [3.39e-04;1.99e-01], 1:0.517s, tol=1.81e-01 20.00% -> 20.68s
12: 0.337mm, [3.21e-04;3.42e-01], 1:0.515s, tol=1.65e-01 22.00% -> 20.08s
13: 0.367mm, [2.96e-04;5.27e-01], 1:0.510s, tol=1.52e-01 24.00% -> 19.37s
14: 0.398mm, [2.68e-04;7.44e-01], 1:0.513s, tol=1.59e-01 26.00% -> 18.98s
15: 0.429mm, [2.44e-04;9.33e-01], 1:0.513s, tol=2.32e-01 28.00% -> 18.45s
16: 0.459mm, [2.15e-04;1.02e+00], 1:0.512s, tol=3.72e-01 30.00% -> 17.91s
17: 0.490mm, [1.82e-04;1.02e+00], 1:0.511s, tol=4.13e-01 32.00% -> 17.39s
18: 0.520mm, [1.52e-04;1.01e+00], 1:0.510s, tol=3.98e-01 34.00% -> 16.82s
19: 0.551mm, [1.24e-04;1.01e+00], 1:0.510s, tol=2.99e-01 36.00% -> 16.34s
20: 0.582mm, [1.03e-04;1.01e+00], 1:0.511s, tol=2.73e-01 38.00% -> 15.84s
21: 0.612mm, [8.73e-05;1.01e+00], 1:0.511s, tol=2.61e-01 40.00% -> 15.33s
22: 0.643mm, [7.52e-05;1.01e+00], 1:0.512s, tol=2.47e-01 42.00% -> 14.84s
23: 0.673mm, [6.58e-05;1.01e+00], 1:0.517s, tol=2.42e-01 44.00% -> 14.47s
24: 0.704mm, [5.74e-05;1.00e+00], 1:0.519s, tol=2.07e-01 46.00% -> 14.01s
25: 0.735mm, [5.02e-05;1.00e+00], 1:0.533s, tol=2.11e-01 48.00% -> 13.86s
26: 0.765mm, [4.43e-05;1.00e+00], 1:0.541s, tol=1.77e-01 50.00% -> 13.52s
27: 0.796mm, [3.96e-05;1.01e+00], 1:0.537s, tol=1.67e-01 52.00% -> 12.89s
28: 0.827mm, [3.54e-05;1.01e+00], 1:0.548s, tol=1.74e-01 54.00% -> 12.61s
29: 0.857mm, [3.16e-05;1.00e+00], 1:0.547s, tol=1.49e-01 56.00% -> 12.03s
30: 0.888mm, [2.87e-05;1.00e+00], 1:0.541s, tol=1.39e-01 58.00% -> 11.36s
31: 0.918mm, [2.61e-05;1.00e+00], 1:0.542s, tol=1.29e-01 60.00% -> 10.83s
32: 0.949mm, [2.38e-05;1.00e+00], 1:0.543s, tol=1.31e-01 62.00% -> 10.32s
33: 0.980mm, [2.20e-05;1.00e+00], 1:0.568s, tol=1.26e-01 64.00% -> 10.23s
34: 1.010mm, [2.02e-05;1.00e+00], 1:0.572s, tol=1.12e-01 66.00% -> 9.72s
35: 1.041mm, [1.86e-05;1.00e+00], 1:0.568s, tol=1.11e-01 68.00% -> 9.09s
36: 1.071mm, [1.75e-05;1.00e+00], 1:0.670s, tol=1.09e-01 70.00% -> 10.06s
37: 1.102mm, [1.64e-05;1.01e+00], 1:0.671s, tol=1.08e-01 72.00% -> 9.39s
38: 1.133mm, [1.54e-05;1.01e+00], 1:0.680s, tol=1.09e-01 74.00% -> 8.83s
39: 1.163mm, [1.48e-05;1.01e+00], 1:0.675s, tol=1.13e-01 76.00% -> 8.10s
40: 1.194mm, [1.45e-05;1.01e+00], 1:0.684s, tol=1.24e-01 78.00% -> 7.52s
41: 1.224mm, [1.46e-05;1.01e+00], 1:0.692s, tol=1.55e-01 80.00% -> 6.92s
42: 1.255mm, [1.48e-05;1.00e+00], 1:0.686s, tol=4.20e-01 82.00% -> 6.17s
43: 1.286mm, [1.54e-05;1.00e+00], 1:0.704s, tol=7.43e-01 84.00% -> 5.63s
44: 1.316mm, [1.59e-05;1.00e+00], 1:0.721s, tol=1.25e-01 86.00% -> 5.05s
45: 1.347mm, [1.64e-05;1.01e+00], 1:0.701s, tol=1.33e-01 88.00% -> 4.20s
46: 1.378mm, [1.71e-05;1.00e+00], 1:0.716s, tol=1.92e-02 90.00% -> 3.58s
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.4/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.4/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
Generate movie 01/23 (4.35 %) 3.06 s
Generate movie 02/23 (8.70 %) 2.64 s
Generate movie 03/23 (13.04 %) 2.41 s
Generate movie 04/23 (17.39 %) 2.30 s
Generate movie 05/23 (21.74 %) 2.16 s
Generate movie 06/23 (26.09 %) 2.05 s
Generate movie 07/23 (30.43 %) 1.96 s
Generate movie 08/23 (34.78 %) 1.85 s
Generate movie 09/23 (39.13 %) 1.69 s
Generate movie 10/23 (43.48 %) 1.57 s
Generate movie 11/23 (47.83 %) 1.45 s
Generate movie 12/23 (52.17 %) 1.35 s
Generate movie 13/23 (56.52 %) 1.24 s
Generate movie 14/23 (60.87 %) 1.09 s
Generate movie 15/23 (65.22 %) 970.83 ms
Generate movie 16/23 (69.57 %) 850.10 ms
Generate movie 17/23 (73.91 %) 744.58 ms
Generate movie 18/23 (78.26 %) 609.21 ms
Generate movie 19/23 (82.61 %) 490.69 ms
Generate movie 20/23 (86.96 %) 363.55 ms
Generate movie 21/23 (91.30 %) 242.63 ms
Generate movie 22/23 (95.65 %) 121.66 ms
Generate movie 23/23 (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, PyVista, Paraview
17 from EasyFEA.Geoms import Point, Points, Circle, Line, Contour, Domain
18
19 if __name__ == "__main__":
20 Display.Clear()
21
22 # ----------------------------------------------
23 # Config
24 # ----------------------------------------------
25
26 dim = 2
27
28 # simu options
29 doSimu = True
30 meshTest = True
31 optimMesh = True
32
33 # outputs
34 folder = Folder.Results_Dir()
35 plotGeom = False
36 plotIter = False
37 makeParaview = False
38 makeMovie = True
39
40 # geom
41 L = 60 # mm
42 e = 4
43 t = 30
44 r = 2
45 t2 = 15
46 diam = 8
47 thickness = 8
48
49 # model
50 Gc = 100
51 nL = 100
52 l0 = L / nL
53 split = Models.PhaseField.SplitType.Miehe
54 regu = Models.PhaseField.ReguType.AT1
55
56 folder_save = Simulations.PhaseField.Folder(
57 f"{folder}{dim}D", "", split, regu, "", 1, "", meshTest, optimMesh, nL=nL
58 )
59
60 Display.MyPrint(folder_save, "green", end="\n")
61
62 # ----------------------------------------------
63 # Geom
64 # ----------------------------------------------
65 clC = l0 if meshTest else l0 / 2 # meshSize on the crack
66 clD = l0 * 2 if optimMesh else clC
67 mS = l0 / 2
68
69 pt1 = Point(0, -L / 2)
70 pt2 = Point(L, -L / 2)
71 pt3 = Point(L, L / 2)
72 pt4 = Point(0, L / 2)
73 pt5 = Point(0, e / 2)
74 pt6 = Point(t + r, e / 2, r=r)
75 pt7 = Point(t + r, -e / 2, r=r)
76 pt8 = Point(0, -e / 2)
77 points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], clD)
78
79 contour = points.Get_Contour()
80
81 circle1 = Circle(Point(t2, -L / 2 + t2), diam, clD)
82 circle2 = Circle(Point(t2, L / 2 - t2), diam, clD)
83
84 if plotGeom:
85 ax = contour.Plot()
86 contour.Plot_Geoms([circle1, circle2], ax=ax)
87
88 # ----------------------------------------------
89 # Mesh
90 # ----------------------------------------------
91
92 refineGeom = (
93 Domain(Point(t, -e * 1.5), Point(L, e * 1.5, thickness), clC)
94 if optimMesh
95 else None
96 )
97
98 if plotGeom:
99 refineGeom.Plot(ax=ax)
100
101 if dim == 2:
102 crack = Line(
103 Point(t + r, 0, isOpen=True), Point(t + r + 6, 0), clC, isOpen=True
104 )
105
106 mesh = contour.Mesh_2D(
107 [circle1, circle2],
108 ElemType.TRI3,
109 cracks=[crack],
110 refineGeoms=[refineGeom],
111 )
112 else:
113 elemType = ElemType.PRISM6
114
115 pc1 = Point(t + r, 0, 0, True)
116 pc2 = Point(t + r + 6, 0, 0)
117 pc3 = pc2 + [0, 0, thickness]
118 pc4 = Point(t + r, 0, thickness, True)
119 #
120 line1 = Line(pc1, pc2, clC, True)
121 line2 = Line(pc2, pc3, clC, False)
122 line3 = Line(pc3, pc4, clC, True)
123 line4 = Line(pc4, pc1, clC, True)
124 #
125 crack = Contour([line1, line2, line3, line4], isOpen=True)
126 cracks = [crack]
127
128 if plotGeom:
129 crack.Plot(ax=ax)
130
131 mesh = contour.Mesh_Extrude(
132 [circle1, circle2],
133 [0, 0, thickness],
134 [4],
135 elemType,
136 cracks=cracks,
137 additionalLines=[line1],
138 refineGeoms=[refineGeom],
139 )
140
141 # PyVista.Plot_Mesh(mesh).show()
142 # PyVista.Plot_Nodes(mesh, mesh.orphanNodes).show()
143
144 nodes_1 = mesh.Nodes_Cylinder(circle1)
145 nodes_2 = mesh.Nodes_Cylinder(circle2)
146
147 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
148
149 # ----------------------------------------------
150 # Simu
151 # ----------------------------------------------
152
153 mat = Models.Elastic.Isotropic(dim, thickness=thickness, planeStress=True)
154
155 pfm = Models.PhaseField(mat, split, regu, Gc, l0)
156
157 if doSimu:
158 displacements = np.linspace(0, L / 40, 50)
159
160 config = """
161 displacements = np.linspace(0, L/40, 50)
162
163 for i, dep in enumerate(displacements):
164
165 if dim == 2:
166 simu.add_dirichlet(nodes_1, [0,-dep], ["x","y"])
167 simu.add_dirichlet(nodes_2, [0,dep], ["x","y"])
168 else:
169 simu.add_dirichlet(nodes_1, [0,-dep, 0], ["x", "y", "z"])
170 simu.add_dirichlet(nodes_2, [0,dep, 0], ["x", "y", "z"])
171 """
172
173 simu = Simulations.PhaseField(mesh, pfm, folder=folder_save)
174 simu.Results_Set_Bc_Summary(config)
175
176 if plotIter:
177 ax = Display.Plot_Result(simu, "damage", 1, plotMesh=True)
178
179 for i, dep in enumerate(displacements):
180 simu.Bc_Init()
181
182 if dim == 2:
183 simu.add_dirichlet(nodes_1, [0, -dep], ["x", "y"])
184 simu.add_dirichlet(nodes_2, [0, dep], ["x", "y"])
185 else:
186 simu.add_dirichlet(nodes_1, [0, -dep, 0], ["x", "y", "z"])
187 simu.add_dirichlet(nodes_2, [0, dep, 0], ["x", "y", "z"])
188
189 # PyVista.Plot_BoundaryConditions(simu).show()
190
191 u, d, converg = simu.Solve()
192
193 simu.Results_Set_Iteration_Summary(
194 i, dep, "mm", i / displacements.size, remove=True
195 )
196
197 assert converg
198
199 simu.Save_Iter()
200
201 if simu.Detect_Damage(nodes_xL, 1):
202 break
203
204 if plotIter:
205 Display.Plot_Result(simu, "damage", 1, plotMesh=True, ax=ax)
206 plt.pause(1e-12)
207
208 simu.Save(folder_save)
209
210 else:
211 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
212
213 # ----------------------------------------------
214 # Results
215 # ----------------------------------------------
216
217 if makeParaview:
218 Paraview.Save_simu(simu, folder_save)
219
220 if makeMovie:
221 simu.Set_Iter(-1)
222 deformFactor = L * 0.05 / simu.Result("displacement_norm").max()
223
224 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
225
226 def Func(plotter, iter):
227 simu.Set_Iter(iterations[iter])
228 thresh = PyVista._pvMesh(simu, "damage", deformFactor).threshold((0, 0.8))
229 PyVista.Plot(thresh, "damage", plotMesh=True, plotter=plotter, clim=(0, 1))
230
231 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
232
233 Display.Plot_Result(simu, "damage", folder=folder_save)
234 Display.Plot_Result(simu, "uy", deformFactor=1)
235 Display.Plot_Mesh(mesh)
236 Display.Plot_Tags(mesh, folder=folder_save)
237
238 plt.show()
Total running time of the script: (0 minutes 32.357 seconds)