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