{ "cells": [ { "cell_type": "markdown", "id": "1cf5ba7b-acbd-4019-96ba-f5c4679c41f8", "metadata": {}, "source": [ "# Two inclusions mechanical analysis" ] }, { "cell_type": "markdown", "id": "cdd08d7d-4080-4ebc-89d5-b37d37407963", "metadata": {}, "source": [ "Here we present a study case of a thick plate with 2 inclusions. One softer and the other stiffer than the base material. Then we compute the energy deformation on one inclusion.\n" ] }, { "cell_type": "markdown", "id": "f9510b12-57b2-489f-b2a4-32dcb96d8c83", "metadata": {}, "source": [ "## Includes" ] }, { "cell_type": "code", "execution_count": null, "id": "5449da9b-d2fb-4e19-9c39-91eb09ac39e2", "metadata": {}, "outputs": [], "source": [ "from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym\n", "from Muscat.FE.KR.KRBlock import KRBlockVector\n", "from Muscat.IO.GmshReader import ReadGmsh\n", "from Muscat.FE.Fields.FieldTools import GetPointRepresentation\n", "from Muscat.FE.Fields.FEField import FEField\n", "from Muscat.FE.SymPhysics import MechPhysics\n", "import Muscat.FE.SymWeakForm as SWF\n", "from Muscat.FE.MaterialHelp import HookeIso\n", "from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter\n", "from Muscat.FE.DofNumbering import ComputeDofNumbering\n", "from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0\n", "from Muscat.FE.Integration import IntegrateGeneral\n", "from Muscat.IO.XdmfWriter import XdmfWriter\n", "from Muscat.MeshTools.MeshInspectionTools import ExtractElementsByElementFilter\n", "from Muscat.FE.Fields.FieldTools import VectorToFEFieldsData\n", "from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory\n", "\n", "import numpy as np\n", "import pyvista\n", "pyvista.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots\n", "\n", "import pathlib" ] }, { "cell_type": "markdown", "id": "27f510c9-54a2-4aa0-94e1-f2d796441e4f", "metadata": {}, "source": [ "## Unstructured FEA with mechanical " ] }, { "cell_type": "code", "execution_count": null, "id": "bab328da-7c0e-4c2b-a09c-14f456b4e4f3", "metadata": {}, "outputs": [], "source": [ "problem = UnstructuredFeaSym()\n", "\n", "# the mechanical problem\n", "mechPhysics = MechPhysics()\n", "\n", "# Definition of the degree of the spaces [1 or 2]\n", "mechPhysics.SetSpaceToLagrange(P=1)\n", "\n", "# add weak form terms to the tanget matrix\n", "youngModulusInclusionBase = 1.0\n", "mechPhysics.AddBFormulation( \"Bulk\",mechPhysics.GetBulkFormulation(youngModulusInclusionBase,0.3) )\n", "mechPhysics.AddBFormulation( \"Inclusion1\",mechPhysics.GetBulkFormulation(5*youngModulusInclusionBase,0.3) )\n", "youngModulusInclusionII = youngModulusInclusionBase/2\n", "mechPhysics.AddBFormulation( \"Inclusion2\",mechPhysics.GetBulkFormulation(youngModulusInclusionII,0.3) )\n", "\n", "# add weak form term to the rhs\n", "mechPhysics.AddLFormulation( \"Right\", mechPhysics.GetForceFormulation([1,0,0],1) )\n", "problem.physics.append(mechPhysics)" ] }, { "cell_type": "markdown", "id": "84c86dcc-dac8-431f-8ef7-681af1d2f4ec", "metadata": {}, "source": [ "## Boundary conditions" ] }, { "cell_type": "code", "execution_count": null, "id": "3b4869f3-53e7-4164-8fff-9659b296adc4", "metadata": {}, "outputs": [], "source": [ "# the boundary conditions\n", "dirichlet = KRBlockVector()\n", "dirichlet.AddArg(\"u\").On('Left').Fix0().Fix1().Fix2().To()\n", "\n", "problem.solver.constraints.AddConstraint(dirichlet)" ] }, { "cell_type": "markdown", "id": "3800838f-8873-407b-9af2-07c8deae1770", "metadata": {}, "source": [ "## Numerical part (mesh)" ] }, { "cell_type": "code", "execution_count": null, "id": "5fa151bd-7b7b-46a8-aca2-02a561f4bbb7", "metadata": {}, "outputs": [], "source": [ "# Read The mesh\n", "\n", "mesh = ReadGmsh(pathlib.Path(\"\").resolve() / pathlib.Path(\"TwoInclusions.msh\"))\n", "problem.SetMesh(mesh)\n", "\n", "print(mesh)" ] }, { "cell_type": "markdown", "id": "371827b4-6155-4eff-96a7-1b031438d467", "metadata": {}, "source": [ "## Resolution" ] }, { "cell_type": "code", "execution_count": null, "id": "9449b1ee-a3e7-466d-b1c6-21da947464b8", "metadata": {}, "outputs": [], "source": [ "# we compute the numbering (maping from the mesh to the linear system)\n", "problem.ComputeDofNumbering()\n", "\n", "from Muscat.Helpers.Timer import Timer\n", "with Timer(\"Assembly \"):\n", " k,f = problem.GetLinearProblem()\n", "\n", "#compute the constraints to add to the system\n", "problem.ComputeConstraintsEquations()\n", "\n", "\n", "with Timer(\"Solve\"):\n", " problem.Solve(k,f)\n", "\n", "problem.PushSolutionVectorToUnknownFields()\n", "\n", "# Recover a point representation of the displacement\n", "mesh.nodeFields[\"sol\"] = GetPointRepresentation(problem.unknownFields)\n", "print(mesh)" ] }, { "cell_type": "markdown", "id": "181ead61-426b-4ef2-9737-1783227a9f08", "metadata": {}, "source": [ "## Post Processing " ] }, { "cell_type": "code", "execution_count": null, "id": "3a8e6731-6cac-4940-afbc-e890185be49d", "metadata": {}, "outputs": [], "source": [ "#Creation of a fake fields to export the rhs member\n", "rhsFields = [ FEField(mesh=mesh,space=None,numbering=problem.numberings[i]) for i in range(3) ]\n", "VectorToFEFieldsData(f,rhsFields)\n", "problem.mesh.nodeFields[\"RHS\"] = GetPointRepresentation(rhsFields)" ] }, { "cell_type": "markdown", "id": "029ab6e9-623c-4871-8f62-453b6ea57496", "metadata": {}, "source": [ "### Compute of the strain energy only on the second inclusion (integral in each element)" ] }, { "cell_type": "code", "execution_count": null, "id": "e792f38a-ac5d-462f-bf5f-dadda4eed0eb", "metadata": {}, "outputs": [], "source": [ "symdep = SWF.GetField(\"u\",3,3)\n", "K = HookeIso(youngModulusInclusionII,0.3,dim=3)\n", "symCellData = SWF.GetField(\"cellData\",1,3)\n", "symCellDataT = SWF.GetTestField(\"cellData\",1,3)\n", "\n", "print(\"Post process\")\n", "\n", "EnerForm = SWF.ToVoigtEpsilon(SWF.Strain(symdep)).T*K*SWF.ToVoigtEpsilon(SWF.Strain(symdep))*symCellDataT\n", "\n", "print(\"Post process Eval\")\n", "ff = ElementFilter(eTag=\"Inclusion2\")\n", "\n", "p0Numbering = ComputeDofNumbering(mesh,LagrangeSpaceP0)\n", "\n", "energyDensityField = FEField(name=\"cellData\",mesh=problem.mesh,numbering=p0Numbering,space=LagrangeSpaceP0)\n", "\n", "m,energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm, constants={},\n", " fields=problem.unknownFields, unknownFields = [ energyDensityField ], elementFilter=ff)\n", "print(\"energyDensity\",energyDensity)\n", "energyDensityField.data = energyDensity\n", "\n", "problem.mesh.elemFields[\"Energy\"] = energyDensity\n", "print(\"Strain energy on the second inclusion:\", np.sum(energyDensity) )" ] }, { "cell_type": "markdown", "id": "c0d95438-795d-40f2-b68e-fafbb468d085", "metadata": {}, "source": [ "## Graphical Output " ] }, { "cell_type": "code", "execution_count": null, "id": "ca649100-b447-4214-b58e-51dc0a1c3392", "metadata": {}, "outputs": [], "source": [ "# generate a element field with un index per material\n", "mesh.elemFields[\"inclussions\"] = np.zeros(mesh.GetNumberOfElements())\n", "for i, name in enumerate((\"Inclusion1\",\"Inclusion2\")):\n", " for selection in ElementFilter(eTag=name)(mesh):\n", " mesh.elemFields[\"inclussions\"][selection.GetMeshSlice()]= 1 + i\n", "\n", "import pyvista as pv\n", "from Muscat.Bridges.PyVistaBridge import MeshToPyVista\n", "\n", "plotter = pv.Plotter(shape=(2, 2))\n", "\n", "camera = [(1.4464393507963245, 1.760895123462273, 2.9110300316336843),\n", " (-0.22723828652575934, -0.08454810909441325, 0.21480857439407788),\n", " (-0.25495941474725625, 0.8643967380870753, -0.43337509851877826)]\n", "\n", "plotter.subplot(0, 0)\n", "plotter.add_text(\"Two Inclusions\", font_size=10)\n", "plotter.add_mesh(MeshToPyVista(mesh), scalars=\"inclussions\", show_edges=True, show_scalar_bar=False, edge_color=\"grey\", cmap=\"rainbow_r\")\n", "plotter.camera_position = camera\n", "\n", "plotter.subplot(0, 1)\n", "plotter.add_text(\"Deformation\", font_size=10)\n", "deformedMesh = mesh.View()\n", "deformedMesh.nodes = mesh.nodes + problem.mesh.nodeFields[\"sol\"]*.3\n", "deformedMesh.nodeFields[\"sol\"] = problem.mesh.nodeFields[\"sol\"]\n", "plotter.add_mesh(MeshToPyVista(deformedMesh), scalars=\"sol\", show_edges=False, show_scalar_bar=False, cmap=\"coolwarm_r\")\n", "plotter.camera_position = camera\n", "rightMesh = ExtractElementsByElementFilter(deformedMesh,ElementFilter(eTag=\"Right\"))\n", "rightMesh.nodeFields = deformedMesh.nodeFields\n", "glyphs = MeshToPyVista(rightMesh).glyph(orient=\"RHS\", scale=\"RHS\", factor=10, geom=pv.Arrow())\n", "plotter.add_mesh(glyphs, show_scalar_bar=False)\n", "\n", "plotter.subplot(1, 0)\n", "plotter.add_text(\"Inclusion 2 energy dencity\", font_size=10)\n", "plotter.add_mesh(MeshToPyVista(ExtractElementsByElementFilter(mesh,ElementFilter(eTag=\"Bulk\")) ).extract_feature_edges(45) )\n", "in1Mesh = ExtractElementsByElementFilter(mesh,ElementFilter(eTag=\"Inclusion2\"))\n", "in1Mesh.elemFields[\"Energy\"] = mesh.elemFields[\"Energy\"][in1Mesh.GetElementsOriginalIDs()]\n", "plotter.add_mesh(MeshToPyVista(in1Mesh), scalars=\"Energy\", style=\"surface\" , cmap=\"coolwarm\" )\n", "plotter.camera_position = camera\n", "\n", "plotter.subplot(1, 1)\n", "plotter.add_text(\"Inclusion deformation\", font_size=10)\n", "deformedMesh = mesh.View()\n", "deformedMesh.nodes = mesh.nodes + problem.mesh.nodeFields[\"sol\"]*.3\n", "buldkDeformedMesh = ExtractElementsByElementFilter(deformedMesh,ElementFilter(eTag=\"Bulk\"))\n", "buldkDeformedMesh.nodeFields[\"sol\"] = problem.mesh.nodeFields[\"sol\"]\n", "plotter.add_mesh(MeshToPyVista(buldkDeformedMesh).extract_feature_edges(45),scalars=\"sol\", show_edges=False, show_scalar_bar=False)\n", "in1Mesh.nodes = mesh.nodes + problem.mesh.nodeFields[\"sol\"]*.3\n", "plotter.add_mesh(MeshToPyVista(in1Mesh), scalars=\"Energy\", style=\"surface\" , cmap=\"coolwarm\" )\n", "plotter.add_mesh(glyphs)\n", "plotter.camera_position = camera\n", "\n", "plotter.show()" ] }, { "cell_type": "markdown", "id": "c7da2233-1e5b-4736-8825-1ba4f97e3853", "metadata": {}, "source": [ "## Export data to file" ] }, { "cell_type": "code", "execution_count": null, "id": "ee1474af-47b7-46c1-a634-c150f5e53600", "metadata": {}, "outputs": [], "source": [ "# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)\n", "writer = XdmfWriter(TemporaryDirectory.GetTempPath() + 'TwoInclusions_Output.xdmf')\n", "writer.SetHdf5(False)\n", "writer.Open()\n", "writer.Write(mesh,PointFields=list(mesh.nodeFields.values()), PointFieldsNames=list(mesh.nodeFields.keys()),\n", " CellFields=list(mesh.elemFields.values()), CellFieldsNames=list(mesh.elemFields.keys()))\n", "writer.Close()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.10" } }, "nbformat": 4, "nbformat_minor": 5 }