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