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 : 1.847 s
Boundary Conditions : 1.372 ms
Matrix : 3.136 s
Solver : 2.180 s
Resolutions : 4.183 s
PostProcessing : 505.707 ms
Display : 1.490 s
12 import matplotlib.pyplot as plt
13 import numpy as np
14
15 from EasyFEA import (
16 Display,
17 Folder,
18 Models,
19 Tic,
20 Mesher,
21 ElemType,
22 Simulations,
23 Paraview,
24 )
25 from EasyFEA.Geoms import Domain, Point
26
27 if __name__ == "__main__":
28 Display.Clear()
29
30 # ----------------------------------------------
31 # Configuration
32 # ----------------------------------------------
33 dim = 2 # Define the dimension of the problem (2D or 3D)
34
35 # outputs
36 folder = Folder.Results_Dir() + f"{dim}D"
37 plotResult = True
38 makeParaview = False
39
40 # geom
41 L = 120 # mm
42 h = 13 # Height
43 b = 13 # Width
44
45 # model
46 E = 210000 # MPa (Young's modulus)
47 v = 0.25 # Poisson's ratio
48 material = Models.Elastic.Isotropic(dim, thickness=b, E=E, v=v, planeStress=True)
49
50 # load
51 P = 800 # N
52
53 # expected energy
54 WdefRef = 2 * P**2 * L / E / h / b * (L**2 / h / b + (1 + v) * 3 / 5)
55
56 # ----------------------------------------------
57 # Mesh
58 # ----------------------------------------------
59 isOrganised = True
60
61 # List of mesh sizes (number of elements) to investigate convergence
62 if dim == 2:
63 list_N = np.arange(1, 10, 1)
64 else:
65 list_N = np.arange(1, 8, 2)
66
67 # Lists to store data for plotting
68 times_elem_N = [] # times for element type and N size
69 wDef_elem_N = [] # energy
70 dofs_elem_N = [] # dofs
71 zz1_elem_N = [] # zz1
72
73 # ----------------------------------------------
74 # Simulations
75 # ----------------------------------------------
76
77 # Loop over each element type for both 2D and 3D simulations
78 elemTypes = ElemType.Get_2D()[:] if dim == 2 else ElemType.Get_3D()
79
80 # elemTypes = [elem.name for elem in elemTypes.copy()]
81
82 mesher = Mesher()
83
84 for e, elemType in enumerate(elemTypes):
85 times_N = []
86 wDef_N = []
87 dofs_N = []
88 zz1_N = []
89
90 # Loop over each mesh size (number of elements)
91 for N in list_N:
92 meshSize = b / N
93
94 # Define the domain for the mesh
95 domain = Domain(Point(), Point(x=L, y=h), meshSize)
96
97 # Generate the mesh using Gmsh
98 if dim == 2:
99 mesh = mesher.Mesh_2D(domain, [], elemType, isOrganised=isOrganised)
100 volume = mesh.area * material.thickness
101 else:
102 mesh = mesher.Mesh_Extrude(
103 domain,
104 [],
105 elemType=elemType,
106 extrude=[0, 0, b],
107 layers=[4],
108 isOrganised=isOrganised,
109 )
110 volume = mesh.volume
111 # Ensure that the volume matches the expected value (L * h * b)
112 assert np.abs(volume - (L * h * b)) / volume <= 1e-10
113
114 # Define nodes on the left boundary (x=0) and right boundary (x=L)
115 nodes_x0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
116 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
117
118 # Create or update the simulation object with the current mesh
119 if e == 0 and N == list_N[0]:
120 simu = Simulations.Elastic(mesh, material)
121 else:
122 simu.Bc_Init()
123 simu.mesh = mesh
124
125 # Set displacement boundary conditions
126 simu.add_dirichlet(nodes_x0, [0] * dim, simu.Get_unknowns())
127 # Set surface load on the right boundary (y-direction)
128 simu.add_surfLoad(nodes_xL, [-P / h / b], ["y"])
129
130 tic = Tic()
131
132 # Solve the simulation
133 simu.Solve()
134 simu.Save_Iter()
135
136 time = tic.Tac("Resolutions", "Total time", False)
137
138 # Get the computed deformation energy
139 Wdef = simu.Result("Wdef")
140
141 # Store the results for the current mesh size
142 times_N.append(time)
143 wDef_N.append(Wdef)
144 dofs_N.append(mesh.Nn * dim)
145 zz1_N.append(simu.Result("ZZ1"))
146
147 if elemType != mesh.elemType:
148 print("Error in mesh generation")
149
150 print(
151 f"Elem: {mesh.elemType}, nby: {N:2}, Wdef = {np.round(Wdef, 3)}, "
152 f"error = {np.abs(WdefRef - Wdef) / WdefRef:.2e}"
153 )
154
155 # Store the results for the current element type
156 times_elem_N.append(times_N)
157 wDef_elem_N.append(wDef_N)
158 dofs_elem_N.append(dofs_N)
159 zz1_elem_N.append(zz1_N)
160
161 # ----------------------------------------------
162 # Results
163 # ----------------------------------------------
164 # Display the convergence of deformation energy
165 ax_Wdef = Display.Init_Axes()
166 ax_error = Display.Init_Axes()
167 ax_times = Display.Init_Axes()
168 ax_zz1 = Display.Init_Axes()
169
170 print(f"\nWSA = {np.round(WdefRef, 4)} mJ")
171
172 for e, elemType in enumerate(elemTypes):
173 # Convergence of deformation energy
174 ax_Wdef.plot(dofs_elem_N[e], wDef_elem_N[e])
175
176 # Error in deformation energy
177 Wdef = np.array(wDef_elem_N[e])
178 error = (WdefRef - Wdef) / WdefRef * 100
179 ax_error.loglog(dofs_elem_N[e], error)
180
181 # Computation time
182 ax_times.loglog(dofs_elem_N[e], times_elem_N[e])
183 # ax_Times.plot(listDofs_e_nb[e], listTimes_e_nb[e])
184 # ax_Times.set_xscale('log')
185
186 # ZZ1
187 if elemType == elemTypes[0]:
188 last = ax_zz1.loglog(dofs_elem_N[e], error, label=f"{elemType}")
189 ax_zz1.loglog(
190 dofs_elem_N[e],
191 zz1_elem_N[e],
192 ls="--",
193 color=last[0]._color,
194 label=f"{elemType} (ZZ1)",
195 )
196
197 WdefRefArray = np.ones_like(dofs_elem_N[0]) * WdefRef
198 WdefRefArray5 = WdefRefArray * 0.95
199 # WdefRefArray5 = WdefRefArray * 1
200
201 # Deformation energy
202 ax_Wdef.grid()
203 ax_Wdef.set_xlim([-10, np.max(dofs_elem_N[0])])
204 ax_Wdef.set_xlabel("Degrees of Freedom (DOF)")
205 ax_Wdef.set_ylabel("Strain energy W [mJ]")
206 ax_Wdef.legend(elemTypes)
207 # ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color='red')
208 ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color="red")
209 plt.figure(ax_Wdef.figure)
210 Display.Save_fig(folder, "Energy")
211
212 # Error in deformation energy
213 ax_error.grid()
214 ax_error.set_xlabel("Degrees of Freedom (DOF)")
215 ax_error.set_ylabel("Error W [%]")
216 ax_error.legend(elemTypes)
217 plt.figure(ax_error.figure)
218 Display.Save_fig(folder, "Error")
219
220 # Error in deformation energy
221 ax_zz1.grid()
222 ax_zz1.set_xlabel("Degrees of Freedom (DOF)")
223 ax_zz1.set_ylabel("Error [%]")
224 ax_zz1.legend()
225 plt.figure(ax_zz1.figure)
226 Display.Save_fig(folder, "Error ZZ1")
227
228 # Computation time
229 ax_times.grid()
230 ax_times.set_xlabel("Degrees of Freedom (DOF)")
231 ax_times.set_ylabel("Computation Time [s]")
232 ax_times.legend(elemTypes)
233 plt.figure(ax_times.figure)
234 Display.Save_fig(folder, "Time")
235
236 # Plot the von Mises stress result using 20 color levels
237 Display.Plot_Result(simu, "Svm", ncolors=20)
238
239 if makeParaview:
240 # Generate Paraview files for visualization
241 Paraview.Save_simu(simu, folder, details=True)
242
243 # Show the total computation time
244 print()
245 Tic.Resume()
246
247 # Display the computation time history
248 # Tic.Plot_History(folder)
249
250 # Show all plots
251 plt.show()
Total running time of the script: (0 minutes 9.205 seconds)