Note
Go to the end to download the full example code.
MeshOptim1#
Mesh optimization using the ZZ1 criterion for a bending bracket.
error = 23.670 %
error = 10.163 %
error = 6.232 %
error = 3.046 %
error = 1.459 %
error = 0.434 %
Generate movie 1/6 (16.67 %) 1.63 s
Generate movie 2/6 (33.33 %) 825.63 ms
Generate movie 3/6 (50.00 %) 536.98 ms
Generate movie 4/6 (66.67 %) 376.18 ms
Generate movie 5/6 (83.33 %) 194.30 ms
Generate movie 6/6 (100.00 %) 0.00 µ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 Mesh,
23 Simulations,
24 Paraview,
25 PyVista,
26 )
27 from EasyFEA.Geoms import Point, Points
28 from EasyFEA.FEM import Calc_projector
29
30 if __name__ == "__main__":
31 Display.Clear()
32
33 # ----------------------------------------------
34 # Configuration
35 # ----------------------------------------------
36 dim = 2
37
38 # outputs
39 folder = Folder.Results_Dir()
40 plotProj = False
41 makeMovie = True
42 makeParaview = False
43
44 # geom
45 L = 120 # mm
46 h = L * 0.3
47 b = h
48
49 # load
50 P = 800 # N
51 lineLoad = P / h # N/mm
52 surfLoad = P / h / b # N/mm2
53
54 # criteria
55 threshold = (
56 1 / 100 if dim == 2 else 0.04
57 ) # Target error for the optimization process
58 iterMax = 20 # Maximum number of iterations
59 coef = 1 / 10 # Scaling coefficient for the optimization process
60
61 # ----------------------------------------------
62 # Mesh
63 # ----------------------------------------------
64 if dim == 2:
65 elemType = ElemType.TRI3
66 else:
67 elemType = ElemType.TETRA4
68 N = 2
69
70 pt1 = Point(isOpen=True, r=-10)
71 pt2 = Point(x=L)
72 pt3 = Point(x=L, y=h)
73 pt4 = Point(x=h, y=h, r=10)
74 pt5 = Point(x=h, y=L)
75 pt6 = Point(y=L)
76 pt7 = Point(x=h, y=h)
77
78 points = Points([pt1, pt2, pt3, pt4, pt5, pt6], h / N)
79
80 PyVista.Plot_Geoms(points.Get_Contour()).show()
81
82 def DoMesh(refineGeom=None) -> Mesh:
83 """Function used to generate the mesh."""
84 if dim == 2:
85 return Mesher().Mesh_2D(points, [], elemType, [], [refineGeom])
86 else:
87 return Mesher().Mesh_Extrude(
88 points, [], [0, 0, b], [], elemType, [], [refineGeom]
89 )
90
91 # Construct the initial mesh
92 mesh = DoMesh()
93 PyVista.Plot_Mesh(mesh).show()
94
95 # ----------------------------------------------
96 # Material and Simulation
97 # ----------------------------------------------
98 material = Models.Elastic.Isotropic(dim, E=210000, v=0.3, thickness=b)
99 simu = Simulations.Elastic(mesh, material)
100 simu.rho = 8100 * 1e-9
101
102 def DoSimu(refineGeom: str):
103 simu.mesh = DoMesh(refineGeom)
104
105 # get the nodes
106 nodes_Fixed = simu.mesh.Nodes_Conditions(lambda x, y, z: x == 0)
107 nodes_Load = simu.mesh.Nodes_Conditions(lambda x, y, z: x == L)
108
109 # do the simulation
110 simu.Bc_Init()
111 simu.add_dirichlet(
112 nodes_Fixed, [0] * dim, simu.Get_unknowns(), description="Fixed"
113 )
114 simu.add_surfLoad(nodes_Load, [-surfLoad], ["y"])
115
116 simu.Solve()
117
118 simu.Save_Iter()
119
120 return simu
121
122 simu = Simulations.Mesh_Optim_ZZ1(DoSimu, folder, threshold, iterMax, 1 / 10)
123 PyVista.Plot_BoundaryConditions(simu).show()
124
125 # ----------------------------------------------
126 # Results
127 # ----------------------------------------------
128 PyVista.Plot_Mesh(simu.mesh).show()
129 PyVista.Plot(simu, "ZZ1_e", nodeValues=False, nColors=11).show()
130
131 if plotProj:
132 simu.Set_Iter(0)
133 mesh0 = simu.mesh
134 u0 = np.reshape(simu.displacement, (mesh0.Nn, -1))
135
136 simu.Set_Iter(1)
137 mesh1 = simu.mesh
138
139 proj = Calc_projector(mesh0, mesh1)
140 uProj = np.zeros((mesh1.Nn, dim), dtype=float)
141 for d in range(dim):
142 uProj[:, d] = proj @ u0[:, d]
143
144 ax = Display.Plot_Result(
145 mesh0, np.linalg.norm(u0, axis=1), plotMesh=True, title="u0"
146 )
147 ax.plot(*mesh1.coord[:, :dim].T, ls="", marker="+", c="k", label="new nodes")
148 ax.legend()
149 Display.Plot_Result(
150 mesh1, np.linalg.norm(uProj, axis=1), plotMesh=True, title="uProj"
151 )
152
153 if makeParaview:
154 Paraview.Save_simu(simu, folder, nodeFields=["ZZ1_e"])
155
156 if makeMovie:
157
158 def func(plotter, n):
159 simu.Set_Iter(n)
160
161 PyVista.Plot_Mesh(simu, plotter=plotter)
162
163 zz1 = simu._Calc_ZZ1()[0]
164
165 plotter.add_title(f"ZZ1 = {zz1 * 100:.2f} %")
166
167 PyVista.Movie_func(func, len(simu.results), folder, "bracket.gif")
168
169 Tic.Plot_History()
170 plt.show()
Total running time of the script: (0 minutes 5.806 seconds)