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 : 703.348 ms
Boundary Conditions : 707.388 µs
Matrix : 1.008 s
Solver : 696.789 ms
Resolutions : 1.379 s
PostProcessing : 170.383 ms
Display : 825.457 ms
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import (
17 Display,
18 Folder,
19 Models,
20 Tic,
21 Mesher,
22 ElemType,
23 Simulations,
24 Paraview,
25 )
26 from EasyFEA.Geoms import Domain, Point
27
28 if __name__ == "__main__":
29 Display.Clear()
30
31 # ----------------------------------------------
32 # Configuration
33 # ----------------------------------------------
34 dim = 2 # Define the dimension of the problem (2D or 3D)
35
36 # outputs
37 folder = Folder.Results_Dir() + f"{dim}D"
38 plotResult = True
39 makeParaview = False
40
41 # geom
42 L = 120 # mm
43 h = 13 # Height
44 b = 13 # Width
45
46 # model
47 E = 210000 # MPa (Young's modulus)
48 v = 0.25 # Poisson's ratio
49 material = Models.Elastic.Isotropic(dim, thickness=b, E=E, v=v, planeStress=True)
50
51 # load
52 P = 800 # N
53
54 # expected energy
55 WdefRef = 2 * P**2 * L / E / h / b * (L**2 / h / b + (1 + v) * 3 / 5)
56
57 # ----------------------------------------------
58 # Mesh
59 # ----------------------------------------------
60 isOrganised = True
61
62 # List of mesh sizes (number of elements) to investigate convergence
63 if dim == 2:
64 list_N = np.arange(1, 10, 1)
65 else:
66 list_N = np.arange(1, 8, 2)
67
68 # Lists to store data for plotting
69 times_elem_N = [] # times for element type and N size
70 wDef_elem_N = [] # energy
71 dofs_elem_N = [] # dofs
72 zz1_elem_N = [] # zz1
73
74 # ----------------------------------------------
75 # Simulations
76 # ----------------------------------------------
77
78 # Loop over each element type for both 2D and 3D simulations
79 elemTypes = ElemType.Get_2D()[:] if dim == 2 else ElemType.Get_3D()
80
81 # elemTypes = [elem.name for elem in elemTypes.copy()]
82
83 mesher = Mesher()
84
85 for e, elemType in enumerate(elemTypes):
86 times_N = []
87 wDef_N = []
88 dofs_N = []
89 zz1_N = []
90
91 # Loop over each mesh size (number of elements)
92 for N in list_N:
93 meshSize = b / N
94
95 # Define the domain for the mesh
96 domain = Domain(Point(), Point(x=L, y=h), meshSize)
97
98 # Generate the mesh using Gmsh
99 if dim == 2:
100 mesh = mesher.Mesh_2D(domain, [], elemType, isOrganised=isOrganised)
101 volume = mesh.area * material.thickness
102 else:
103 mesh = mesher.Mesh_Extrude(
104 domain,
105 [],
106 elemType=elemType,
107 extrude=[0, 0, b],
108 layers=[4],
109 isOrganised=isOrganised,
110 )
111 volume = mesh.volume
112 # Ensure that the volume matches the expected value (L * h * b)
113 assert np.abs(volume - (L * h * b)) / volume <= 1e-10
114
115 # Define nodes on the left boundary (x=0) and right boundary (x=L)
116 nodes_x0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
117 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
118
119 # Create or update the simulation object with the current mesh
120 if e == 0 and N == list_N[0]:
121 simu = Simulations.Elastic(mesh, material)
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.549 seconds)