.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/LinearizedElasticity/MeshConvergence.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_LinearizedElasticity_MeshConvergence.py: MeshConvergence =============== Verification of energy convergence for a bending beam for all available elements. .. GENERATED FROM PYTHON SOURCE LINES 12-241 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_001.png :alt: MeshConvergence :srcset: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_001.png :class: sphx-glr-multi-img * .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_002.png :alt: MeshConvergence :srcset: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_003.png :alt: MeshConvergence :srcset: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_004.png :alt: MeshConvergence :srcset: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_005.png :alt: $\sigma_{vm}$ :srcset: /examples/LinearizedElasticity/images/sphx_glr_MeshConvergence_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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.143 s Boundary Conditions: 1.036 ms Matrix: 1.644 s Solver: 1.004 s Resolutions: 2.173 s PostProcessing: 278.500 ms Display: 1.166 s | .. code-block:: Python :lineno-start: 13 import matplotlib.pyplot as plt import numpy as np from EasyFEA import Display, Folder, Models, Tic, ElemType, Simulations, Paraview from EasyFEA.Geoms import Domain, Point if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- dim = 2 # Define the dimension of the problem (2D or 3D) # outputs folder = Folder.Results_Dir() + f"{dim}D" plotResult = True makeParaview = False # geom L = 120 # mm h = 13 # Height b = 13 # Width # model E = 210000 # MPa (Young's modulus) v = 0.25 # Poisson's ratio material = Models.Elastic.Isotropic(dim, thickness=b, E=E, v=v, planeStress=True) # load P = 800 # N # expected energy WdefRef = 2 * P**2 * L / E / h / b * (L**2 / h / b + (1 + v) * 3 / 5) # ---------------------------------------------- # Mesh # ---------------------------------------------- isOrganised = True # List of mesh sizes (number of elements) to investigate convergence if dim == 2: list_N = np.arange(1, 10, 1) else: list_N = np.arange(1, 8, 2) # Lists to store data for plotting times_elem_N = [] # times for element type and N size wDef_elem_N = [] # energy dofs_elem_N = [] # dofs zz1_elem_N = [] # zz1 # ---------------------------------------------- # Simulations # ---------------------------------------------- # Loop over each element type for both 2D and 3D simulations elemTypes = ElemType.Get_2D()[:] if dim == 2 else ElemType.Get_3D() # elemTypes = [elem.name for elem in elemTypes.copy()] for e, elemType in enumerate(elemTypes): times_N = [] wDef_N = [] dofs_N = [] zz1_N = [] # Loop over each mesh size (number of elements) for N in list_N: meshSize = b / N # Define the domain for the mesh domain = Domain(Point(), Point(x=L, y=h), meshSize) # Generate the mesh using Gmsh if dim == 2: mesh = domain.Mesh_2D([], elemType, isOrganised=isOrganised) volume = mesh.area * material.thickness else: mesh = domain.Mesh_Extrude( [], elemType=elemType, extrude=[0, 0, b], layers=[4], isOrganised=isOrganised, ) volume = mesh.volume # Ensure that the volume matches the expected value (L * h * b) assert np.abs(volume - (L * h * b)) / volume <= 1e-10 # Define nodes on the left boundary (x=0) and right boundary (x=L) nodes_x0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0) nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L) # Create or update the simulation object with the current mesh if e == 0 and N == list_N[0]: simu = Simulations.Elastic(mesh, material) else: simu.Bc_Init() simu.mesh = mesh # Set displacement boundary conditions simu.add_dirichlet(nodes_x0, [0] * dim, simu.Get_unknowns()) # Set surface load on the right boundary (y-direction) simu.add_surfLoad(nodes_xL, [-P / h / b], ["y"]) tic = Tic() # Solve the simulation simu.Solve() simu.Save_Iter() time = tic.Tac("Resolutions", "Total time", False) # Get the computed deformation energy Wdef = simu.Result("Wdef") # Store the results for the current mesh size times_N.append(time) wDef_N.append(Wdef) dofs_N.append(mesh.Nn * dim) zz1_N.append(simu.Result("ZZ1")) if elemType != mesh.elemType: print("Error in mesh generation") print( f"Elem: {mesh.elemType}, nby: {N:2}, Wdef = {np.round(Wdef, 3)}, " f"error = {np.abs(WdefRef - Wdef) / WdefRef:.2e}" ) # Store the results for the current element type times_elem_N.append(times_N) wDef_elem_N.append(wDef_N) dofs_elem_N.append(dofs_N) zz1_elem_N.append(zz1_N) # ---------------------------------------------- # Results # ---------------------------------------------- # Display the convergence of deformation energy ax_Wdef = Display.Init_Axes() ax_error = Display.Init_Axes() ax_times = Display.Init_Axes() ax_zz1 = Display.Init_Axes() print(f"\nWSA = {np.round(WdefRef, 4)} mJ") for e, elemType in enumerate(elemTypes): # Convergence of deformation energy ax_Wdef.plot(dofs_elem_N[e], wDef_elem_N[e]) # Error in deformation energy Wdef = np.array(wDef_elem_N[e]) error = (WdefRef - Wdef) / WdefRef * 100 ax_error.loglog(dofs_elem_N[e], error) # Computation time ax_times.loglog(dofs_elem_N[e], times_elem_N[e]) # ax_Times.plot(listDofs_e_nb[e], listTimes_e_nb[e]) # ax_Times.set_xscale('log') # ZZ1 if elemType == elemTypes[0]: last = ax_zz1.loglog(dofs_elem_N[e], error, label=f"{elemType}") ax_zz1.loglog( dofs_elem_N[e], zz1_elem_N[e], ls="--", color=last[0]._color, label=f"{elemType} (ZZ1)", ) WdefRefArray = np.ones_like(dofs_elem_N[0]) * WdefRef WdefRefArray5 = WdefRefArray * 0.95 # WdefRefArray5 = WdefRefArray * 1 # Deformation energy ax_Wdef.grid() ax_Wdef.set_xlim([-10, np.max(dofs_elem_N[0])]) ax_Wdef.set_xlabel("Degrees of Freedom (DOF)") ax_Wdef.set_ylabel("Strain energy W [mJ]") ax_Wdef.legend(elemTypes) # ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color='red') ax_Wdef.fill_between(dofs_N, WdefRefArray, WdefRefArray5, alpha=0.5, color="red") plt.figure(ax_Wdef.figure) Display.Save_fig(folder, "Energy") # Error in deformation energy ax_error.grid() ax_error.set_xlabel("Degrees of Freedom (DOF)") ax_error.set_ylabel("Error W [%]") ax_error.legend(elemTypes) plt.figure(ax_error.figure) Display.Save_fig(folder, "Error") # Error in deformation energy ax_zz1.grid() ax_zz1.set_xlabel("Degrees of Freedom (DOF)") ax_zz1.set_ylabel("Error [%]") ax_zz1.legend() plt.figure(ax_zz1.figure) Display.Save_fig(folder, "Error ZZ1") # Computation time ax_times.grid() ax_times.set_xlabel("Degrees of Freedom (DOF)") ax_times.set_ylabel("Computation Time [s]") ax_times.legend(elemTypes) plt.figure(ax_times.figure) Display.Save_fig(folder, "Time") # Plot the von Mises stress result using 20 color levels Display.Plot_Result(simu, "Svm", ncolors=20) if makeParaview: # Generate Paraview files for visualization Paraview.Save_simu(simu, folder, details=True) # Show the total computation time print() Tic.Resume() # Display the computation time history # Tic.Plot_History(folder) # Show all plots plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.435 seconds) .. _sphx_glr_download_examples_LinearizedElasticity_MeshConvergence.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: MeshConvergence.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: MeshConvergence.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: MeshConvergence.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_