{ "cells": [ { "cell_type": "markdown", "id": "a69faa87-ef81-4a54-8179-205bb045bb7c", "metadata": {}, "source": [ "# Mesh adaption\n", "\n", "In this section, we propose to illustrate on a 2d example the various tools available in Muscat for mesh adaption." ] }, { "cell_type": "markdown", "id": "b92e0301-6648-4dda-b093-f1c81e4cdf7a", "metadata": {}, "source": [ "## Modeling the problem\n", "\n", "Consider the pierced plate shown in the following figure (length 4 m, width 2m):\n", "\n", "\n", "\n", "We want to calculate the plate deformation when it is subjected to a tension of F = 1 GPa at both ends, as depicted in the figure. \n", "The plate is assumed to be made of a homogeneous, isotropic linear elastic material with Poisson ratio ν = 0.3 and Young modulus E = 10 GPa.\n", "By virtue of symmetry, we can solve the problem on any quarter of the plate, imposing zero transverse displacements along the symmetry axes, shown in blue. \n", "\n", "## Numerical solution using Muscat\n", "\n", "Muscat's integrated finite element solver is used to solve the problem.\n", "\n", "First we initialize a new variational problem." ] }, { "cell_type": "code", "execution_count": null, "id": "acbfcf41", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym\n", "\n", "problem = UnstructuredFeaSym(spaceDim=2)" ] }, { "cell_type": "markdown", "id": "8e1baad0", "metadata": {}, "source": [ "Then we choose a linear elastic formulation with isoparametric finite elements" ] }, { "cell_type": "code", "execution_count": null, "id": "ee94418c", "metadata": {}, "outputs": [], "source": [ "from Muscat.FE.SymPhysics import MechPhysics\n", "mechPhysics = MechPhysics(dim=2)\n", "mechPhysics.planeStress = True\n", "mechPhysics.SetSpaceToLagrange(isoParam=True)" ] }, { "cell_type": "markdown", "id": "a4241a63", "metadata": {}, "source": [ "and we define the material properties and the loading." ] }, { "cell_type": "code", "execution_count": null, "id": "0b0b48b5", "metadata": {}, "outputs": [], "source": [ "mechPhysics.AddBFormulation(\"bulkA\", mechPhysics.GetBulkFormulation(10., 0.3))\n", "mechPhysics.AddBFormulation(\"2D\", mechPhysics.GetBulkFormulation(10., 0.3))\n", "mechPhysics.AddLFormulation(\"Y1\", mechPhysics.GetPressureFormulation(1))\n", "\n", "problem.physics.append(mechPhysics)" ] }, { "cell_type": "markdown", "id": "cf68d089", "metadata": {}, "source": [ "We impose the Dirichlet boundary conditions." ] }, { "cell_type": "code", "execution_count": null, "id": "93f583b5", "metadata": {}, "outputs": [], "source": [ "from Muscat.FE.KR.KRBlock import KRBlockVector\n", "\n", "dirichlet_lock_x = KRBlockVector().AddArg(\"u\").On('x0').On('X0').BlockDirection0()\n", "problem.solver.constraints.AddConstraint(dirichlet_lock_x)\n", "\n", "dirichlet_lock_y = KRBlockVector().AddArg(\"u\").On('y0').BlockDirection1()\n", "problem.solver.constraints.AddConstraint(dirichlet_lock_y)" ] }, { "cell_type": "markdown", "id": "a9ad6f34", "metadata": {}, "source": [ "We read the mesh and link it to the problem." ] }, { "cell_type": "code", "execution_count": null, "id": "1a032bd0", "metadata": {}, "outputs": [], "source": [ "from Muscat.IO.XdmfReader import ReadXdmf\n", "from Muscat.TestData import GetTestDataPath\n", "\n", "mesh = ReadXdmf(GetTestDataPath()+\"holeplatequad.xdmf\")\n", "\n", "print(\"Initial\", mesh)\n", "\n", "problem.SetMesh(mesh)\n", "problem.ComputeDofNumbering()" ] }, { "cell_type": "markdown", "id": "3b82c388", "metadata": {}, "source": [ "Finally we solve the variational formulation" ] }, { "cell_type": "code", "execution_count": null, "id": "13bcc663", "metadata": {}, "outputs": [], "source": [ "k,f = problem.GetLinearProblem()\n", "problem.ComputeConstraintsEquations()\n", "problem.Solve(k,f)\n", "problem.PushSolutionVectorToUnknownFields()\n", "\n", "from Muscat.FE.Fields.FieldTools import GetPointRepresentation\n", "\n", "sol = GetPointRepresentation(problem.unknownFields)\n", "mesh.nodeFields[\"solution\"] = sol" ] }, { "cell_type": "markdown", "id": "233e5f9d", "metadata": {}, "source": [ "and plot the solution." ] }, { "cell_type": "code", "execution_count": null, "id": "d7825ae7", "metadata": {}, "outputs": [], "source": [ "import pyvista as pv\n", "pv.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots\n", "from Muscat.Bridges.PyVistaBridge import MeshToPyVista\n", "\n", "pyVistaMesh = MeshToPyVista(mesh)\n", "edges = pyVistaMesh.separate_cells().extract_feature_edges()\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, scalars=\"solution\", scalar_bar_args=dict(vertical=True))\n", "pl.add_mesh(edges, color=\"black\")\n", "pl.show(cpos= \"xy\")" ] }, { "cell_type": "markdown", "id": "0fb1b677-520f-4d1c-89cf-1bf231ea97d4", "metadata": {}, "source": [ "## Mesh adaption\n", "\n", "### Linear simplex meshes\n", "\n", "There is no reliable remesher for meshes other than linear simplex meshes, i.e. those composed of points, 2-point bars, 3-point triangles and 4-point tetrahedrons. Muscat offers tools for transforming any mesh (with pyramids, prisms, quads, hexa, etc.) into a simplex mesh, while guaranteeing the transfer of tags and associated fields.\n", "\n", "First we can linearize the mesh:" ] }, { "cell_type": "code", "execution_count": null, "id": "1610e8b3", "metadata": {}, "outputs": [], "source": [ "from Muscat.MeshTools.MeshTetrahedrization import FromQuadToLinear\n", "\n", "mesh = FromQuadToLinear(mesh)\n", "\n", "print(\"Linear\", mesh)" ] }, { "cell_type": "markdown", "id": "54386643", "metadata": {}, "source": [ "Then we can convert it to a simplex mesh:" ] }, { "cell_type": "code", "execution_count": null, "id": "ee90fcdc", "metadata": {}, "outputs": [], "source": [ "from Muscat.MeshTools.MeshTetrahedrization import Tetrahedrization\n", "\n", "mesh = Tetrahedrization(mesh)\n", "\n", "print(\"Linear simplex\", mesh)" ] }, { "cell_type": "markdown", "id": "649739e7", "metadata": {}, "source": [ "And plot the resulting mesh with the solution." ] }, { "cell_type": "code", "execution_count": null, "id": "1497cf84", "metadata": {}, "outputs": [], "source": [ "pyVistaMesh = MeshToPyVista(mesh)\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, scalars=\"solution\", show_edges=True, scalar_bar_args=dict(vertical=True))\n", "pl.show(cpos= \"xy\")" ] }, { "cell_type": "markdown", "id": "da4a44dd", "metadata": {}, "source": [ "### Remeshing from a mesh size map\n", "\n", "The first step is to select a criterion on which to adapt the mesh. This criterion can be any scalar field defined at mesh nodes. Here, for example, we will choose the displacement norm, and we want to obtain a mesh which guarantees that the relative error committed on this quantity of interest is lower than 1e-5 at all points.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fff1722b", "metadata": {}, "outputs": [], "source": [ "psi = np.linalg.norm(mesh.nodeFields[\"solution\"], axis=1)\n", "err = 1e-5" ] }, { "cell_type": "markdown", "id": "3f5ba62d", "metadata": {}, "source": [ "The second step is to compute thank to Muscat the ideal mesh size associated with this criterion. " ] }, { "cell_type": "code", "execution_count": null, "id": "72aa353d", "metadata": {}, "outputs": [], "source": [ "from Muscat.MeshTools.AnisotropicMetricComputation import ComputeMetric\n", "\n", "siz = ComputeMetric(mesh, psi, err=err, isotropic=True) \n", "\n", "mesh.nodeFields[\"size\"] = siz\n", "\n", "pyVistaMesh = MeshToPyVista(mesh)\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, scalars=\"size\", show_edges=True, scalar_bar_args=dict(vertical=True), clim=[0.01, 0.1])\n", "pl.show(cpos= \"xy\")" ] }, { "cell_type": "markdown", "id": "4ddd8a90", "metadata": {}, "source": [ "The last step consists in using a remesher to build the desired mesh from the mesh size map." ] }, { "cell_type": "code", "execution_count": null, "id": "40d22fec-ad0c-429f-936b-8e4467ccaae6", "metadata": {}, "outputs": [], "source": [ "from Muscat.MeshTools.Remesh import Remesh\n", "\n", "iso_mesh = Remesh(mesh, metric=siz)\n", "\n", "print(\"Isotropic\", iso_mesh)\n", "\n", "pyVistaMesh = MeshToPyVista(iso_mesh)\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, show_edges=True)\n", "pl.show(cpos= \"xy\")" ] }, { "cell_type": "markdown", "id": "87f46ab3", "metadata": {}, "source": [ "### Remeshing from a metric computation\n", "\n", "From the same criterion, it is also possible to compute an anisotropic metric, that gives at each node the principal directions and the associated ideal mesh sizes." ] }, { "cell_type": "code", "execution_count": null, "id": "e3a690da", "metadata": {}, "outputs": [], "source": [ "met = ComputeMetric(mesh, psi, err=err) " ] }, { "cell_type": "markdown", "id": "da6a12f4", "metadata": {}, "source": [ "The remesher will then produce an anisotropic mesh." ] }, { "cell_type": "code", "execution_count": null, "id": "8d8abda0", "metadata": {}, "outputs": [], "source": [ "aniso_mesh = Remesh(mesh, metric=met)\n", "\n", "print(\"Anisotropic\", aniso_mesh)\n", "\n", "pyVistaMesh = MeshToPyVista(aniso_mesh)\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, show_edges=True)\n", "pl.show(cpos= \"xy\")" ] }, { "cell_type": "markdown", "id": "38cc2c48", "metadata": {}, "source": [ "### Intersection of metrics\n", "\n", "In Muscat, it is possible to remesh with respect to several criteria simultaneously. In the following example, we request a mesh that is not only optimal with respect to the previously calculated solution, but also with the analytical function |y - x - 2|. \n", "\n", "We first compute the metric associated with the analytical function." ] }, { "cell_type": "code", "execution_count": null, "id": "265eadea", "metadata": {}, "outputs": [], "source": [ "psi2 = np.abs(mesh.nodes[:, 1]-mesh.nodes[:, 0]-2)\n", "met2 = ComputeMetric(mesh, psi2, err=1e-4) " ] }, { "cell_type": "markdown", "id": "08bb19f3", "metadata": {}, "source": [ "And we intersect it with the previous metric." ] }, { "cell_type": "code", "execution_count": null, "id": "356964b3", "metadata": {}, "outputs": [], "source": [ "from Muscat.MeshTools.AnisotropicMetricComputation import IntersectionOfMetrics\n", "\n", "met3 = IntersectionOfMetrics(metrics=[met, met2])" ] }, { "cell_type": "markdown", "id": "f618653d", "metadata": {}, "source": [ "We can use this new metric to obtain a mesh that is adapted to both criteria." ] }, { "cell_type": "code", "execution_count": null, "id": "1e2c8cdc", "metadata": {}, "outputs": [], "source": [ "newMesh = Remesh(mesh, metric=met3)\n", "\n", "print(\"Intersected\", newMesh)\n", "\n", "pyVistaMesh = MeshToPyVista(newMesh)\n", "pl = pv.Plotter()\n", "pl.add_mesh(pyVistaMesh, show_edges=True)\n", "pl.show(cpos=\"xy\")" ] } ], "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.11" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "0ad6b6964bc94e04b178fa6bdff4306e": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "0ca0741127ca4f00b281d6209585f3c5": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "115ff364d54f48459a3e5cdd933221bc": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_402619bb773a4638bea12afc9ed31ddd", "style": "IPY_MODEL_2e7f521310b64f9c84e8b61472f9de29", "value": "" } }, "1274e3d6d52f4f5da650b7e12c48cddf": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "29ad252eecf94b5b8cd96b2d98d10f1c": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "2c2cf398cf9a4e5cad4ddb41ab052a8d": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "2e7f521310b64f9c84e8b61472f9de29": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "3081904fc8af4fd588dfc219e0e79ff3": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "402619bb773a4638bea12afc9ed31ddd": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "41e8004cc50547f2ad6b232850da0248": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "4784afb52eea4934a781147b02853520": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "4847f791c69e42979bae07ea1f897c87": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "4c71aa0209124e4a9aac93fa1b11064b": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_9133b0a085f34c35824ed61ebcf02d9a", "style": "IPY_MODEL_a0f0220c23ca41529695926531aa1326", "value": "" } }, "664f24b62c0a4603a10079dd6af11035": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "6fa299cad8784b438ecbaca9cdf6a9ab": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_4784afb52eea4934a781147b02853520", "style": "IPY_MODEL_3081904fc8af4fd588dfc219e0e79ff3", "value": "" } }, "8ea3971d0131494b9e908dae0a4743d8": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_1274e3d6d52f4f5da650b7e12c48cddf", "style": "IPY_MODEL_0ca0741127ca4f00b281d6209585f3c5", "value": "" } }, "911a3cf0298d41c3945a8654545ae98f": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_a5f7b0acc4464e9fad7666f38f22d505", "style": "IPY_MODEL_edbcaba1de234543b4eae485fa286675", "value": "" } }, "9133b0a085f34c35824ed61ebcf02d9a": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "97deadeca22b4908afd91708a8cf2749": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_2c2cf398cf9a4e5cad4ddb41ab052a8d", "style": "IPY_MODEL_664f24b62c0a4603a10079dd6af11035", "value": "" } }, "a0f0220c23ca41529695926531aa1326": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "a5f7b0acc4464e9fad7666f38f22d505": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "bc2dbb06483b4622844935307612d90b": { "model_module": "@jupyter-widgets/base", "model_module_version": "2.0.0", "model_name": "LayoutModel", "state": {} }, "d596fab2890444cbbfb875dafe62123e": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_0ad6b6964bc94e04b178fa6bdff4306e", "style": "IPY_MODEL_4847f791c69e42979bae07ea1f897c87", "value": "" } }, "eb4be0ff2bce47a3b144bc641af7dabf": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "edbcaba1de234543b4eae485fa286675": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLStyleModel", "state": { "description_width": "", "font_size": null, "text_color": null } }, "f828e2467c7147a68839ff8a8840e724": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_41e8004cc50547f2ad6b232850da0248", "style": "IPY_MODEL_29ad252eecf94b5b8cd96b2d98d10f1c", "value": "" } }, "fe427b2912844e01b761ed5be3eff1e4": { "model_module": "@jupyter-widgets/controls", "model_module_version": "2.0.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_bc2dbb06483b4622844935307612d90b", "style": "IPY_MODEL_eb4be0ff2bce47a3b144bc641af7dabf", "value": "" } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }