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.6.5/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100 1 : 0.000 mm, [0.00e+00; 0.00e+00], 1:0.443 s, tol=0.00e+00
2 : 0.031 mm, [0.00e+00; 0.00e+00], 1:0.417 s, tol=1.00e+00 2.00 % -> 20.43 s
3 : 0.061 mm, [0.00e+00; 0.00e+00], 1:0.424 s, tol=1.00e+00 4.00 % -> 20.34 s
4 : 0.092 mm, [0.00e+00; 0.00e+00], 1:0.419 s, tol=5.56e-01 6.00 % -> 19.67 s
5 : 0.122 mm, [0.00e+00; 0.00e+00], 1:0.420 s, tol=4.37e-01 8.00 % -> 19.33 s
6 : 0.153 mm, [0.00e+00; 0.00e+00], 1:0.418 s, tol=3.60e-01 10.00 % -> 18.83 s
7 : 0.184 mm, [0.00e+00; 0.00e+00], 1:0.419 s, tol=3.06e-01 12.00 % -> 18.43 s
8 : 0.214 mm, [0.00e+00; 0.00e+00], 1:0.417 s, tol=2.65e-01 14.00 % -> 17.91 s
9 : 0.245 mm, [1.53e-04; 2.46e-02], 1:0.420 s, tol=2.32e-01 16.00 % -> 17.63 s
10 : 0.276 mm, [2.60e-04; 8.92e-02], 1:0.417 s, tol=2.04e-01 18.00 % -> 17.12 s
11 : 0.306 mm, [3.24e-04; 1.99e-01], 1:0.417 s, tol=1.81e-01 20.00 % -> 16.67 s
12 : 0.337 mm, [3.13e-04; 3.42e-01], 1:0.419 s, tol=1.65e-01 22.00 % -> 16.35 s
13 : 0.367 mm, [2.86e-04; 5.27e-01], 1:0.419 s, tol=1.52e-01 24.00 % -> 15.92 s
14 : 0.398 mm, [2.60e-04; 7.44e-01], 1:0.417 s, tol=1.59e-01 26.00 % -> 15.45 s
15 : 0.429 mm, [2.37e-04; 9.33e-01], 1:0.418 s, tol=2.32e-01 28.00 % -> 15.05 s
16 : 0.459 mm, [2.09e-04; 1.02e+00], 1:0.416 s, tol=3.72e-01 30.00 % -> 14.55 s
17 : 0.490 mm, [1.78e-04; 1.02e+00], 1:0.420 s, tol=4.13e-01 32.00 % -> 14.28 s
18 : 0.520 mm, [1.49e-04; 1.01e+00], 1:0.418 s, tol=3.98e-01 34.00 % -> 13.80 s
19 : 0.551 mm, [1.22e-04; 1.01e+00], 1:0.418 s, tol=2.98e-01 36.00 % -> 13.37 s
20 : 0.582 mm, [1.01e-04; 1.01e+00], 1:0.417 s, tol=2.73e-01 38.00 % -> 12.93 s
21 : 0.612 mm, [8.60e-05; 1.01e+00], 1:0.417 s, tol=2.61e-01 40.00 % -> 12.51 s
22 : 0.643 mm, [7.42e-05; 1.01e+00], 1:0.420 s, tol=2.47e-01 42.00 % -> 12.19 s
23 : 0.673 mm, [6.50e-05; 1.01e+00], 1:0.417 s, tol=2.42e-01 44.00 % -> 11.67 s
24 : 0.704 mm, [5.68e-05; 1.00e+00], 1:0.421 s, tol=2.07e-01 46.00 % -> 11.36 s
25 : 0.735 mm, [4.99e-05; 1.00e+00], 1:0.442 s, tol=2.11e-01 48.00 % -> 11.50 s
26 : 0.765 mm, [4.41e-05; 1.00e+00], 1:0.440 s, tol=1.77e-01 50.00 % -> 11.01 s
27 : 0.796 mm, [3.95e-05; 1.01e+00], 1:0.440 s, tol=1.67e-01 52.00 % -> 10.55 s
28 : 0.827 mm, [3.54e-05; 1.01e+00], 1:0.447 s, tol=1.74e-01 54.00 % -> 10.29 s
29 : 0.857 mm, [3.16e-05; 1.00e+00], 1:0.450 s, tol=1.49e-01 56.00 % -> 9.90 s
30 : 0.888 mm, [2.86e-05; 1.00e+00], 1:0.442 s, tol=1.39e-01 58.00 % -> 9.29 s
31 : 0.918 mm, [2.61e-05; 1.00e+00], 1:0.442 s, tol=1.29e-01 60.00 % -> 8.85 s
32 : 0.949 mm, [2.38e-05; 1.00e+00], 1:0.443 s, tol=1.31e-01 62.00 % -> 8.42 s
33 : 0.980 mm, [2.20e-05; 1.00e+00], 1:0.471 s, tol=1.26e-01 64.00 % -> 8.48 s
34 : 1.010 mm, [2.01e-05; 1.00e+00], 1:0.471 s, tol=1.12e-01 66.00 % -> 8.01 s
35 : 1.041 mm, [1.86e-05; 1.00e+00], 1:0.471 s, tol=1.11e-01 68.00 % -> 7.54 s
36 : 1.071 mm, [1.74e-05; 1.00e+00], 1:0.558 s, tol=1.09e-01 70.00 % -> 8.36 s
37 : 1.102 mm, [1.64e-05; 1.01e+00], 1:0.561 s, tol=1.08e-01 72.00 % -> 7.86 s
38 : 1.133 mm, [1.54e-05; 1.01e+00], 1:0.558 s, tol=1.09e-01 74.00 % -> 7.26 s
39 : 1.163 mm, [1.48e-05; 1.01e+00], 1:0.563 s, tol=1.13e-01 76.00 % -> 6.75 s
40 : 1.194 mm, [1.45e-05; 1.01e+00], 1:0.564 s, tol=1.24e-01 78.00 % -> 6.21 s
41 : 1.224 mm, [1.46e-05; 1.01e+00], 1:0.563 s, tol=1.55e-01 80.00 % -> 5.63 s
42 : 1.255 mm, [1.47e-05; 1.00e+00], 1:0.565 s, tol=4.20e-01 82.00 % -> 5.08 s
43 : 1.286 mm, [1.54e-05; 1.00e+00], 1:0.581 s, tol=7.43e-01 84.00 % -> 4.64 s
44 : 1.316 mm, [1.59e-05; 1.00e+00], 1:0.597 s, tol=1.25e-01 86.00 % -> 4.18 s
45 : 1.347 mm, [1.64e-05; 1.01e+00], 1:0.590 s, tol=1.33e-01 88.00 % -> 3.54 s
46 : 1.378 mm, [1.71e-05; 1.00e+00], 1:0.597 s, tol=1.92e-02 90.00 % -> 2.98 s
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100/simulation.pickle
Saved:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.6.5/examples/PhaseField/results/CT2D/Test/Isot_Miehe_AT1_optimMesh nL=100/summary.txt
Generate movie 01/23 (4.35 %) 2.53 s
Generate movie 02/23 (8.70 %) 2.18 s
Generate movie 03/23 (13.04 %) 2.01 s
Generate movie 04/23 (17.39 %) 1.89 s
Generate movie 05/23 (21.74 %) 1.80 s
Generate movie 06/23 (26.09 %) 1.77 s
Generate movie 07/23 (30.43 %) 1.63 s
Generate movie 08/23 (34.78 %) 1.53 s
Generate movie 09/23 (39.13 %) 1.39 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.12 s
Generate movie 13/23 (56.52 %) 992.03 ms
Generate movie 14/23 (60.87 %) 898.98 ms
Generate movie 15/23 (65.22 %) 792.22 ms
Generate movie 16/23 (69.57 %) 708.47 ms
Generate movie 17/23 (73.91 %) 616.77 ms
Generate movie 18/23 (78.26 %) 498.99 ms
Generate movie 19/23 (82.61 %) 400.58 ms
Generate movie 20/23 (86.96 %) 296.64 ms
Generate movie 21/23 (91.30 %) 196.33 ms
Generate movie 22/23 (95.65 %) 101.16 ms
Generate movie 23/23 (100.00 %) 0.00 µs
12 from EasyFEA import (
13 Display,
14 Folder,
15 Models,
16 plt,
17 np,
18 ElemType,
19 Simulations,
20 PyVista,
21 Paraview,
22 )
23 from EasyFEA.Geoms import Point, Points, Circle, Line, Contour, Domain
24
25 if __name__ == "__main__":
26 Display.Clear()
27
28 # ----------------------------------------------
29 # Config
30 # ----------------------------------------------
31
32 dim = 2
33
34 # simu options
35 doSimu = True
36 meshTest = True
37 optimMesh = True
38
39 # outputs
40 folder = Folder.Results_Dir()
41 plotGeom = False
42 plotIter = False
43 makeParaview = False
44 makeMovie = True
45
46 # geom
47 L = 60 # mm
48 e = 4
49 t = 30
50 r = 2
51 t2 = 15
52 diam = 8
53 thickness = 8
54
55 # model
56 Gc = 100
57 nL = 100
58 l0 = L / nL
59 split = "Miehe"
60 regu = "AT1"
61
62 folder_save = Simulations.PhaseField.Folder(
63 f"{folder}{dim}D", "Isot", split, regu, "", 1, "", meshTest, optimMesh, nL=nL
64 )
65
66 Display.MyPrint(folder_save, "green")
67
68 # ----------------------------------------------
69 # Geom
70 # ----------------------------------------------
71 clC = l0 if meshTest else l0 / 2 # meshSize on the crack
72 clD = l0 * 2 if optimMesh else clC
73 mS = l0 / 2
74
75 pt1 = Point(0, -L / 2)
76 pt2 = Point(L, -L / 2)
77 pt3 = Point(L, L / 2)
78 pt4 = Point(0, L / 2)
79 pt5 = Point(0, e / 2)
80 pt6 = Point(t + r, e / 2, r=r)
81 pt7 = Point(t + r, -e / 2, r=r)
82 pt8 = Point(0, -e / 2)
83 points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], clD)
84
85 contour = points.Get_Contour()
86
87 circle1 = Circle(Point(t2, -L / 2 + t2), diam, clD)
88 circle2 = Circle(Point(t2, L / 2 - t2), diam, clD)
89
90 if plotGeom:
91 ax = contour.Plot()
92 contour.Plot_Geoms([circle1, circle2], ax=ax)
93
94 # ----------------------------------------------
95 # Mesh
96 # ----------------------------------------------
97
98 refineGeom = (
99 Domain(Point(t, -e * 1.5), Point(L, e * 1.5, thickness), clC)
100 if optimMesh
101 else None
102 )
103
104 if plotGeom:
105 refineGeom.Plot(ax=ax)
106
107 if dim == 2:
108 crack = Line(
109 Point(t + r, 0, isOpen=True), Point(t + r + 6, 0), clC, isOpen=True
110 )
111
112 mesh = contour.Mesh_2D(
113 [circle1, circle2],
114 ElemType.TRI3,
115 cracks=[crack],
116 refineGeoms=[refineGeom],
117 )
118 else:
119 elemType = ElemType.PRISM6
120
121 pc1 = Point(t + r, 0, 0, True)
122 pc2 = Point(t + r + 6, 0, 0)
123 pc3 = pc2 + [0, 0, thickness]
124 pc4 = Point(t + r, 0, thickness, True)
125 #
126 line1 = Line(pc1, pc2, clC, True)
127 line2 = Line(pc2, pc3, clC, False)
128 line3 = Line(pc3, pc4, clC, True)
129 line4 = Line(pc4, pc1, clC, True)
130 #
131 crack = Contour([line1, line2, line3, line4], isOpen=True)
132 cracks = [crack]
133
134 if plotGeom:
135 crack.Plot(ax=ax)
136
137 mesh = contour.Mesh_Extrude(
138 [circle1, circle2],
139 [0, 0, thickness],
140 [4],
141 elemType,
142 cracks=cracks,
143 additionalLines=[line1],
144 refineGeoms=[refineGeom],
145 )
146
147 # PyVista.Plot_Mesh(mesh).show()
148 # PyVista.Plot_Nodes(mesh, mesh.orphanNodes).show()
149
150 nodes_1 = mesh.Nodes_Cylinder(circle1)
151 nodes_2 = mesh.Nodes_Cylinder(circle2)
152
153 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
154
155 # ----------------------------------------------
156 # Simu
157 # ----------------------------------------------
158
159 mat = Models.Elastic.Isotropic(dim, thickness=thickness, planeStress=True)
160
161 pfm = Models.PhaseField(mat, split, regu, Gc, l0)
162
163 if doSimu:
164 displacements = np.linspace(0, L / 40, 50)
165
166 config = """
167 displacements = np.linspace(0, L/40, 50)
168
169 for i, dep in enumerate(displacements):
170
171 if dim == 2:
172 simu.add_dirichlet(nodes_1, [0,-dep], ["x","y"])
173 simu.add_dirichlet(nodes_2, [0,dep], ["x","y"])
174 else:
175 simu.add_dirichlet(nodes_1, [0,-dep, 0], ["x", "y", "z"])
176 simu.add_dirichlet(nodes_2, [0,dep, 0], ["x", "y", "z"])
177 """
178
179 simu = Simulations.PhaseField(mesh, pfm)
180 simu.Results_Set_Bc_Summary(config)
181
182 if plotIter:
183 ax = Display.Plot_Result(simu, "damage", 1, plotMesh=True)
184
185 for i, dep in enumerate(displacements):
186 simu.Bc_Init()
187
188 if dim == 2:
189 simu.add_dirichlet(nodes_1, [0, -dep], ["x", "y"])
190 simu.add_dirichlet(nodes_2, [0, dep], ["x", "y"])
191 else:
192 simu.add_dirichlet(nodes_1, [0, -dep, 0], ["x", "y", "z"])
193 simu.add_dirichlet(nodes_2, [0, dep, 0], ["x", "y", "z"])
194
195 # PyVista.Plot_BoundaryConditions(simu).show()
196
197 u, d, K, converg = simu.Solve()
198
199 simu.Results_Set_Iteration_Summary(
200 i, dep, "mm", i / displacements.size, remove=True
201 )
202
203 assert converg
204
205 simu.Save_Iter()
206
207 if np.any(d[nodes_xL] >= 1):
208 break
209
210 if plotIter:
211 Display.Plot_Result(simu, "damage", 1, plotMesh=True, ax=ax)
212 plt.pause(1e-12)
213
214 simu.Save(folder_save)
215
216 else:
217 simu: Simulations.PhaseField = Simulations.Load_Simu(folder_save)
218
219 # ----------------------------------------------
220 # Results
221 # ----------------------------------------------
222
223 if makeParaview:
224 Paraview.Save_simu(simu, folder_save)
225
226 if makeMovie:
227 simu.Set_Iter(-1)
228 depMax = simu.Result("displacement_norm").max()
229 deformFactor = L * 0.05 / depMax
230
231 iterations = np.arange(0, simu.Niter, simu.Niter // 20)
232
233 def Func(plotter, iter):
234 simu.Set_Iter(iterations[iter])
235
236 grid = PyVista._pvMesh(simu, "damage", deformFactor)
237
238 tresh = grid.threshold((0, 0.8))
239
240 PyVista.Plot(
241 tresh,
242 "damage",
243 deformFactor,
244 plotMesh=True,
245 plotter=plotter,
246 clim=(0, 1),
247 )
248
249 PyVista.Movie_func(Func, iterations.size, folder_save, "damage.gif")
250
251 Display.Plot_Result(simu, "damage", folder=folder_save)
252 Display.Plot_Result(simu, "uy", deformFactor=1)
253 Display.Plot_Mesh(mesh)
254 Display.Plot_Tags(mesh, folder=folder_save)
255
256 plt.show()
Total running time of the script: (0 minutes 27.447 seconds)