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 -1.655e-15 1.123e-15 1.990e-06
4.564e+00 6.776e+00 4.385e+00 4.705e-16 -3.070e-15 -6.731e-06
4.385e+00 4.385e+00 2.400e+01 -2.263e-16 -5.685e-16 -1.367e-06
-1.090e-16 7.578e-16 6.303e-17 1.569e+00 -8.944e-08 1.438e-16
4.238e-17 -4.182e-16 -1.587e-14 -8.944e-08 1.569e+00 5.258e-17
1.990e-06 -6.731e-06 -1.367e-06 3.698e-16 1.535e-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, isFilled=True)
211 # contour.Plot_Geoms([contour, inclusion])
212
213 elemType = ElemType.PRISM15
214 mesh = contour.Mesh_Extrude([inclusion], [0, 0, 1], [1 / meshSize], elemType)
215 mesh.Translate(*-mesh.center) # center mesh on 0,0,0
216
217 plotter = PyVista.Plot_Mesh(mesh)
218 plotter.show_grid()
219 plotter.add_title("RVE")
220 plotter.show()
221
222 # ----------------------------------------------
223 # Get paired nodes
224 # ----------------------------------------------
225
226 tuple_nodes = Get_nodes(mesh, plotSurfaces=plotSurfaces)
227 if usePBC:
228 nodesKUBC = None
229 pairedNodes = Get_pairedNodes(mesh, *tuple_nodes, plotPBC=plotPBC)
230 else:
231 nodesKUBC = set(np.concatenate(tuple_nodes))
232 nodesKUBC = list(nodesKUBC)
233 pairedNodes = None
234
235 # ----------------------------------------------
236 # Material and Simulation
237 # ----------------------------------------------
238 elements_inclusion = (
239 np.array([]) if not inclusion.isFilled else mesh.Elements_Tags(["V1"])
240 )
241 elements_matrix = mesh.Elements_Tags(["V0"])
242
243 E = np.zeros_like(mesh.groupElem.elements, dtype=float)
244 v = np.zeros_like(mesh.groupElem.elements, dtype=float)
245
246 E[elements_matrix] = 1 # MPa
247 v[elements_matrix] = 0.45
248
249 if elements_inclusion.size > 0:
250 E[elements_inclusion] = 50
251 v[elements_inclusion] = 0.3
252
253 material = Models.Elastic.Isotropic(3, E, v)
254
255 simu = Simulations.Elastic(mesh, material)
256
257 PyVista.Plot(simu, E, nodeValues=False, colorbarTitle="E [MPa]").show()
258 PyVista.Plot(simu, v, nodeValues=False, colorbarTitle="v").show()
259
260 # ----------------------------------------------
261 # Homogenization
262 # ----------------------------------------------
263 r2 = np.sqrt(2)
264 E1 = np.array(
265 [
266 [1, 0, 0],
267 [0, 0, 0],
268 [0, 0, 0],
269 ]
270 )
271 E2 = np.array(
272 [
273 [0, 0, 0],
274 [0, 1, 0],
275 [0, 0, 0],
276 ]
277 )
278 E3 = np.array(
279 [
280 [0, 0, 0],
281 [0, 0, 0],
282 [0, 0, 1],
283 ]
284 )
285 E12 = np.array(
286 [
287 [0, 1 / r2, 0],
288 [1 / r2, 0, 0],
289 [0, 0, 0],
290 ]
291 )
292 E13 = np.array(
293 [
294 [0, 0, 1 / r2],
295 [0, 0, 0],
296 [1 / r2, 0, 0],
297 ]
298 )
299 E23 = np.array(
300 [
301 [0, 0, 0],
302 [0, 0, 1 / r2],
303 [0, 1 / r2, 0],
304 ]
305 )
306
307 u11 = Compute_ukl(simu, E1, nodesKUBC, pairedNodes, True)
308 u22 = Compute_ukl(simu, E2, nodesKUBC, pairedNodes)
309 u33 = Compute_ukl(simu, E3, nodesKUBC, pairedNodes)
310 u12 = Compute_ukl(simu, E12, nodesKUBC, pairedNodes, True)
311 u13 = Compute_ukl(simu, E13, nodesKUBC, pairedNodes)
312 u23 = Compute_ukl(simu, E23, nodesKUBC, pairedNodes)
313
314 u11_e = mesh.Locates_sol_e(u11, asFeArray=True)
315 u22_e = mesh.Locates_sol_e(u22, asFeArray=True)
316 u33_e = mesh.Locates_sol_e(u33, asFeArray=True)
317 u12_e = mesh.Locates_sol_e(u12, asFeArray=True)
318 u13_e = mesh.Locates_sol_e(u13, asFeArray=True)
319 u23_e = mesh.Locates_sol_e(u23, asFeArray=True)
320
321 # ----------------------------------------------
322 # Effective elasticity tensor (C_hom)
323 # ----------------------------------------------
324 U_e = FeArray.zeros(*u11_e.shape, 6)
325
326 U_e[..., 0] = u11_e
327 U_e[..., 1] = u22_e
328 U_e[..., 2] = u33_e
329 U_e[..., 3] = u23_e
330 U_e[..., 4] = u13_e
331 U_e[..., 5] = u12_e
332
333 matrixType = "mass"
334 wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType)
335 B_e_pg = mesh.Get_B_e_pg(matrixType)
336
337 C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2])
338
339 xMin, yMin, zMin = mesh.coord.min(axis=0)
340 xMax, yMax, zMax = mesh.coord.max(axis=0)
341 volume = (xMax - xMin) * (yMax - yMin) * (zMax - zMin)
342
343 C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / volume
344
345 formatted_array = ""
346 for i in range(6):
347 formatted_array += "\n"
348 for j in range(6):
349 formatted_array += f"{C_hom[i,j]:10.3e} "
350
351 print("C_hom =", formatted_array)
352
353 plt.show()
Total running time of the script: (0 minutes 7.755 seconds)