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