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 : 2.749 s
Boundary Conditions : 2.066 ms
Matrix : 7.258 s
Solver : 2.166 m
Display : 2.060 s
PostProcessing : 565.774 ms
Resolution hyperelastic : 3.740 s
PyVista_Interface : 22.822 s
Paraview : 835.035 ms
Resolutions : 1.451 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.Results_Dir() + f"{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.Elastic.Isotropic(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.Elastic(mesh, material)
120 else:
121 simu.Bc_Init()
122 simu.mesh = mesh
123
124 # Set displacement boundary conditions
125 simu.add_dirichlet(nodes_x0, [0] * dim, simu.Get_unknowns())
126 # Set surface load on the right boundary (y-direction)
127 simu.add_surfLoad(nodes_xL, [-P / h / b], ["y"])
128
129 tic = Tic()
130
131 # Solve the simulation
132 simu.Solve()
133 simu.Save_Iter()
134
135 time = tic.Tac("Resolutions", "Total time", False)
136
137 # Get the computed deformation energy
138 Wdef = simu.Result("Wdef")
139
140 # Store the results for the current mesh size
141 times_N.append(time)
142 wDef_N.append(Wdef)
143 dofs_N.append(mesh.Nn * dim)
144 zz1_N.append(simu.Result("ZZ1"))
145
146 if elemType != mesh.elemType:
147 print("Error in mesh generation")
148
149 print(
150 f"Elem: {mesh.elemType}, nby: {N:2}, Wdef = {np.round(Wdef, 3)}, "
151 f"error = {np.abs(WdefRef - Wdef) / WdefRef:.2e}"
152 )
153
154 # Store the results for the current element type
155 times_elem_N.append(times_N)
156 wDef_elem_N.append(wDef_N)
157 dofs_elem_N.append(dofs_N)
158 zz1_elem_N.append(zz1_N)
159
160 # ----------------------------------------------
161 # Results
162 # ----------------------------------------------
163 # Display the convergence of deformation energy
164 ax_Wdef = Display.Init_Axes()
165 ax_error = Display.Init_Axes()
166 ax_times = Display.Init_Axes()
167 ax_zz1 = Display.Init_Axes()
168
169 print(f"\nWSA = {np.round(WdefRef, 4)} mJ")
170
171 for e, elemType in enumerate(elemTypes):
172 # Convergence of deformation energy
173 ax_Wdef.plot(dofs_elem_N[e], wDef_elem_N[e])
174
175 # Error in deformation energy
176 Wdef = np.array(wDef_elem_N[e])
177 error = (WdefRef - Wdef) / WdefRef * 100
178 ax_error.loglog(dofs_elem_N[e], error)
179
180 # Computation time
181 ax_times.loglog(dofs_elem_N[e], times_elem_N[e])
182 # ax_Times.plot(listDofs_e_nb[e], listTimes_e_nb[e])
183 # ax_Times.set_xscale('log')
184
185 # ZZ1
186 if elemType == elemTypes[0]:
187 last = ax_zz1.loglog(dofs_elem_N[e], error, label=f"{elemType}")
188 ax_zz1.loglog(
189 dofs_elem_N[e],
190 zz1_elem_N[e],
191 ls="--",
192 color=last[0]._color,
193 label=f"{elemType} (ZZ1)",
194 )
195
196 WdefRefArray = np.ones_like(dofs_elem_N[0]) * WdefRef
197 WdefRefArray5 = WdefRefArray * 0.95
198 # WdefRefArray5 = WdefRefArray * 1
199
200 # Deformation energy
201 ax_Wdef.grid()
202 ax_Wdef.set_xlim([-10, np.max(dofs_elem_N[0])])
203 ax_Wdef.set_xlabel("Degrees of Freedom (DOF)")
204 ax_Wdef.set_ylabel("Strain energy W [mJ]")
205 ax_Wdef.legend(elemTypes)
206 # ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color='red')
207 ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color="red")
208 plt.figure(ax_Wdef.figure)
209 Display.Save_fig(folder, "Energy")
210
211 # Error in deformation energy
212 ax_error.grid()
213 ax_error.set_xlabel("Degrees of Freedom (DOF)")
214 ax_error.set_ylabel("Error W [%]")
215 ax_error.legend(elemTypes)
216 plt.figure(ax_error.figure)
217 Display.Save_fig(folder, "Error")
218
219 # Error in deformation energy
220 ax_zz1.grid()
221 ax_zz1.set_xlabel("Degrees of Freedom (DOF)")
222 ax_zz1.set_ylabel("Error [%]")
223 ax_zz1.legend()
224 plt.figure(ax_zz1.figure)
225 Display.Save_fig(folder, "Error ZZ1")
226
227 # Computation time
228 ax_times.grid()
229 ax_times.set_xlabel("Degrees of Freedom (DOF)")
230 ax_times.set_ylabel("Computation Time [s]")
231 ax_times.legend(elemTypes)
232 plt.figure(ax_times.figure)
233 Display.Save_fig(folder, "Time")
234
235 # Plot the von Mises stress result using 20 color levels
236 Display.Plot_Result(simu, "Svm", ncolors=20)
237
238 if makeParaview:
239 # Generate Paraview files for visualization
240 Paraview.Save_simu(simu, folder, details=True)
241
242 # Show the total computation time
243 print()
244 Tic.Resume()
245
246 # Display the computation time history
247 # Tic.Plot_History(folder)
248
249 # Show all plots
250 plt.show()
Total running time of the script: (0 minutes 3.627 seconds)