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