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