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.10.1/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
1: 0.000mm, [0.00e+00;0.00e+00], 1:0.117s, tol=0.00e+00
2: 0.031mm, [0.00e+00;0.00e+00], 1:0.089s, tol=1.00e+00 2.00% -> 4.38s
3: 0.061mm, [0.00e+00;0.00e+00], 1:0.091s, tol=7.50e-01 4.00% -> 4.37s
4: 0.092mm, [0.00e+00;0.00e+00], 1:0.092s, tol=5.56e-01 6.00% -> 4.33s
5: 0.122mm, [0.00e+00;0.00e+00], 1:0.091s, tol=4.38e-01 8.00% -> 4.20s
6: 0.153mm, [0.00e+00;0.00e+00], 1:0.091s, tol=3.60e-01 10.00% -> 4.11s
7: 0.184mm, [0.00e+00;0.00e+00], 1:0.091s, tol=3.06e-01 12.00% -> 3.99s
8: 0.214mm, [0.00e+00;0.00e+00], 1:0.091s, tol=2.65e-01 14.00% -> 3.90s
9: 0.245mm, [1.63e-04;2.46e-02], 1:0.091s, tol=2.32e-01 16.00% -> 3.81s
10: 0.276mm, [2.76e-04;8.92e-02], 1:0.093s, tol=2.04e-01 18.00% -> 3.81s
11: 0.306mm, [3.39e-04;1.99e-01], 1:0.092s, tol=1.81e-01 20.00% -> 3.66s
12: 0.337mm, [3.21e-04;3.42e-01], 1:0.093s, tol=1.65e-01 22.00% -> 3.62s
13: 0.367mm, [2.96e-04;5.27e-01], 1:0.092s, tol=1.52e-01 24.00% -> 3.50s
14: 0.398mm, [2.68e-04;7.44e-01], 1:0.091s, tol=1.59e-01 26.00% -> 3.37s
15: 0.429mm, [2.44e-04;9.33e-01], 1:0.091s, tol=2.32e-01 28.00% -> 3.29s
16: 0.459mm, [2.15e-04;1.02e+00], 1:0.093s, tol=3.72e-01 30.00% -> 3.26s
17: 0.490mm, [1.82e-04;1.02e+00], 1:0.091s, tol=4.13e-01 32.00% -> 3.10s
18: 0.520mm, [1.52e-04;1.01e+00], 1:0.091s, tol=3.98e-01 34.00% -> 3.00s
19: 0.551mm, [1.24e-04;1.01e+00], 1:0.092s, tol=2.99e-01 36.00% -> 2.95s
20: 0.582mm, [1.03e-04;1.01e+00], 1:0.092s, tol=2.73e-01 38.00% -> 2.84s
21: 0.612mm, [8.73e-05;1.01e+00], 1:0.093s, tol=2.61e-01 40.00% -> 2.79s
22: 0.643mm, [7.52e-05;1.01e+00], 1:0.093s, tol=2.47e-01 42.00% -> 2.70s
23: 0.673mm, [6.58e-05;1.01e+00], 1:0.092s, tol=2.42e-01 44.00% -> 2.58s
24: 0.704mm, [5.74e-05;1.00e+00], 1:0.092s, tol=2.07e-01 46.00% -> 2.49s
25: 0.735mm, [5.02e-05;1.00e+00], 1:0.092s, tol=2.11e-01 48.00% -> 2.38s
26: 0.765mm, [4.43e-05;1.00e+00], 1:0.093s, tol=1.77e-01 50.00% -> 2.32s
27: 0.796mm, [3.96e-05;1.01e+00], 1:0.091s, tol=1.67e-01 52.00% -> 2.18s
28: 0.827mm, [3.54e-05;1.01e+00], 1:0.091s, tol=1.74e-01 54.00% -> 2.09s
29: 0.857mm, [3.16e-05;1.00e+00], 1:0.091s, tol=1.49e-01 56.00% -> 2.00s
30: 0.888mm, [2.87e-05;1.00e+00], 1:0.092s, tol=1.39e-01 58.00% -> 1.93s
31: 0.918mm, [2.61e-05;1.00e+00], 1:0.094s, tol=1.29e-01 60.00% -> 1.88s
32: 0.949mm, [2.38e-05;1.00e+00], 1:0.094s, tol=1.31e-01 62.00% -> 1.78s
33: 0.980mm, [2.20e-05;1.00e+00], 1:0.094s, tol=1.26e-01 64.00% -> 1.69s
34: 1.010mm, [2.02e-05;1.00e+00], 1:0.094s, tol=1.12e-01 66.00% -> 1.60s
35: 1.041mm, [1.86e-05;1.00e+00], 1:0.093s, tol=1.11e-01 68.00% -> 1.49s
36: 1.071mm, [1.75e-05;1.00e+00], 1:0.093s, tol=1.09e-01 70.00% -> 1.39s
37: 1.102mm, [1.64e-05;1.01e+00], 1:0.093s, tol=1.08e-01 72.00% -> 1.31s
38: 1.133mm, [1.54e-05;1.01e+00], 1:0.093s, tol=1.09e-01 74.00% -> 1.21s
39: 1.163mm, [1.48e-05;1.01e+00], 1:0.092s, tol=1.13e-01 76.00% -> 1.11s
40: 1.194mm, [1.45e-05;1.01e+00], 1:0.092s, tol=1.24e-01 78.00% -> 1.01s
41: 1.224mm, [1.46e-05;1.01e+00], 1:0.093s, tol=1.55e-01 80.00% -> 929.61ms
42: 1.255mm, [1.48e-05;1.00e+00], 1:0.094s, tol=4.20e-01 82.00% -> 844.84ms
43: 1.286mm, [1.54e-05;1.00e+00], 1:0.092s, tol=7.43e-01 84.00% -> 732.66ms
44: 1.316mm, [1.59e-05;1.00e+00], 1:0.092s, tol=1.25e-01 86.00% -> 643.71ms
45: 1.347mm, [1.64e-05;1.01e+00], 1:0.093s, tol=1.33e-01 88.00% -> 558.85ms
46: 1.378mm, [1.71e-05;1.00e+00], 1:0.093s, tol=1.92e-02 90.00% -> 466.22ms
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.10.1/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.10.1/examples/PhaseField/results/CT2D/Test/Miehe_AT1_optimMesh nL=100
Generate movie 01/23 (4.35 %) 3.49 s
Generate movie 02/23 (8.70 %) 3.22 s
Generate movie 03/23 (13.04 %) 3.02 s
Generate movie 04/23 (17.39 %) 2.89 s
Generate movie 05/23 (21.74 %) 2.72 s
Generate movie 06/23 (26.09 %) 2.62 s
Generate movie 07/23 (30.43 %) 2.42 s
Generate movie 08/23 (34.78 %) 2.26 s
Generate movie 09/23 (39.13 %) 2.11 s
Generate movie 10/23 (43.48 %) 1.97 s
Generate movie 11/23 (47.83 %) 1.81 s
Generate movie 12/23 (52.17 %) 1.67 s
Generate movie 13/23 (56.52 %) 1.53 s
Generate movie 14/23 (60.87 %) 1.38 s
Generate movie 15/23 (65.22 %) 1.21 s
Generate movie 16/23 (69.57 %) 1.06 s
Generate movie 17/23 (73.91 %) 922.51 ms
Generate movie 18/23 (78.26 %) 758.30 ms
Generate movie 19/23 (82.61 %) 608.67 ms
Generate movie 20/23 (86.96 %) 449.87 ms
Generate movie 21/23 (91.30 %) 299.93 ms
Generate movie 22/23 (95.65 %) 149.90 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 12.411 seconds)