.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/LinearizedElasticity/Homog4.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_LinearizedElasticity_Homog4.py: Homog4 ====== Conduct 3d homogenization. WARNING ------- Verify that the periodic boundary conditions have been correctly applied. .. GENERATED FROM PYTHON SOURCE LINES 16-393 .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_001.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_001.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_001.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_002.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_002.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_003.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_003.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_003.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_004.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_004.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_004.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_005.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_005.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_005.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_006.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_006.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_006.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_007.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_007.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_007.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_008.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_008.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_008.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_009.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_009.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_009.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_010.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_010.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_010.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_011.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_011.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_011.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_012.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_012.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_012.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_013.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_013.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_013.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_014.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_014.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_014.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_015.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_015.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_015.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/LinearizedElasticity/images/sphx_glr_Homog4_016.png :alt: Homog4 :srcset: /examples/LinearizedElasticity/images/sphx_glr_Homog4_016.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.7.2/docs/examples/LinearizedElasticity/images/sphx_glr_Homog4_016.vtksz .. rst-class:: sphx-glr-script-out .. code-block:: none 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 | .. code-block:: Python :lineno-start: 17 import matplotlib.pyplot as plt import numpy as np from EasyFEA import Display, Models, ElemType, Simulations, PyVista from EasyFEA.Geoms import Points, Circle from EasyFEA.FEM import Mesh, LagrangeCondition, FeArray from typing import Optional def Compute_ukl( simu: Simulations.Elastic, nodes_border: np.ndarray, Ekl: np.ndarray, paired_nodes: Optional[np.ndarray] = None, pltSol=False, ): simu.Bc_Init() mesh = simu.mesh coord = mesh.coord usePER = paired_nodes is not None def func_ux(x, y, z): return Ekl.dot([x, y, z])[0] def func_uy(x, y, z): return Ekl.dot([x, y, z])[1] def func_uz(x, y, z): return Ekl.dot([x, y, z])[2] directions = ["x", "y", "z"] simu.add_dirichlet(nodes_border, [func_ux, func_uy, func_uz], directions) if usePER: for n0, n1 in paired_nodes: nodes = np.array([n0, n1]) for d, direction in enumerate(directions): dofs = simu.Bc_dofs_nodes(nodes, [direction]) values = Ekl @ [ coord[n0, 0] - coord[n1, 0], coord[n0, 1] - coord[n1, 1], coord[n0, 2] - coord[n1, 2], ] value = values[d] condition = LagrangeCondition( "elastic", nodes, dofs, [direction], [value], [1, -1] ) simu._Bc_Add_Lagrange(condition) ukl = simu.Solve() simu.Save_Iter() if pltSol: PyVista.Plot(simu, "ux", deformFactor=0.3).show() PyVista.Plot(simu, "uy", deformFactor=0.3).show() PyVista.Plot(simu, "Sxx", deformFactor=0.3).show() PyVista.Plot(simu, "Syy", deformFactor=0.3).show() PyVista.Plot(simu, "Sxy", deformFactor=0.3).show() return ukl def Get_nodes(mesh: Mesh, dimElem: int): """Returns\n nodesLeft, nodesRight, nodesUpper, nodesLower, nodesFront, nodesBack""" coord = mesh.coord conditions = { "left": lambda x, y, z: x == coord[:, 0].min(), "right": lambda x, y, z: x == coord[:, 0].max(), "lower": lambda x, y, z: y == coord[:, 1].min(), "upper": lambda x, y, z: y == coord[:, 1].max(), "back": lambda x, y, z: z == coord[:, 2].min(), "front": lambda x, y, z: z == coord[:, 2].max(), } # get nodes nodes = {key: [] for key in conditions} for groupElem in mesh.Get_list_groupElem(dimElem): for key, condition in conditions.items(): nodes[key].extend(groupElem.Get_Nodes_Conditions(condition)) return tuple(np.asarray(nodes[key]) for key in conditions) def Get_kubc_nodes(mesh: Mesh, dimElem: int = 0): """Returns nodes_kubc""" coord = mesh.coord # conditions to get nodes on x and y faces conditions = [ lambda x, y, z: x == coord[:, 0].min(), lambda x, y, z: x == coord[:, 0].max(), lambda x, y, z: y == coord[:, 1].min(), lambda x, y, z: y == coord[:, 1].max(), # lambda x, y, z: z == coord[:, 2].min(), # lambda x, y, z: z == coord[:, 2].max(), ] # get nodes nodes_kubc: set[int] = set() for groupElem in mesh.Get_list_groupElem(dimElem): for condition in conditions: nodes = groupElem.Get_Nodes_Conditions(condition) nodes_kubc = nodes_kubc.union(nodes) return np.asarray(list(nodes_kubc), dtype=int) if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- # use Periodic boundary conditions ? usePBC = True plotPBC = True plotSurfaces = False # ---------------------------------------------- # Mesh # ---------------------------------------------- meshSize = 1 / 10 p0 = (-1 / 2, -1 / 2) p1 = (1 / 2, -1 / 2) p2 = (1 / 2, 1 / 2) p3 = (-1 / 2, 1 / 2) pts = [p0, p1, p2, p3] contour = Points(pts, meshSize) # inclusion f = 0.4 r = 1 * np.sqrt(f / np.pi) inclusion = Circle((0, 0), 2 * r, meshSize, isHollow=False) # contour.Plot_Geoms([contour, inclusion]) elemType = ElemType.PRISM6 mesh = contour.Mesh_Extrude([inclusion], [0, 0, 1], [1 / meshSize], elemType) plotter = PyVista.Plot_Mesh(mesh) plotter.add_title("RVE") plotter.show() # ---------------------------------------------- # Get paired nodes # ---------------------------------------------- # PyVista.Plot_Tags(mesh).show() nodesLeft, nodesRight, nodesLower, nodesUpper, nodesFront, nodesBack = Get_nodes( mesh, dimElem=2 ) if plotSurfaces: plotter = PyVista.Plot(mesh, alpha=0.1) colors = plt.get_cmap("tab10").colors dict_nodes = { f"{nodesLeft=}".split("=")[0]: nodesLeft, f"{nodesRight=}".split("=")[0]: nodesRight, f"{nodesLower=}".split("=")[0]: nodesLower, f"{nodesUpper=}".split("=")[0]: nodesUpper, f"{nodesFront=}".split("=")[0]: nodesFront, f"{nodesBack=}".split("=")[0]: nodesBack, } for i, (name, nodes) in enumerate(dict_nodes.items()): PyVista.Plot_Elements( mesh, nodes, 2, color=colors[i], label=name, plotter=plotter ) plotter.add_legend() plotter.show() if usePBC: # Hybrid PBC setup nodes_kubc = Get_kubc_nodes(mesh, dimElem=0) # PyVista.Plot_Nodes(mesh, nodes_kubc).show() # get nodes along x, y and z directions nodes_x = set(nodesLeft).union(nodesRight) nodes_y = set(nodesLower).union(nodesUpper) nodes_z = set(nodesFront).union(nodesBack) # get nodes along ij directions nodes_xy = nodes_x.union(nodes_y) nodes_xz = nodes_x.union(nodes_z) nodes_yz = nodes_y.union(nodes_z) # Get nodes along the y-axis without nodes along x axis. nodesLower = np.array([node for node in nodesLower if node not in nodes_x]) nodesUpper = np.array([node for node in nodesUpper if node not in nodes_x]) # Get nodes along the z-axis without nodes along x and y axis. nodesFront = np.array([node for node in nodesFront if node not in nodes_xy]) nodesBack = np.array([node for node in nodesBack if node not in nodes_xy]) paired_surfaces = [ (nodesLeft, nodesRight), # x (nodesLower, nodesUpper), # y (nodesFront, nodesBack), # z ] # get paired nodes for i, (nodes1, nodes2) in enumerate(paired_surfaces): # remove kubc nodes nodes1 = np.array(list(set(nodes1) - set(nodes_kubc))) nodes2 = np.array(list(set(nodes2) - set(nodes_kubc))) # sort nodes sorted_nodes2 = np.zeros_like(nodes1) coord = mesh.coord for n, node in enumerate(nodes1): dist = np.linalg.norm(coord[nodes2] - coord[node], axis=1) sorted_nodes2[n] = nodes2[dist.argmin()] nodes2 = sorted_nodes2 if plotPBC: plotter = PyVista.Plot_Mesh(mesh, alpha=0.1) PyVista.Plot_Nodes(mesh, nodes1, plotter=plotter) PyVista.Plot_Nodes(mesh, nodes2, plotter=plotter) plotter.add_lines( np.concatenate( (mesh.coord[nodes1], mesh.coord[nodes2]), axis=1 ).reshape(-1, 3), color="k", ) plotter.add_title(f"paired nodes along {['x','y','z'][i]} axis.") plotter.show() # set new paired nodes newPairedNodes = np.array([nodes1, nodes2]).T if i == 0: paired_nodes = newPairedNodes else: paired_nodes = np.concatenate((paired_nodes, newPairedNodes)) else: nodes_kubc = set( np.concatenate( [nodesLeft, nodesRight, nodesUpper, nodesLower, nodesBack, nodesFront] ) ) nodes_kubc = list(nodes_kubc) paired_nodes = None # ---------------------------------------------- # Material and Simulation # ---------------------------------------------- elements_inclusion = mesh.Elements_Tags(["V1"]) elements_matrix = mesh.Elements_Tags(["V0"]) E = np.zeros_like(mesh.groupElem.elements, dtype=float) v = np.zeros_like(mesh.groupElem.elements, dtype=float) E[elements_matrix] = 1 # MPa v[elements_matrix] = 0.45 if elements_inclusion.size > 0: E[elements_inclusion] = 50 v[elements_inclusion] = 0.3 material = Models.Elastic.Isotropic(3, E, v) simu = Simulations.Elastic(mesh, material) PyVista.Plot(simu, E, nodeValues=False, colorbarTitle="E [MPa]").show() PyVista.Plot(simu, v, nodeValues=False, colorbarTitle="v").show() # ---------------------------------------------- # Homogenization # ---------------------------------------------- r2 = np.sqrt(2) E11 = np.array( [ [1, 0, 0], [0, 0, 0], [0, 0, 0], ] ) E22 = np.array( [ [0, 0, 0], [0, 1, 0], [0, 0, 0], ] ) E33 = np.array( [ [0, 0, 0], [0, 0, 0], [0, 0, 1], ] ) E12 = np.array( [ [0, 1 / r2, 0], [1 / r2, 0, 0], [0, 0, 0], ] ) E13 = np.array( [ [0, 0, 1 / r2], [0, 0, 0], [1 / r2, 0, 0], ] ) E23 = np.array( [ [0, 0, 0], [0, 0, 1 / r2], [0, 1 / r2, 0], ] ) u11 = Compute_ukl(simu, nodes_kubc, E11, paired_nodes, True) u22 = Compute_ukl(simu, nodes_kubc, E22, paired_nodes) u33 = Compute_ukl(simu, nodes_kubc, E33, paired_nodes) u12 = Compute_ukl(simu, nodes_kubc, E12, paired_nodes, True) u13 = Compute_ukl(simu, nodes_kubc, E13, paired_nodes) u23 = Compute_ukl(simu, nodes_kubc, E23, paired_nodes) u11_e = mesh.Locates_sol_e(u11, asFeArray=True) u22_e = mesh.Locates_sol_e(u22, asFeArray=True) u33_e = mesh.Locates_sol_e(u33, asFeArray=True) u12_e = mesh.Locates_sol_e(u12, asFeArray=True) u13_e = mesh.Locates_sol_e(u13, asFeArray=True) u23_e = mesh.Locates_sol_e(u23, asFeArray=True) # ---------------------------------------------- # Effective elasticity tensor (C_hom) # ---------------------------------------------- U_e = FeArray.zeros(*u11_e.shape, 6) U_e[..., 0] = u11_e U_e[..., 1] = u22_e U_e[..., 2] = u33_e U_e[..., 3] = u23_e U_e[..., 4] = u13_e U_e[..., 5] = u12_e matrixType = "mass" wJ_e_pg = mesh.Get_weightedJacobian_e_pg(matrixType) B_e_pg = mesh.Get_B_e_pg(matrixType) C_Mat = Models.Reshape_variable(material.C, *B_e_pg.shape[:2]) # Be careful here you have to use all the volume even if there are holes # if you use the mesh volume, multiply C_hom by the porosity (1-f) C_hom = (wJ_e_pg * C_Mat @ B_e_pg @ U_e).sum((0, 1)) / mesh.volume if inclusion.isHollow: C_hom *= 1 - f formatted_array = "" for i in range(6): formatted_array += "\n" for j in range(6): formatted_array += f"{C_hom[i,j]:10.3e} " print("C_hom =", formatted_array) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.152 seconds) .. _sphx_glr_download_examples_LinearizedElasticity_Homog4.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Homog4.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Homog4.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Homog4.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_