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