{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mesh\n",
    "\n",
    "See [Mesh in MIKE IO Documentation](https://dhi.github.io/mikeio/mesh.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import mikeio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple mesh\n",
    "\n",
    "Let's consider a simple mesh consisting of 2 triangular elements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = \"data/two_elements.mesh\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(fn, \"r\") as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh = mikeio.open(fn)\n",
    "msh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot(show_mesh=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.node_coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.element_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.element_coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.get_element_area()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the node and element coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xn, yn = msh.node_coordinates[:,0], msh.node_coordinates[:,1]\n",
    "xe, ye = msh.element_coordinates[:,0], msh.element_coordinates[:,1]\n",
    "\n",
    "ax = msh.plot(show_mesh=True)\n",
    "ax.plot(xn, yn, 'ro', markersize=10)\n",
    "ax.plot(xe, ye, 'bx', markersize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundary polylines\n",
    "\n",
    "It can sometimes be convenient to have mesh boundary as a polyline (or multiple in case of more complex meshes). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bxy = msh.boundary_polylines.exteriors[0].xy\n",
    "plt.plot(bxy[:,0], bxy[:,1])\n",
    "plt.axis(\"equal\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inside domain?\n",
    "\n",
    "MIKE IO has a method for determining if a point (or a list of points) is inside the domain: \n",
    "\n",
    "* contains()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt_1 = [2.0, 1.2]\n",
    "msh.contains(pt_1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or multiple points at the same time\n",
    "pt_2 = [4.0, 1.2]\n",
    "pts = np.array([pt_1, pt_2])\n",
    "msh.contains(pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(bxy[:,0], bxy[:,1], label='boundary')\n",
    "plt.plot(xe[0], ye[0], 'b*', markersize=10, label=\"center, elem 0\")\n",
    "plt.plot(xe[1], ye[1], 'c*', markersize=10, label=\"center, elem 1\")\n",
    "plt.plot(*pt_1, 'go', markersize=10, label=\"pt_1\")\n",
    "plt.plot(*pt_2, 'rs', markersize=10, label=\"pt_2\")\n",
    "plt.axis(\"equal\")\n",
    "plt.legend(loc=\"upper right\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find element containing point\n",
    "\n",
    "MIKE IO has a method for obtaining the index of the element *containing* a point: \n",
    "\n",
    "* find_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = msh.geometry\n",
    "g.find_index(coords=pt_1)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MIKE IO also has a method for obtaining a list of the n *closest* element centers: \n",
    "\n",
    "* find_nearest_elements()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.find_nearest_elements(pt_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.find_nearest_elements(pt_1, return_distances=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.find_nearest_elements(pt_1, n_nearest=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for multiple points\n",
    "g.find_nearest_elements(pts, return_distances=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A larger mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs = mikeio.open(\"data/FakeLake.dfsu\")\n",
    "g = dfs.geometry\n",
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.plot();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inline Exercise\n",
    "\n",
    "1) please check if the point A: (x,y)=(-0.5, 0.0) is inside the mesh\n",
    "2) please check if the point B: (x,y)=(-0.5, 0.4) is inside the mesh\n",
    "3) find index of the 5 closest points to B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.max_nodes_per_element"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh = mikeio.open(\"data/FakeLake.dfsu\").geometry\n",
    "msh.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.node_coordinates[:,2] = np.clip(msh.node_coordinates[:,2], -15, 0) # clip depth to interval [-15,0]\n",
    "\n",
    "\n",
    "msh.plot(title=\"No change??\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del msh.element_coordinates # remove cached element coords calculated based on original node coords)\n",
    "msh.plot(title=\"Updated\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.to_mesh('Fake_lake_clip15.mesh')   # save to a new file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh = mikeio.open(\"data/southern_north_sea.mesh\")\n",
    "msh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default is to plot the elements and color them according to the bathymetry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot.outline();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot.mesh();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we would like to higlight the bathymetric variations in some range, in this case in the -40, -20m range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot(vmin=-40, vmax=-20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are other options as well, such as explicit specification of which contour lines to show or choosing a specific colormap ([matplotlib colormaps](https://matplotlib.org/stable/tutorials/colors/colormaps.html))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msh.plot.contour(show_mesh=True, \n",
    "         levels=[-50,-30,-20,-10,-5], cmap=\"tab10\",\n",
    "         figsize=(12,12), title=\"Coarse North Sea model\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the [MIKE IO Mesh Example notebook](https://nbviewer.jupyter.org/github/DHI/mikeio/blob/main/notebooks/Mesh.ipynb) for more Mesh operations (including shapely operations)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.5 ('base')",
   "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.10.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7aec4f91c09090e98e6ae8203c38529831bb4a3ce54cd1b69639b53cb01a6aa9"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}