Note
Go to the end to download the full example code.
Homog4#
Conduct 3d homogenization. # Warning: Verify that the periodic boundary conditions have been correctly applied.
C_hom =
5.793e+00 3.569e+00 3.512e+00 3.532e-16 -5.546e-15 -1.519e-04
3.569e+00 5.857e+00 3.537e+00 3.884e-15 -7.533e-16 1.315e-04
3.512e+00 3.537e+00 2.300e+01 1.235e-15 -1.707e-15 -7.543e-06
5.712e-17 2.133e-16 6.665e-15 1.539e+00 3.430e-07 7.187e-18
2.110e-16 -7.865e-17 2.774e-14 3.430e-07 1.557e+00 -4.734e-17
-1.519e-04 1.315e-04 -7.543e-06 -3.702e-15 -1.614e-15 1.282e+00
13 from EasyFEA import Display, Models, plt, np, ElemType, Simulations, PyVista
14 from EasyFEA.Geoms import Points, Circle
15 from EasyFEA.fem import Mesh, LagrangeCondition, FeArray
16 from typing import Optional
17
18
19 def Compute_ukl(
20 simu: Simulations.ElasticSimu,
21 nodes_border: np.ndarray,
22 Ekl: np.ndarray,
23 paired_nodes: Optional[np.ndarray] = None,
24 pltSol=False,
25 ):
26 simu.Bc_Init()
27 mesh = simu.mesh
28 coord = mesh.coord
29
30 usePER = paired_nodes is not None
31
32 def func_ux(x, y, z):
33 return Ekl.dot([x, y, z])[0]
34
35 def func_uy(x, y, z):
36 return Ekl.dot([x, y, z])[1]
37
38 def func_uz(x, y, z):
39 return Ekl.dot([x, y, z])[2]
40
41 directions = ["x", "y", "z"]
42
43 simu.add_dirichlet(nodes_border, [func_ux, func_uy, func_uz], directions)
44
45 if usePER:
46 for n0, n1 in paired_nodes:
47 nodes = np.array([n0, n1])
48
49 for d, direction in enumerate(directions):
50 dofs = simu.Bc_dofs_nodes(nodes, [direction])
51
52 values = Ekl @ [
53 coord[n0, 0] - coord[n1, 0],
54 coord[n0, 1] - coord[n1, 1],
55 coord[n0, 2] - coord[n1, 2],
56 ]
57 value = values[d]
58
59 condition = LagrangeCondition(
60 "elastic", nodes, dofs, [direction], [value], [1, -1]
61 )
62 simu._Bc_Add_Lagrange(condition)
63
64 ukl = simu.Solve()
65
66 simu.Save_Iter()
67
68 if pltSol:
69 PyVista.Plot(simu, "ux", deformFactor=0.3).show()
70 PyVista.Plot(simu, "uy", deformFactor=0.3).show()
71
72 PyVista.Plot(simu, "Sxx", deformFactor=0.3).show()
73 PyVista.Plot(simu, "Syy", deformFactor=0.3).show()
74 PyVista.Plot(simu, "Sxy", deformFactor=0.3).show()
75
76 return ukl
77
78
79 def Get_nodes(mesh: Mesh, dimElem: int):
80 """Returns\n
81 nodesLeft, nodesRight, nodesUpper, nodesLower, nodesFront, nodesBack"""
82
83 coord = mesh.coord
84
85 conditions = {
86 "left": lambda x, y, z: x == coord[:, 0].min(),
87 "right": lambda x, y, z: x == coord[:, 0].max(),
88 "lower": lambda x, y, z: y == coord[:, 1].min(),
89 "upper": lambda x, y, z: y == coord[:, 1].max(),
90 "back": lambda x, y, z: z == coord[:, 2].min(),
91 "front": lambda x, y, z: z == coord[:, 2].max(),
92 }
93
94 # get nodes
95 nodes = {key: [] for key in conditions}
96 for groupElem in mesh.Get_list_groupElem(dimElem):
97 for key, condition in conditions.items():
98 nodes[key].extend(groupElem.Get_Nodes_Conditions(condition))
99
100 return tuple(np.asarray(nodes[key]) for key in conditions)
101
102
103 def Get_kubc_nodes(mesh: Mesh, dimElem: int = 0):
104 """Returns nodes_kubc"""
105
106 coord = mesh.coord
107
108 # conditions to get nodes on x and y faces
109 conditions = [
110 lambda x, y, z: x == coord[:, 0].min(),
111 lambda x, y, z: x == coord[:, 0].max(),
112 lambda x, y, z: y == coord[:, 1].min(),
113 lambda x, y, z: y == coord[:, 1].max(),
114 # lambda x, y, z: z == coord[:, 2].min(),
115 # lambda x, y, z: z == coord[:, 2].max(),
116 ]
117
118 # get nodes
119 nodes_kubc: set[int] = set()
120 for groupElem in mesh.Get_list_groupElem(dimElem):
121 for condition in conditions:
122 nodes = groupElem.Get_Nodes_Conditions(condition)
123 nodes_kubc = nodes_kubc.union(nodes)
124
125 return np.asarray(list(nodes_kubc), dtype=int)
126
127
128 if __name__ == "__main__":
129 Display.Clear()
130
131 # ----------------------------------------------
132 # Configuration
133 # ----------------------------------------------
134
135 # use Periodic boundary conditions ?
136 usePBC = True
137 plotPBC = True
138 plotSurfaces = False
139
140 # ----------------------------------------------
141 # Mesh
142 # ----------------------------------------------
143 meshSize = 1 / 10
144
145 p0 = (-1 / 2, -1 / 2)
146 p1 = (1 / 2, -1 / 2)
147 p2 = (1 / 2, 1 / 2)
148 p3 = (-1 / 2, 1 / 2)
149 pts = [p0, p1, p2, p3]
150 contour = Points(pts, meshSize)
151
152 # inclusion
153 f = 0.4
154 r = 1 * np.sqrt(f / np.pi)
155 inclusion = Circle((0, 0), 2 * r, meshSize, isHollow=False)
156 # contour.Plot_Geoms([contour, inclusion])
157
158 elemType = ElemType.PRISM6
159 mesh = contour.Mesh_Extrude([inclusion], [0, 0, 1], [1 / meshSize], elemType)
160
161 plotter = PyVista.Plot_Mesh(mesh)
162 plotter.add_title("RVE")
163 plotter.show()
164
165 # ----------------------------------------------
166 # Get paired nodes
167 # ----------------------------------------------
168
169 # PyVista.Plot_Tags(mesh).show()
170
171 nodesLeft, nodesRight, nodesLower, nodesUpper, nodesFront, nodesBack = Get_nodes(
172 mesh, dimElem=2
173 )
174
175 if plotSurfaces:
176 plotter = PyVista.Plot(mesh, alpha=0.1)
177 colors = plt.get_cmap("tab10").colors
178
179 dict_nodes = {
180 f"{nodesLeft=}".split("=")[0]: nodesLeft,
181 f"{nodesRight=}".split("=")[0]: nodesRight,
182 f"{nodesLower=}".split("=")[0]: nodesLower,
183 f"{nodesUpper=}".split("=")[0]: nodesUpper,
184 f"{nodesFront=}".split("=")[0]: nodesFront,
185 f"{nodesBack=}".split("=")[0]: nodesBack,
186 }
187
188 for i, (name, nodes) in enumerate(dict_nodes.items()):
189 PyVista.Plot_Elements(
190 mesh, nodes, 2, color=colors[i], label=name, plotter=plotter
191 )
192 plotter.add_legend()
193 plotter.show()
194
195 if usePBC:
196 # Hybrid PBC setup
197 nodes_kubc = Get_kubc_nodes(mesh, dimElem=0)
198 # PyVista.Plot_Nodes(mesh, nodes_kubc).show()
199
200 # get nodes along x, y and z directions
201 nodes_x = set(nodesLeft).union(nodesRight)
202 nodes_y = set(nodesLower).union(nodesUpper)
203 nodes_z = set(nodesFront).union(nodesBack)
204
205 # get nodes along ij directions
206 nodes_xy = nodes_x.union(nodes_y)
207 nodes_xz = nodes_x.union(nodes_z)
208 nodes_yz = nodes_y.union(nodes_z)
209
210 # Get nodes along the y-axis without nodes along x axis.
211 nodesLower = np.array([node for node in nodesLower if node not in nodes_x])
212 nodesUpper = np.array([node for node in nodesUpper if node not in nodes_x])
213
214 # Get nodes along the z-axis without nodes along x and y axis.
215 nodesFront = np.array([node for node in nodesFront if node not in nodes_xy])
216 nodesBack = np.array([node for node in nodesBack if node not in nodes_xy])
217
218 paired_surfaces = [
219 (nodesLeft, nodesRight), # x
220 (nodesLower, nodesUpper), # y
221 (nodesFront, nodesBack), # z
222 ]
223
224 # get paired nodes
225 for i, (nodes1, nodes2) in enumerate(paired_surfaces):
226
227 # remove kubc nodes
228 nodes1 = np.array(list(set(nodes1) - set(nodes_kubc)))
229 nodes2 = np.array(list(set(nodes2) - set(nodes_kubc)))
230
231 # sort nodes
232 sorted_nodes2 = np.zeros_like(nodes1)
233 coord = mesh.coord
234 for n, node in enumerate(nodes1):
235 dist = np.linalg.norm(coord[nodes2] - coord[node], axis=1)
236 sorted_nodes2[n] = nodes2[dist.argmin()]
237 nodes2 = sorted_nodes2
238
239 if plotPBC:
240 plotter = PyVista.Plot_Mesh(mesh, alpha=0.1)
241 PyVista.Plot_Nodes(mesh, nodes1, plotter=plotter)
242 PyVista.Plot_Nodes(mesh, nodes2, plotter=plotter)
243 plotter.add_lines(
244 np.concatenate(
245 (mesh.coord[nodes1], mesh.coord[nodes2]), axis=1
246 ).reshape(-1, 3),
247 color="k",
248 )
249 plotter.add_title(f"paired nodes along {['x','y','z'][i]} axis.")
250 plotter.show()
251
252 # set new paired nodes
253 newPairedNodes = np.array([nodes1, nodes2]).T
254 if i == 0:
255 paired_nodes = newPairedNodes
256 else:
257 paired_nodes = np.concatenate((paired_nodes, newPairedNodes))
258
259 else:
260 nodes_kubc = set(
261 np.concatenate(
262 [nodesLeft, nodesRight, nodesUpper, nodesLower, nodesBack, nodesFront]
263 )
264 )
265 nodes_kubc = list(nodes_kubc)
266 paired_nodes = None
267
268 # ----------------------------------------------
269 # Material and Simulation
270 # ----------------------------------------------
271 elements_inclusion = mesh.Elements_Tags(["V1"])
272 elements_matrix = mesh.Elements_Tags(["V0"])
273
274 E = np.zeros_like(mesh.groupElem.elements, dtype=float)
275 v = np.zeros_like(mesh.groupElem.elements, dtype=float)
276
277 E[elements_matrix] = 1 # MPa
278 v[elements_matrix] = 0.45
279
280 if elements_inclusion.size > 0:
281 E[elements_inclusion] = 50
282 v[elements_inclusion] = 0.3
283
284 material = Models.ElasIsot(3, E, v)
285
286 simu = Simulations.ElasticSimu(mesh, material)
287
288 PyVista.Plot(simu, E, nodeValues=False, colorbarTitle="E [MPa]").show()
289 PyVista.Plot(simu, v, nodeValues=False, colorbarTitle="v").show()
290
291 # ----------------------------------------------
292 # Homogenization
293 # ----------------------------------------------
294 r2 = np.sqrt(2)
295 E11 = np.array(
296 [
297 [1, 0, 0],
298 [0, 0, 0],
299 [0, 0, 0],
300 ]
301 )
302 E22 = np.array(
303 [
304 [0, 0, 0],
305 [0, 1, 0],
306 [0, 0, 0],
307 ]
308 )
309 E33 = np.array(
310 [
311 [0, 0, 0],
312 [0, 0, 0],
313 [0, 0, 1],
314 ]
315 )
316 E12 = np.array(
317 [
318 [0, 1 / r2, 0],
319 [1 / r2, 0, 0],
320 [0, 0, 0],
321 ]
322 )
323 E13 = np.array(
324 [
325 [0, 0, 1 / r2],
326 [0, 0, 0],
327 [1 / r2, 0, 0],
328 ]
329 )
330 E23 = np.array(
331 [
332 [0, 0, 0],
333 [0, 0, 1 / r2],
334 [0, 1 / r2, 0],
335 ]
336 )
337
338 u11 = Compute_ukl(simu, nodes_kubc, E11, paired_nodes, True)
339 u22 = Compute_ukl(simu, nodes_kubc, E22, paired_nodes)
340 u33 = Compute_ukl(simu, nodes_kubc, E33, paired_nodes)
341 u12 = Compute_ukl(simu, nodes_kubc, E12, paired_nodes, True)
342 u13 = Compute_ukl(simu, nodes_kubc, E13, paired_nodes)
343 u23 = Compute_ukl(simu, nodes_kubc, E23, paired_nodes)
344
345 u11_e = mesh.Locates_sol_e(u11, asFeArray=True)
346 u22_e = mesh.Locates_sol_e(u22, asFeArray=True)
347 u33_e = mesh.Locates_sol_e(u33, asFeArray=True)
348 u12_e = mesh.Locates_sol_e(u12, asFeArray=True)
349 u13_e = mesh.Locates_sol_e(u13, asFeArray=True)
350 u23_e = mesh.Locates_sol_e(u23, asFeArray=True)
351
352 # ----------------------------------------------
353 # Effective elasticity tensor (C_hom)
354 # ----------------------------------------------
355 U_e = FeArray.zeros(*u11_e.shape, 6)
356
357 U_e[..., 0] = u11_e
358 U_e[..., 1] = u22_e
359 U_e[..., 2] = u33_e
360 U_e[..., 3] = u23_e
361 U_e[..., 4] = u13_e
362 U_e[..., 5] = u12_e
363
364 matrixType = "mass"
365 wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType)
366 B_e_pg = mesh.Get_B_e_pg(matrixType)
367
368 C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2])
369
370 # Be careful here you have to use all the volume even if there are holes
371 # if you use the mesh volume, multiply C_hom by the porosity (1-f)
372 C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / mesh.volume
373
374 if inclusion.isHollow:
375 C_hom *= 1 - f
376
377 formatted_array = ""
378 for i in range(6):
379 formatted_array += "\n"
380 for j in range(6):
381 formatted_array += f"{C_hom[i,j]:10.3e} "
382
383 print("C_hom =", formatted_array)
384
385 plt.show()
Total running time of the script: (0 minutes 7.707 seconds)