{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2daad93c-8820-47dd-8419-2640531b4c78",
   "metadata": {},
   "source": [
    "# Examples of usage of filters\n",
    "\n",
    "\n",
    "Here we show how to use node and element filters. The main Muscat modules used in this example are listed below:\n",
    "\n",
    "- `ElementFilter`: allows for filtering the elements within your mesh.\n",
    "- `NodeFilter`: allows for filtering the nodes within your mesh.\n",
    "- `IntersectionFilter`: useful for combining intersecting filters together.\n",
    "- `UnionFilter`: useful for combining filters together.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72888afe-6d57-4f58-a811-956072ab996b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.Containers.Filters.FilterTools import ElementFilter, NodeFilter\n",
    "from Muscat.Containers.Filters.FilterOperators import UnionFilter, IntersectionFilter\n",
    "import numpy as np\n",
    "from scipy.spatial import Delaunay\n",
    "from Muscat.Containers import MeshCreationTools as UMCT\n",
    "import pyvista\n",
    "pyvista.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots\n",
    "\n",
    "def PlotnTag(mesh, tag):\n",
    "    plotter = pv.Plotter()\n",
    "    actor = plotter.add_mesh(MeshToPyVista(mesh,TagsAsFields=True),show_edges=True, edge_color=\"grey\", style=\"wireframe\")\n",
    "    actor = plotter.add_mesh(MeshToPyVista(mesh,TagsAsFields=True), scalars=tag,show_edges=True, edge_color=\"grey\", style=\"points\")\n",
    "    plotter.enable_parallel_projection()\n",
    "    plotter.show(cpos=\"xy\")\n",
    "    \n",
    "def PloteTag(mesh, tag):\n",
    "    plotter = pv.Plotter()\n",
    "    actor = plotter.add_mesh(MeshToPyVista(mesh,TagsAsFields=True),scalars=tag, show_edges=True, edge_color=\"grey\", style=\"surface\")\n",
    "    #actor = plotter.add_mesh(MeshToPyVista(mesh,TagsAsFields=True), scalars=tag,show_edges=True, edge_color=\"grey\", style=\"points\")\n",
    "    plotter.enable_parallel_projection()\n",
    "    plotter.show(cpos=\"xy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7c2028d-c16b-4800-8d87-46c8743d208a",
   "metadata": {},
   "source": [
    "For this simple example, we generate a dummy mesh using `scipy.spatial.Delaunay`,\n",
    "and create a `Mesh` object using `Muscat.Containers.MeshCreationTools`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d25a4187-35a5-4775-a03d-501b8a9f4d8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mesh_example():\n",
    "    \"\"\"This function creates a dummy mesh using Delaunay triangulation.\n",
    "    \"\"\"\n",
    "    # Create a grid of points\n",
    "    x, y = np.meshgrid(np.linspace(0,1,100), np.linspace(0,1,100))\n",
    "    points = np.stack([x.ravel(), y.ravel()], axis=1)\n",
    "    # Generate the triangles with Delaunay\n",
    "    tri = Delaunay(points)\n",
    "    triangles = tri.simplices\n",
    "    # Create a Muscat Mesh using the CreateMeshOfTriangles utility\n",
    "    mesh = UMCT.CreateMeshOfTriangles(points, triangles)\n",
    "    return mesh\n",
    "\n",
    "mesh = create_mesh_example()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef00cb27-db0e-4b2e-9090-ec6450fbf506",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.Bridges.PyVistaBridge import MeshToPyVista\n",
    "import pyvista as pv\n",
    "\n",
    "PloteTag(mesh,\"2D\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ee71ad8-6017-41b3-b762-a9728f2aa73b",
   "metadata": {},
   "source": [
    "Node filters can be used to define sets of nodes satisfying an inequality\n",
    "constraint of the form :math:`f(x) < 0`. \n",
    "\n",
    "In the example below, we select nodes that lie within the \n",
    "ball centered at :math:`(0.5, 0.5)` and with radius :math:`0.25`.\n",
    "\n",
    "The node filter can then be used to get the indices of the nodes satisfying the constraint,\n",
    "which in turn can be used to create a new nodal tag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37bed8a2-8554-465a-bfae-3b5e93768e09",
   "metadata": {},
   "outputs": [],
   "source": [
    "nf = NodeFilter(zone=lambda p: (np.linalg.norm(p - 0.5, axis=1) - 0.25))\n",
    "indices = nf.GetNodesIndices(mesh)\n",
    "mesh.GetNodalTag(\"circle\").AddToTag(indices)\n",
    "\n",
    "PlotnTag(mesh, \"circle\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88638acc-cab2-40d1-a5ac-0f735f782ca4",
   "metadata": {},
   "source": [
    "Multiple filters can be combined together using `UnionFilter`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "130ebfe7-fd39-44c6-8f0d-147592939c86",
   "metadata": {},
   "outputs": [],
   "source": [
    "nf1 = NodeFilter(zone=lambda p: (np.linalg.norm(p - 0.6, axis=1) - 0.25))\n",
    "nf2 = NodeFilter(zone=lambda p: (np.linalg.norm(p - 0.4, axis=1) - 0.25))\n",
    "nf = UnionFilter(filters=[nf1, nf2])\n",
    "\n",
    "indices = nf.GetNodesIndices(mesh)\n",
    "mesh.GetNodalTag(\"two_circles_union\").AddToTag(indices)\n",
    "\n",
    "PlotnTag(mesh, \"two_circles_union\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6992625e-5299-4d48-b9e9-f3b8ee9aa2be",
   "metadata": {},
   "source": [
    "We can also compute the intersection of multiple filters using `IntersectionFilter`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73a9cc91-d46d-4781-b31a-80ca697ea9fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "nf1 = NodeFilter(zone=lambda p: (np.linalg.norm(p - 0.6, axis=1) - 0.25))\n",
    "nf2 = NodeFilter(zone=lambda p: (np.linalg.norm(p - 0.4, axis=1) - 0.25))\n",
    "nf = IntersectionFilter(filters=[nf1, nf2])\n",
    "\n",
    "indices = nf.GetNodesIndices(mesh)\n",
    "mesh.GetNodalTag(\"two_circles_intersection\").AddToTag(indices)\n",
    "PlotnTag(mesh, \"two_circles_intersection\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa44513-50eb-45c1-b8c2-46a88b76ba65",
   "metadata": {},
   "source": [
    "The same logic applies to `ElementFilter` in Muscat. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a509ba7-854e-4d1e-8213-ec2af64898bf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ne = ElementFilter(zone=lambda p: (np.linalg.norm(p - 0.5, axis=1) - 0.25))\n",
    "for selection in ne(mesh):\n",
    "    selection.elements.tags.CreateTag(\"circle_element_filter\").SetIds(selection.indices)\n",
    "#indices = ne.GetNodesIndices(mesh)\n",
    "#mesh.GetNodalTag(\"circle_element_filter\").AddToTag(indices)\n",
    "PloteTag(mesh, \"circle_element_filter\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fe877c9-0e30-47d6-92f5-d476a17c037f",
   "metadata": {},
   "source": [
    "Element filters can also be combined together in a similar fashion.\n",
    "\n",
    "Elements and nodes filters can be used together to define new sets of nodes. \n",
    "As an illustration, applying the `ComputeSkin` method to our mesh \n",
    "creates a new mesh with an additional element tag `Skin`. The latter\n",
    "corresponds to the set of elements that lie on the boundary of the \n",
    "mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a71ce3d-2eb1-464f-bcea-7c4b656b6fd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.Containers import MeshModificationTools as MMT\n",
    "\n",
    "new_mesh = MMT.ComputeSkin(mesh, inPlace=True)\n",
    "#Printing the mesh shows all the tags that we added previsouly and the new `Skin` tag for elements.\n",
    "print(new_mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d2fc21e-9f7d-4fde-b648-9ddba367b4c7",
   "metadata": {},
   "source": [
    "We can now easily determine the nodes that lie on the top of the mesh\n",
    "using a `NodeFilter` that only considers the elements belonging to `Skin`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee6275f3-2dba-47b9-a7f7-cc0bcbe4c483",
   "metadata": {},
   "outputs": [],
   "source": [
    "nf = NodeFilter(eTag=\"Skin\", zone = lambda p: -p[:, 1] + 0.99)\n",
    "indices = nf.GetNodesIndices(new_mesh)\n",
    "new_mesh.GetNodalTag(\"top\").AddToTag(indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4991a1b6-8f8a-4e1d-be06-d2e9a03c7d00",
   "metadata": {},
   "source": [
    "The mesh now has a new nodal tag `top`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be5c1043-0ef4-4c62-8376-3093452b51dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(new_mesh)\n",
    "PlotnTag(new_mesh,'top')"
   ]
  }
 ],
 "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
}