Note
Go to the end to download the full example code.
MeshOptim3#
Mesh optimization using the ZZ1 criterion for a letter weigher.
error = 80.031 %
error = 63.003 %
error = 34.483 %
error = 19.716 %
error = 9.033 %
error = 2.937 %
error = 1.279 %
error = 0.694 %
Generate movie 1/8 (12.50 %) 884.93 ms
Generate movie 2/8 (25.00 %) 675.95 ms
Generate movie 3/8 (37.50 %) 548.39 ms
Generate movie 4/8 (50.00 %) 447.77 ms
Generate movie 5/8 (62.50 %) 355.53 ms
Generate movie 6/8 (75.00 %) 264.46 ms
Generate movie 7/8 (87.50 %) 159.34 ms
Generate movie 8/8 (100.00 %) 0.00 µs
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 Mesh,
24 Simulations,
25 Paraview,
26 PyVista,
27 )
28 from EasyFEA.Geoms import Point, Points
29 from EasyFEA.FEM import Calc_projector
30
31 if __name__ == "__main__":
32 Display.Clear()
33
34 # ----------------------------------------------
35 # Configuration
36 # ----------------------------------------------
37 dim = 2
38
39 # outputs
40 folder = Folder.Results_Dir()
41 plotProj = False
42 makeMovie = True
43 makeParaview = False
44
45 # geom
46 L = 80
47 h1 = L / 4
48 e1 = h1 * 0.1
49 h2 = (h1 - e1) * 0.95
50 e2 = (h1 - e1 - h2) / 2
51 r = h2 / 4
52 l = L / 2
53 b = h1
54
55 # load
56 P = 800 # N
57 lineLoad = P / h1 # N/mm
58 surfLoad = P / h1 / b # N/mm2
59
60 # criteria
61 threshold = (
62 1 / 100 if dim == 2 else 0.04
63 ) # Target error for the optimization process
64 iterMax = 20 # Maximum number of iterations
65 coef = 1 / 10 # Scaling coefficient for the optimization process
66
67 # ----------------------------------------------
68 # Mesh
69 # ----------------------------------------------
70 if dim == 2:
71 elemType = ElemType.TRI3
72 else:
73 elemType = ElemType.TETRA4
74 meshSize = r / 3
75
76 pt1 = Point()
77 pt2 = Point(L - h1)
78 pt3 = pt2 + [e1, -e1]
79 pt4 = Point(L, -e1)
80 pt5 = pt4 + [0, h1]
81 pt6 = Point(h1, h1 - e1)
82 pt7 = pt6 + [-e1, e1]
83 pt8 = pt1 + [0, h1]
84
85 points = Points([pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8], meshSize)
86
87 p1 = Point(L / 2 - l / 2, e2, r=r)
88 p2 = Point(L / 2 - l / 2 + 2 * r, e2, r=r)
89 p3 = Point(L / 2 - l / 2 + 2 * r, e2 + r)
90 p4 = p3.copy()
91 p4.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
92 p5 = p2.copy()
93 p5.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
94 p6 = p1.copy()
95 p6.Symmetry((L / 2, (h1 - e1) / 2), (1, 0))
96 p7 = p6.copy()
97 p7.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
98 p8 = p5.copy()
99 p8.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
100 p9 = p4.copy()
101 p9.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
102 p10 = p3.copy()
103 p10.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
104 p11 = p2.copy()
105 p11.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
106 p12 = p1.copy()
107 p12.Symmetry((L / 2, (h1 - e1) / 2), (0, 1))
108
109 inclusion = Points(
110 [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], meshSize, True
111 )
112 inclusions = [inclusion]
113
114 PyVista.Plot_Geoms([points.Get_Contour(), inclusion.Get_Contour()]).show()
115
116 def DoMesh(refineGeom=None) -> Mesh:
117 """Function used to generate the mesh."""
118 if dim == 2:
119 return Mesher().Mesh_2D(points, inclusions, elemType, [], [refineGeom])
120 else:
121 return Mesher().Mesh_Extrude(
122 points, inclusions, [0, 0, b], [], elemType, [], [refineGeom]
123 )
124
125 # Construct the initial mesh
126 mesh = DoMesh()
127 PyVista.Plot_Mesh(mesh).show()
128
129 # ----------------------------------------------
130 # Material and Simulation
131 # ----------------------------------------------
132 material = Models.Elastic.Isotropic(dim, E=210000, v=0.3, thickness=b)
133 simu = Simulations.Elastic(mesh, material)
134 simu.rho = 8100 * 1e-9
135
136 def DoSimu(refineGeom: str):
137 simu.mesh = DoMesh(refineGeom)
138
139 # get the nodes
140 nodes_Fixed = simu.mesh.Nodes_Conditions(lambda x, y, z: x == 0)
141 nodes_Load = simu.mesh.Nodes_Conditions(lambda x, y, z: x == L)
142
143 # do the simulation
144 simu.Bc_Init()
145 simu.add_dirichlet(
146 nodes_Fixed, [0] * dim, simu.Get_unknowns(), description="Fixed"
147 )
148 simu.add_surfLoad(nodes_Load, [-surfLoad], ["y"])
149
150 simu.Solve()
151
152 simu.Save_Iter()
153
154 return simu
155
156 simu = Simulations.Mesh_Optim_ZZ1(DoSimu, folder, threshold, iterMax, 1 / 10)
157 PyVista.Plot_BoundaryConditions(simu).show()
158
159 # ----------------------------------------------
160 # Plot
161 # ----------------------------------------------
162 PyVista.Plot_Mesh(simu.mesh).show()
163 PyVista.Plot(simu, "ZZ1_e", nodeValues=False, nColors=11).show()
164
165 if plotProj:
166 simu.Set_Iter(0)
167 mesh0 = simu.mesh
168 u0 = np.reshape(simu.displacement, (mesh0.Nn, -1))
169
170 simu.Set_Iter(1)
171 mesh1 = simu.mesh
172
173 proj = Calc_projector(mesh0, mesh1)
174 uProj = np.zeros((mesh1.Nn, dim), dtype=float)
175 for d in range(dim):
176 uProj[:, d] = proj @ u0[:, d]
177
178 ax = Display.Plot_Result(
179 mesh0, np.linalg.norm(u0, axis=1), plotMesh=True, title="u0"
180 )
181 ax.plot(*mesh1.coord[:, :dim].T, ls="", marker="+", c="k", label="new nodes")
182 ax.legend()
183 Display.Plot_Result(
184 mesh1, np.linalg.norm(uProj, axis=1), plotMesh=True, title="uProj"
185 )
186
187 if makeParaview:
188 Paraview.Save_simu(simu, folder, nodeFields=["ZZ1_e"])
189
190 if makeMovie:
191
192 def func(plotter, n):
193 simu.Set_Iter(n)
194
195 PyVista.Plot_Mesh(simu, plotter=plotter)
196
197 zz1 = simu._Calc_ZZ1()[0]
198
199 plotter.add_title(f"ZZ1 = {zz1 * 100:.2f} %")
200
201 PyVista.Movie_func(func, len(simu.results), folder, "letterWeigher.gif")
202
203 Tic.Plot_History()
204 plt.show()
Total running time of the script: (0 minutes 9.989 seconds)