{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "63228282-76b3-420d-90a4-d4ed2852202e",
   "metadata": {},
   "source": [
    "# Simulation éléments finis d'un poutre en flexion 3 points\n",
    "\n",
    "Le code éléments finis est `wombat`, un code FEM pédagogique développé par [Jérémy Bleyer (ENPC)](https://bleyerj.github.io/) pour le cours d'éléments finis de l'École des Ponts et Chaussées.\n",
    "\n",
    "Les quelques cellules ci-dessous contiennent des fonctions servant à mailler la poutre et post-traiter les résultats du calcul par éléments finis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49284038-e011-4739-8b28-c0826329c1ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\", category=SyntaxWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b9eb781-01ed-4c48-a1c5-0bcd3b86db00",
   "metadata": {},
   "outputs": [],
   "source": [
    "from wombat import *\n",
    "\n",
    "import numpy as np\n",
    "import pyvista as pv\n",
    "\n",
    "from scipy.spatial import Delaunay\n",
    "\n",
    "pv.set_jupyter_backend(\"static\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46b15412-e6f1-4032-896f-25c67bd7ea28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def maillage(L, h, taille):\n",
    "    \"\"\"Crée un maillage 2D de poutre avec une taille des éléments\"\"\"\n",
    "    assert (taille < L) and (taille < h)\n",
    "\n",
    "    nx = int(L / taille) + 1\n",
    "    ny = int(h / taille) + 1\n",
    "\n",
    "    x = np.linspace(0, L / 2, nx)\n",
    "    y = np.linspace(0, h, ny)\n",
    "    xx, yy = np.meshgrid(x, y, indexing='ij')\n",
    "\n",
    "    X = np.array([xx.ravel(), yy.ravel()]).T\n",
    "\n",
    "    tri = Delaunay(X)\n",
    "\n",
    "    nodes = np.array([Node2D(x) for x in X])\n",
    "    elem_list = [SolidT3(nodes[el]) for el in tri.simplices]\n",
    "    mesh = Mesh(elem_list, points=X)\n",
    "\n",
    "    def make_boundary_group(label, predicate, dim):\n",
    "        group_nodes = np.array([node for node in nodes if predicate(node)])\n",
    "        node_coords = np.array([node.coor for node in group_nodes])\n",
    "        sort = np.argsort(node_coords[:, dim])\n",
    "        group_nodes = group_nodes[sort]\n",
    "        N = group_nodes.shape[0]\n",
    "\n",
    "        return [TraceSolidT3([group_nodes[i], group_nodes[i + 1]], tag=label) for i in range(N-1)]\n",
    "\n",
    "    boundary = ElementGroup(\n",
    "        make_boundary_group(\"bottom\", lambda n: n.coor[1] == 0., dim=0)\n",
    "        + make_boundary_group(\"left\", lambda n: n.coor[0] == 0., dim=1)\n",
    "        + make_boundary_group(\"right\", lambda n: n.coor[0] == L, dim=1)\n",
    "        + make_boundary_group(\"top\", lambda n: n.coor[1] == h, dim=0)\n",
    "    )       \n",
    "    \n",
    "    return mesh, boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f9fd475-2683-4265-8191-62e8d14cbc92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_mesh(mesh, appuis=None):\n",
    "    fig = Figure(1, \"Maillage\")\n",
    "    fig.plot(mesh)\n",
    "    if appuis:\n",
    "        fig.plot_bc(mesh, appuis)\n",
    "\n",
    "def plot_deformee(mesh, U, facteur_echelle=1.):\n",
    "    pvmesh = pv.UnstructuredGrid({pv.CellType.TRIANGLE: mesh.connec}, np.hstack([mesh.coor, np.zeros((mesh.coor.shape[0], 1))]))\n",
    "    pvmesh.point_data['deplacement'] = np.array([U.get(\"ux\"), U.get(\"uy\"), np.zeros_like(U.get('ux'))]).T * facteur_echelle\n",
    "    warp = pvmesh.warp_by_vector()\n",
    "    pl = pv.Plotter()\n",
    "    pl.add_mesh(pvmesh.outline(), color='black')\n",
    "    pl.add_mesh(warp, color='white', opacity=0.5)\n",
    "    pl.show_bounds(padding=0.1)\n",
    "    pl.camera_position = 'xy'\n",
    "    pl.show()\n",
    "\n",
    "def plot_contrainte(mesh, sigma, nom_contrainte, contours=5, cmap=\"viridis\"):\n",
    "    pvmesh = pv.UnstructuredGrid({pv.CellType.TRIANGLE: mesh.connec}, np.hstack([mesh.coor, np.zeros((mesh.coor.shape[0], 1))]))\n",
    "\n",
    "    pvmesh.cell_data[nom_contrainte] = sigma\n",
    "    pvmesh = pvmesh.cell_data_to_point_data()\n",
    "\n",
    "    pl = pv.Plotter()\n",
    "    pl.add_mesh(pvmesh, n_colors=contours, cmap=cmap)\n",
    "    pl.camera_position = 'xy'\n",
    "    pl.show()\n",
    "    return pvmesh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7e851cf-9768-4c31-b2d3-c7b9d49d0e97",
   "metadata": {},
   "source": [
    "## Géométrie\n",
    "\n",
    "On définie ici la géométrie de la poutre. **Attention aux unités choisies** !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5d575a5-8c5b-4f18-a578-79ec20ed2101",
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################\n",
    "# À remplir\n",
    "#############################\n",
    "\n",
    "# Géométrie\n",
    "L = ... # longueur totale de la poutre\n",
    "h = ... # hauteur de la poutre\n",
    "\n",
    "# Taille de maillage\n",
    "taille = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bae13ceb-c95c-4d90-b7bd-d4770f80bf88",
   "metadata": {},
   "source": [
    "## Création du maillage et affichage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1a893df-218d-44a0-9a18-84c3cbe30cad",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, boundary = maillage(L, h, taille)\n",
    "\n",
    "bottom = boundary.get_elem_from_tag(\"bottom\")\n",
    "right = boundary.get_elem_from_tag(\"right\")\n",
    "top = boundary.get_elem_from_tag(\"top\")\n",
    "left = boundary.get_elem_from_tag(\"left\")\n",
    "top_left = left.get_nodes()[-1]\n",
    "bottom_right = right.get_nodes()[0]\n",
    "\n",
    "plot_mesh(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70c94518-a3ab-4551-8acb-5dd8f5ea8667",
   "metadata": {},
   "source": [
    "## Caractéristiques du matériau\n",
    "\n",
    "On défini le module d'Young et le coefficient de Poisson du matériau. **Attention aux unités** !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e217aa5-5d3a-4a4c-b72c-06b9644118a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = LinearElastic(E=..., nu=...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73554d76-a843-4a0c-b73e-a2e6c6fc7d63",
   "metadata": {},
   "source": [
    "## Conditions limites\n",
    "\n",
    "Il faut utiliser `appuis.add_imposed_displ(group, ux=..., uy=...)` pour imposer les déplacements bloqués sur les bons groupes de nœuds.\n",
    "\n",
    "Il faut utiliser `forces.add_concentrated_forces(group, Fx=..., Fy=...)` pour imposer une force ponctuelle sur un nœud. **Attention aux unités** !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36e39a39-e592-478b-bd0c-6c6c238c575e",
   "metadata": {},
   "outputs": [],
   "source": [
    "appuis = Connections()\n",
    "appuis.add_imposed_displ(left, ux=0)\n",
    "appuis.add_imposed_displ(bottom_right, ux=0, uy=0)\n",
    "\n",
    "forces = ExtForce()\n",
    "forces.add_concentrated_forces(top_left, Fy=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad4b43ac-278e-455d-9776-5cad0dc04460",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(mesh,mat)\n",
    "\n",
    "plot_mesh(mesh, appuis)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fd5751f-62e8-4bf3-b54e-584c3339f597",
   "metadata": {},
   "source": [
    "## Résolution du système éléments finis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2faea11-87a6-4da4-8ea1-8ffb501ae4a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "K = model.assembl_stiffness_matrix()\n",
    "Lmat, Ud = model.assembl_connections(appuis)\n",
    "F = model.assembl_external_forces(forces)\n",
    "\n",
    "U, lamb = model.solve(K, F, Lmat, Ud)\n",
    "Sigma = model.stresses(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb7bdd09-0833-48d3-b50d-c22ca59ac98b",
   "metadata": {},
   "source": [
    "## Post-traitement des résultats\n",
    "\n",
    "- `U` est un objet contenant les déplacements de tous les nœuds du maillage. `U.get(\"ux\")` retourne le déplacement horizontal, `U.get(\"uy\")` le déplacement vertical.\n",
    "- `Sigma` est un objet contenant les composantes du tenseur des contraintes. `Sigma.get(\"sig_xx\")` retourne $\\sigma_{11}$.\n",
    "\n",
    "On souhaite afficher dans l'ordre:\n",
    "\n",
    "1. La déformée, avec `plot_deformee`\n",
    "2. La contrainte longitudinale $\\sigma_{11}$ avec `plot_contrainte`\n",
    "3. La contrainte longitudinale $\\sigma_{11}$ le long de la droite $x_1 = L/2$, comparer avec la solution analytique\n",
    "4. La contrainte de cisaillement $\\sigma_{12}$ sur toute la poutre\n",
    "5. La contrainte de cisaillement $\\sigma_{12}$ le long de la droite $x_1 = L/2$, comparer avec la solution analytique\n",
    "\n",
    "Pour faire les graphes le long d'une droite, on pourra utiliser le code suivant:\n",
    "\n",
    "```python\n",
    "pvmesh = plot_contrainte(...)\n",
    "\n",
    "# Points de départ et d'arriver du segment sur lequel tracer la courbe\n",
    "xa, ya = ..., ...\n",
    "xb, yb = ..., ...\n",
    "\n",
    "pvmesh.plot_over_line([xa, ya, 0], [xb, yb, 0])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68961d4d-0772-4258-8a56-18dc391dfcbf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}