Note
Go to the end to download the full example code.
Homog4#
Conduct 3d homogenization on a simple RVE.
C_hom =
6.776e+00 4.564e+00 4.385e+00 2.205e-16 -1.528e-14 -6.074e-06
4.564e+00 6.776e+00 4.385e+00 4.833e-15 -1.336e-15 5.675e-06
4.385e+00 4.385e+00 2.400e+01 2.155e-15 -4.516e-15 -1.151e-07
-2.564e-16 -1.292e-16 -2.067e-15 1.569e+00 -4.023e-08 -1.083e-17
3.016e-16 -1.885e-16 -4.984e-15 -4.023e-08 1.569e+00 -4.709e-17
-6.074e-06 5.675e-06 -1.151e-07 3.637e-16 1.104e-15 1.262e+00
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Models, Mesh, ElemType, Simulations, PyVista
17 from EasyFEA.Geoms import Points, Circle
18 from EasyFEA.FEM import LagrangeCondition, FeArray
19 from typing import Optional
20
21
22 def Compute_ukl(
23 simu: Simulations.Elastic,
24 Ekl: np.ndarray,
25 nodesKUBC: Optional[np.ndarray] = None,
26 pairedNodes: Optional[np.ndarray] = None,
27 pltSol=False,
28 ):
29 simu.Bc_Init()
30 mesh = simu.mesh
31 coord = mesh.coord
32
33 directions = ["x", "y", "z"]
34
35 if nodesKUBC is not None:
36
37 def func_ux(x, y, z):
38 return Ekl.dot([x, y, z])[0]
39
40 def func_uy(x, y, z):
41 return Ekl.dot([x, y, z])[1]
42
43 def func_uz(x, y, z):
44 return Ekl.dot([x, y, z])[2]
45
46 simu.add_dirichlet(nodesKUBC, [func_ux, func_uy, func_uz], directions)
47
48 if pairedNodes is not None:
49 for n0, n1 in pairedNodes:
50 nodes = np.array([n0, n1])
51
52 delta = coord[n0] - coord[n1]
53 values = Ekl @ delta
54
55 for d, direction in enumerate(directions):
56 dofs = simu.Bc_dofs_nodes(nodes, [direction])
57
58 condition = LagrangeCondition(
59 "elastic", nodes, dofs, [direction], [values[d]], [1, -1]
60 )
61 simu._Bc_Add_Lagrange(condition)
62
63 # remove rigid body motion
64 nodes = mesh.nodes
65 vect = np.ones(mesh.Nn) * 1 / mesh.Nn
66 for direction in directions:
67 # sum d_i / Nn = 0
68 dofs = simu.Bc_dofs_nodes(nodes, [direction])
69 condition = LagrangeCondition(
70 "elastic", nodes, dofs, [direction], [0], [vect]
71 )
72 simu._Bc_Add_Lagrange(condition)
73
74 ukl = simu.Solve()
75
76 simu.Save_Iter()
77
78 if pltSol:
79 PyVista.Plot(simu, "ux", deformFactor=0.3).show()
80 PyVista.Plot(simu, "uy", deformFactor=0.3).show()
81
82 PyVista.Plot(simu, "Sxx", deformFactor=0.3).show()
83 PyVista.Plot(simu, "Syy", deformFactor=0.3).show()
84 PyVista.Plot(simu, "Sxy", deformFactor=0.3).show()
85
86 return ukl
87
88
89 def Get_nodes(mesh: Mesh, tol=1e-6, plotSurfaces=False):
90 """Returns nodesXm, nodesXp, nodesYm, nodesYp, nodesZm, nodesZp"""
91
92 # PyVista.Plot_Tags(mesh).show()
93 coord = mesh.coord
94 nodesXm = mesh.Nodes_Conditions(lambda x, y, z: x <= coord[:, 0].min() + tol)
95 nodesXp = mesh.Nodes_Conditions(lambda x, y, z: x >= coord[:, 0].max() - tol)
96 nodesYm = mesh.Nodes_Conditions(lambda x, y, z: y <= coord[:, 1].min() + tol)
97 nodesYp = mesh.Nodes_Conditions(lambda x, y, z: y >= coord[:, 1].max() - tol)
98 nodesZm = mesh.Nodes_Conditions(lambda x, y, z: z <= coord[:, 2].min() + tol)
99 nodesZp = mesh.Nodes_Conditions(lambda x, y, z: z >= coord[:, 2].max() - tol)
100
101 if plotSurfaces:
102 plotter = PyVista.Plot(mesh, alpha=0.1)
103 colors = plt.get_cmap("tab10").colors
104
105 dict_nodes = {
106 f"{nodesXm=}".split("=")[0]: nodesXm,
107 f"{nodesXp=}".split("=")[0]: nodesXp,
108 f"{nodesYm=}".split("=")[0]: nodesYm,
109 f"{nodesYp=}".split("=")[0]: nodesYp,
110 f"{nodesZm=}".split("=")[0]: nodesZm,
111 f"{nodesZp=}".split("=")[0]: nodesZp,
112 }
113
114 for i, (name, nodes) in enumerate(dict_nodes.items()):
115 PyVista.Plot_Elements(
116 mesh, nodes, 2, color=colors[i], label=name, plotter=plotter
117 )
118 plotter.add_legend()
119 plotter.show()
120
121 return nodesXm, nodesXp, nodesYm, nodesYp, nodesZm, nodesZp
122
123
124 def Get_pairedNodes(
125 mesh: Mesh,
126 nodesXm: np.ndarray,
127 nodesXp: np.ndarray,
128 nodesYm: np.ndarray,
129 nodesYp: np.ndarray,
130 nodesZm: np.ndarray,
131 nodesZp: np.ndarray,
132 plotPBC=False,
133 ):
134
135 # remove Xp nodes on Y nodes
136 nodesYp = np.array(list(set(nodesYp) - set(nodesXp)))
137 nodesYm = np.array(list(set(nodesYm) - set(nodesXp)))
138
139 # remove Xp and Yp nodes on Z nodes
140 nodesZp = np.array(list(set(nodesZp) - set(nodesXp) - set(nodesYp)))
141 nodesZm = np.array(list(set(nodesZm) - set(nodesXp) - set(nodesYp)))
142
143 list_tuple_nodes = [
144 (nodesXm, nodesXp),
145 (nodesYm, nodesYp),
146 (nodesZm, nodesZp),
147 ]
148
149 # get paired nodes
150 for i, (nodes1, nodes2) in enumerate(list_tuple_nodes):
151
152 # sort nodes
153 sorted_nodes2 = np.zeros_like(nodes1)
154 coord = mesh.coord
155 for n, node in enumerate(nodes1):
156 dist = np.linalg.norm(coord[nodes2] - coord[node], axis=1)
157 sorted_nodes2[n] = nodes2[dist.argmin()]
158 nodes2 = sorted_nodes2
159
160 if plotPBC:
161 plotter = PyVista.Plot_Mesh(mesh, alpha=0.1)
162 PyVista.Plot_Nodes(mesh, nodes1, plotter=plotter)
163 PyVista.Plot_Nodes(mesh, nodes2, plotter=plotter)
164 plotter.add_lines(
165 np.concatenate(
166 (mesh.coord[nodes1], mesh.coord[nodes2]), axis=1
167 ).reshape(-1, 3),
168 color="k",
169 )
170 plotter.add_title(f"paired nodes along {['x','y','z'][i]} axis.")
171 plotter.show()
172
173 # set new paired nodes
174 newPairedNodes = np.array([nodes1, nodes2], dtype=int).T
175 if i == 0:
176 pairedNodes = newPairedNodes
177 else:
178 pairedNodes = np.concatenate((pairedNodes, newPairedNodes))
179
180 return pairedNodes
181
182
183 if __name__ == "__main__":
184 Display.Clear()
185
186 # ----------------------------------------------
187 # Configuration
188 # ----------------------------------------------
189
190 # use Periodic boundary conditions ?
191 usePBC = True
192 plotPBC = False
193 plotSurfaces = False
194
195 # ----------------------------------------------
196 # Mesh
197 # ----------------------------------------------
198 meshSize = 1 / 5
199
200 p0 = (-1 / 2, -1 / 2)
201 p1 = (1 / 2, -1 / 2)
202 p2 = (1 / 2, 1 / 2)
203 p3 = (-1 / 2, 1 / 2)
204 pts = [p0, p1, p2, p3]
205 contour = Points(pts, meshSize)
206
207 # inclusion
208 f = 0.4
209 r = 1 * np.sqrt(f / np.pi)
210 inclusion = Circle((0, 0), 2 * r, meshSize, isHollow=False)
211 # contour.Plot_Geoms([contour, inclusion])
212
213 elemType = ElemType.PRISM15
214 mesh = contour.Mesh_Extrude(
215 [inclusion],
216 [0, 0, 1],
217 [1 / meshSize],
218 elemType,
219 additionalPoints=[inclusion.center],
220 )
221 mesh.Translate(*-mesh.center) # center mesh on 0,0,0
222
223 plotter = PyVista.Plot_Mesh(mesh)
224 plotter.show_grid()
225 plotter.add_title("RVE")
226 plotter.show()
227
228 # ----------------------------------------------
229 # Get paired nodes
230 # ----------------------------------------------
231
232 tuple_nodes = Get_nodes(mesh, plotSurfaces=plotSurfaces)
233 if usePBC:
234 nodesKUBC = None
235 pairedNodes = Get_pairedNodes(mesh, *tuple_nodes, plotPBC=plotPBC)
236 else:
237 nodesKUBC = set(np.concatenate(tuple_nodes))
238 nodesKUBC = list(nodesKUBC)
239 pairedNodes = None
240
241 # ----------------------------------------------
242 # Material and Simulation
243 # ----------------------------------------------
244 elements_inclusion = (
245 np.array([]) if inclusion.isHollow else mesh.Elements_Tags(["V1"])
246 )
247 elements_matrix = mesh.Elements_Tags(["V0"])
248
249 E = np.zeros_like(mesh.groupElem.elements, dtype=float)
250 v = np.zeros_like(mesh.groupElem.elements, dtype=float)
251
252 E[elements_matrix] = 1 # MPa
253 v[elements_matrix] = 0.45
254
255 if elements_inclusion.size > 0:
256 E[elements_inclusion] = 50
257 v[elements_inclusion] = 0.3
258
259 material = Models.Elastic.Isotropic(3, E, v)
260
261 simu = Simulations.Elastic(mesh, material)
262
263 PyVista.Plot(simu, E, nodeValues=False, colorbarTitle="E [MPa]").show()
264 PyVista.Plot(simu, v, nodeValues=False, colorbarTitle="v").show()
265
266 # ----------------------------------------------
267 # Homogenization
268 # ----------------------------------------------
269 r2 = np.sqrt(2)
270 E1 = np.array(
271 [
272 [1, 0, 0],
273 [0, 0, 0],
274 [0, 0, 0],
275 ]
276 )
277 E2 = np.array(
278 [
279 [0, 0, 0],
280 [0, 1, 0],
281 [0, 0, 0],
282 ]
283 )
284 E3 = np.array(
285 [
286 [0, 0, 0],
287 [0, 0, 0],
288 [0, 0, 1],
289 ]
290 )
291 E12 = np.array(
292 [
293 [0, 1 / r2, 0],
294 [1 / r2, 0, 0],
295 [0, 0, 0],
296 ]
297 )
298 E13 = np.array(
299 [
300 [0, 0, 1 / r2],
301 [0, 0, 0],
302 [1 / r2, 0, 0],
303 ]
304 )
305 E23 = np.array(
306 [
307 [0, 0, 0],
308 [0, 0, 1 / r2],
309 [0, 1 / r2, 0],
310 ]
311 )
312
313 u11 = Compute_ukl(simu, E1, nodesKUBC, pairedNodes, True)
314 u22 = Compute_ukl(simu, E2, nodesKUBC, pairedNodes)
315 u33 = Compute_ukl(simu, E3, nodesKUBC, pairedNodes)
316 u12 = Compute_ukl(simu, E12, nodesKUBC, pairedNodes, True)
317 u13 = Compute_ukl(simu, E13, nodesKUBC, pairedNodes)
318 u23 = Compute_ukl(simu, E23, nodesKUBC, pairedNodes)
319
320 u11_e = mesh.Locates_sol_e(u11, asFeArray=True)
321 u22_e = mesh.Locates_sol_e(u22, asFeArray=True)
322 u33_e = mesh.Locates_sol_e(u33, asFeArray=True)
323 u12_e = mesh.Locates_sol_e(u12, asFeArray=True)
324 u13_e = mesh.Locates_sol_e(u13, asFeArray=True)
325 u23_e = mesh.Locates_sol_e(u23, asFeArray=True)
326
327 # ----------------------------------------------
328 # Effective elasticity tensor (C_hom)
329 # ----------------------------------------------
330 U_e = FeArray.zeros(*u11_e.shape, 6)
331
332 U_e[..., 0] = u11_e
333 U_e[..., 1] = u22_e
334 U_e[..., 2] = u33_e
335 U_e[..., 3] = u23_e
336 U_e[..., 4] = u13_e
337 U_e[..., 5] = u12_e
338
339 matrixType = "mass"
340 wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType)
341 B_e_pg = mesh.Get_B_e_pg(matrixType)
342
343 C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2])
344
345 # Be careful here you have to use all the volume even if there are holes
346 # if you use the mesh volume, multiply C_hom by the porosity (1-f)
347 C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / mesh.volume
348
349 if inclusion.isHollow:
350 C_hom *= 1 - f
351
352 formatted_array = ""
353 for i in range(6):
354 formatted_array += "\n"
355 for j in range(6):
356 formatted_array += f"{C_hom[i,j]:10.3e} "
357
358 print("C_hom =", formatted_array)
359
360 plt.show()
Total running time of the script: (0 minutes 4.647 seconds)