Note
Go to the end to download the full example code.
MeshConvergence#
Verification of energy convergence for a bending beam for all available elements.
Elem: TRI3, nby: 1, Wdef = 85.307, error = 7.71e-01
Elem: TRI3, nby: 2, Wdef = 199.836, error = 4.63e-01
Elem: TRI3, nby: 3, Wdef = 268.382, error = 2.79e-01
Elem: TRI3, nby: 4, Wdef = 305.294, error = 1.79e-01
Elem: TRI3, nby: 5, Wdef = 326.756, error = 1.22e-01
Elem: TRI3, nby: 6, Wdef = 339.109, error = 8.85e-02
Elem: TRI3, nby: 7, Wdef = 347.06, error = 6.71e-02
Elem: TRI3, nby: 8, Wdef = 352.445, error = 5.26e-02
Elem: TRI3, nby: 9, Wdef = 356.373, error = 4.21e-02
Elem: TRI6, nby: 1, Wdef = 369.006, error = 8.10e-03
Elem: TRI6, nby: 2, Wdef = 370.965, error = 2.84e-03
Elem: TRI6, nby: 3, Wdef = 371.253, error = 2.06e-03
Elem: TRI6, nby: 4, Wdef = 371.354, error = 1.79e-03
Elem: TRI6, nby: 5, Wdef = 371.405, error = 1.65e-03
Elem: TRI6, nby: 6, Wdef = 371.433, error = 1.58e-03
Elem: TRI6, nby: 7, Wdef = 371.451, error = 1.53e-03
Elem: TRI6, nby: 8, Wdef = 371.463, error = 1.50e-03
Elem: TRI6, nby: 9, Wdef = 371.472, error = 1.47e-03
Elem: TRI10, nby: 1, Wdef = 371.008, error = 2.72e-03
Elem: TRI10, nby: 2, Wdef = 371.36, error = 1.78e-03
Elem: TRI10, nby: 3, Wdef = 371.433, error = 1.58e-03
Elem: TRI10, nby: 4, Wdef = 371.463, error = 1.50e-03
Elem: TRI10, nby: 5, Wdef = 371.479, error = 1.45e-03
Elem: TRI10, nby: 6, Wdef = 371.489, error = 1.43e-03
Elem: TRI10, nby: 7, Wdef = 371.495, error = 1.41e-03
Elem: TRI10, nby: 8, Wdef = 371.499, error = 1.40e-03
Elem: TRI10, nby: 9, Wdef = 371.503, error = 1.39e-03
Elem: TRI15, nby: 1, Wdef = 371.309, error = 1.91e-03
Elem: TRI15, nby: 2, Wdef = 371.445, error = 1.55e-03
Elem: TRI15, nby: 3, Wdef = 371.478, error = 1.46e-03
Elem: TRI15, nby: 4, Wdef = 371.493, error = 1.42e-03
Elem: TRI15, nby: 5, Wdef = 371.501, error = 1.40e-03
Elem: TRI15, nby: 6, Wdef = 371.505, error = 1.39e-03
Elem: TRI15, nby: 7, Wdef = 371.508, error = 1.38e-03
Elem: TRI15, nby: 8, Wdef = 371.511, error = 1.37e-03
Elem: TRI15, nby: 9, Wdef = 371.512, error = 1.37e-03
Elem: QUAD4, nby: 1, Wdef = 249.865, error = 3.28e-01
Elem: QUAD4, nby: 2, Wdef = 330.381, error = 1.12e-01
Elem: QUAD4, nby: 3, Wdef = 351.849, error = 5.42e-02
Elem: QUAD4, nby: 4, Wdef = 360.11, error = 3.20e-02
Elem: QUAD4, nby: 5, Wdef = 364.354, error = 2.06e-02
Elem: QUAD4, nby: 6, Wdef = 366.455, error = 1.50e-02
Elem: QUAD4, nby: 7, Wdef = 367.75, error = 1.15e-02
Elem: QUAD4, nby: 8, Wdef = 368.603, error = 9.19e-03
Elem: QUAD4, nby: 9, Wdef = 369.241, error = 7.47e-03
Elem: QUAD8, nby: 1, Wdef = 369.262, error = 7.42e-03
Elem: QUAD8, nby: 2, Wdef = 371.115, error = 2.43e-03
Elem: QUAD8, nby: 3, Wdef = 371.315, error = 1.90e-03
Elem: QUAD8, nby: 4, Wdef = 371.386, error = 1.70e-03
Elem: QUAD8, nby: 5, Wdef = 371.424, error = 1.60e-03
Elem: QUAD8, nby: 6, Wdef = 371.445, error = 1.55e-03
Elem: QUAD8, nby: 7, Wdef = 371.46, error = 1.51e-03
Elem: QUAD8, nby: 8, Wdef = 371.47, error = 1.48e-03
Elem: QUAD8, nby: 9, Wdef = 371.478, error = 1.46e-03
Elem: QUAD9, nby: 1, Wdef = 370.315, error = 4.58e-03
Elem: QUAD9, nby: 2, Wdef = 371.231, error = 2.12e-03
Elem: QUAD9, nby: 3, Wdef = 371.375, error = 1.74e-03
Elem: QUAD9, nby: 4, Wdef = 371.428, error = 1.59e-03
Elem: QUAD9, nby: 5, Wdef = 371.456, error = 1.52e-03
Elem: QUAD9, nby: 6, Wdef = 371.471, error = 1.48e-03
Elem: QUAD9, nby: 7, Wdef = 371.481, error = 1.45e-03
Elem: QUAD9, nby: 8, Wdef = 371.488, error = 1.43e-03
Elem: QUAD9, nby: 9, Wdef = 371.494, error = 1.42e-03
WSA = 372.0206 mJ
Mesh : 5.292 s
Boundary Conditions : 2.213 ms
Matrix : 3.569 s
Solver : 3.817 s
Display : 2.140 s
PostProcessing : 541.164 ms
PyVista_Interface : 29.663 s
Paraview : 818.761 ms
Resolution hyperelastic : 922.252 ms
PyVista : 3.815 s
Resolutions : 1.436 s
12 from EasyFEA import (
13 Display,
14 Folder,
15 Models,
16 Tic,
17 plt,
18 np,
19 Mesher,
20 ElemType,
21 Simulations,
22 Paraview,
23 )
24 from EasyFEA.Geoms import Domain, Point
25
26 if __name__ == "__main__":
27 Display.Clear()
28
29 # ----------------------------------------------
30 # Configuration
31 # ----------------------------------------------
32 dim = 2 # Define the dimension of the problem (2D or 3D)
33
34 # outputs
35 folder = Folder.Join(Folder.RESULTS_DIR, "Meshes", f"Convergence{dim}D")
36 plotResult = True
37 makeParaview = False
38
39 # geom
40 L = 120 # mm
41 h = 13 # Height
42 b = 13 # Width
43
44 # model
45 E = 210000 # MPa (Young's modulus)
46 v = 0.25 # Poisson's ratio
47 material = Models.ElasIsot(dim, thickness=b, E=E, v=v, planeStress=True)
48
49 # load
50 P = 800 # N
51
52 # expected energy
53 WdefRef = 2 * P**2 * L / E / h / b * (L**2 / h / b + (1 + v) * 3 / 5)
54
55 # ----------------------------------------------
56 # Mesh
57 # ----------------------------------------------
58 isOrganised = True
59
60 # List of mesh sizes (number of elements) to investigate convergence
61 if dim == 2:
62 list_N = np.arange(1, 10, 1)
63 else:
64 list_N = np.arange(1, 8, 2)
65
66 # Lists to store data for plotting
67 times_elem_N = [] # times for element type and N size
68 wDef_elem_N = [] # energy
69 dofs_elem_N = [] # dofs
70 zz1_elem_N = [] # zz1
71
72 # ----------------------------------------------
73 # Simulations
74 # ----------------------------------------------
75
76 # Loop over each element type for both 2D and 3D simulations
77 elemTypes = ElemType.Get_2D()[:] if dim == 2 else ElemType.Get_3D()
78
79 # elemTypes = [elem.name for elem in elemTypes.copy()]
80
81 mesher = Mesher()
82
83 for e, elemType in enumerate(elemTypes):
84 times_N = []
85 wDef_N = []
86 dofs_N = []
87 zz1_N = []
88
89 # Loop over each mesh size (number of elements)
90 for N in list_N:
91 meshSize = b / N
92
93 # Define the domain for the mesh
94 domain = Domain(Point(), Point(x=L, y=h), meshSize)
95
96 # Generate the mesh using Gmsh
97 if dim == 2:
98 mesh = mesher.Mesh_2D(domain, [], elemType, isOrganised=isOrganised)
99 volume = mesh.area * material.thickness
100 else:
101 mesh = mesher.Mesh_Extrude(
102 domain,
103 [],
104 elemType=elemType,
105 extrude=[0, 0, b],
106 layers=[4],
107 isOrganised=isOrganised,
108 )
109 volume = mesh.volume
110 # Ensure that the volume matches the expected value (L * h * b)
111 assert np.abs(volume - (L * h * b)) / volume <= 1e-10
112
113 # Define nodes on the left boundary (x=0) and right boundary (x=L)
114 nodes_x0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
115 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
116
117 # Create or update the simulation object with the current mesh
118 if e == 0 and N == list_N[0]:
119 simu = Simulations.ElasticSimu(
120 mesh, material, useIterativeSolvers=False
121 )
122 else:
123 simu.Bc_Init()
124 simu.mesh = mesh
125
126 # Set displacement boundary conditions
127 simu.add_dirichlet(nodes_x0, [0] * dim, simu.Get_unknowns())
128 # Set surface load on the right boundary (y-direction)
129 simu.add_surfLoad(nodes_xL, [-P / h / b], ["y"])
130
131 tic = Tic()
132
133 # Solve the simulation
134 simu.Solve()
135 simu.Save_Iter()
136
137 time = tic.Tac("Resolutions", "Total time", False)
138
139 # Get the computed deformation energy
140 Wdef = simu.Result("Wdef")
141
142 # Store the results for the current mesh size
143 times_N.append(time)
144 wDef_N.append(Wdef)
145 dofs_N.append(mesh.Nn * dim)
146 zz1_N.append(simu.Result("ZZ1"))
147
148 if elemType != mesh.elemType:
149 print("Error in mesh generation")
150
151 print(
152 f"Elem: {mesh.elemType}, nby: {N:2}, Wdef = {np.round(Wdef, 3)}, "
153 f"error = {np.abs(WdefRef - Wdef) / WdefRef:.2e}"
154 )
155
156 # Store the results for the current element type
157 times_elem_N.append(times_N)
158 wDef_elem_N.append(wDef_N)
159 dofs_elem_N.append(dofs_N)
160 zz1_elem_N.append(zz1_N)
161
162 # ----------------------------------------------
163 # Results
164 # ----------------------------------------------
165 # Display the convergence of deformation energy
166 ax_Wdef = Display.Init_Axes()
167 ax_error = Display.Init_Axes()
168 ax_times = Display.Init_Axes()
169 ax_zz1 = Display.Init_Axes()
170
171 print(f"\nWSA = {np.round(WdefRef, 4)} mJ")
172
173 for e, elemType in enumerate(elemTypes):
174 # Convergence of deformation energy
175 ax_Wdef.plot(dofs_elem_N[e], wDef_elem_N[e])
176
177 # Error in deformation energy
178 Wdef = np.array(wDef_elem_N[e])
179 error = (WdefRef - Wdef) / WdefRef * 100
180 ax_error.loglog(dofs_elem_N[e], error)
181
182 # Computation time
183 ax_times.loglog(dofs_elem_N[e], times_elem_N[e])
184 # ax_Times.plot(listDofs_e_nb[e], listTimes_e_nb[e])
185 # ax_Times.set_xscale('log')
186
187 # ZZ1
188 if elemType == elemTypes[0]:
189 last = ax_zz1.loglog(dofs_elem_N[e], error, label=f"{elemType}")
190 ax_zz1.loglog(
191 dofs_elem_N[e],
192 zz1_elem_N[e],
193 ls="--",
194 color=last[0]._color,
195 label=f"{elemType} (ZZ1)",
196 )
197
198 WdefRefArray = np.ones_like(dofs_elem_N[0]) * WdefRef
199 WdefRefArray5 = WdefRefArray * 0.95
200 # WdefRefArray5 = WdefRefArray * 1
201
202 # Deformation energy
203 ax_Wdef.grid()
204 ax_Wdef.set_xlim([-10, np.max(dofs_elem_N[0])])
205 ax_Wdef.set_xlabel("Degrees of Freedom (DOF)")
206 ax_Wdef.set_ylabel("Strain energy W [mJ]")
207 ax_Wdef.legend(elemTypes)
208 # ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color='red')
209 ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color="red")
210 plt.figure(ax_Wdef.figure)
211 Display.Save_fig(folder, "Energy")
212
213 # Error in deformation energy
214 ax_error.grid()
215 ax_error.set_xlabel("Degrees of Freedom (DOF)")
216 ax_error.set_ylabel("Error W [%]")
217 ax_error.legend(elemTypes)
218 plt.figure(ax_error.figure)
219 Display.Save_fig(folder, "Error")
220
221 # Error in deformation energy
222 ax_zz1.grid()
223 ax_zz1.set_xlabel("Degrees of Freedom (DOF)")
224 ax_zz1.set_ylabel("Error [%]")
225 ax_zz1.legend()
226 plt.figure(ax_zz1.figure)
227 Display.Save_fig(folder, "Error ZZ1")
228
229 # Computation time
230 ax_times.grid()
231 ax_times.set_xlabel("Degrees of Freedom (DOF)")
232 ax_times.set_ylabel("Computation Time [s]")
233 ax_times.legend(elemTypes)
234 plt.figure(ax_times.figure)
235 Display.Save_fig(folder, "Time")
236
237 # Plot the von Mises stress result using 20 color levels
238 Display.Plot_Result(simu, "Svm", ncolors=20)
239
240 if makeParaview:
241 # Generate Paraview files for visualization
242 Paraview.Save_simu(simu, folder, details=True)
243
244 # Show the total computation time
245 print()
246 Tic.Resume()
247
248 # Display the computation time history
249 # Tic.Plot_History(folder)
250
251 # Show all plots
252 plt.show()
Total running time of the script: (0 minutes 3.612 seconds)