{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2cf40ca0",
   "metadata": {},
   "source": [
    "# Create mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# -----------------------------------------------------------------------------\n",
    "#\n",
    "#  Gmsh Python extended tutorial 1\n",
    "#\n",
    "#  Geometry and mesh data\n",
    "#  Reference:\n",
    "# https://mail.google.com/mail/u/0/#search/mesh/KtbxLwhGPJBwdSHgPpGvkDnkxlXPmdPpdV\n",
    "# https://gitlab.onelab.info/gmsh/gmsh/blob/gmsh_4_7_1/tutorial/python/x1.py\n",
    "# https://forums.autodesk.com/t5/autocad-forum/adding-points-to-an-existing-geometry/m-p/9505285#M1019586\n",
    "# https://gitlab.onelab.info/gmsh/gmsh/-/issues/1152 - coherence commands\n",
    "# ------------------------------------------------------------------------------\n",
    "# Note to self:\n",
    "# While creating mesh from gmsh GUI using iges file exported from AutoCAD as input, make sure the units in iges file are\n",
    "# are in micrometers and not in centimeters (option 9, 2HUM and not 10, 2HCM)\n",
    "\n",
    "# STEPS:\n",
    "# manual : trial0\n",
    "# 1. create dxf file in cad\n",
    "# 2. import dxf in comsol and create mesh by specifying min/max element size\n",
    "# 3. export the mesh points from Mesh node in comsol to mphtxt\n",
    "# 4. copy and save mesh points from mphtxt to txt file\n",
    "# 5. In autocad, add points from step 4 to the geometry via Points option in autocad; save dxf\n",
    "# 6. import dxf from step 5 into comsol and export the geometry from Geometry node in .step format (issue with iges)\n",
    "# 7. import step in gmsh's GUI, remove duplicate nodes via Coherence -> save .geo; click mesh -> 1D to create elements\n",
    "# 8. save .msh format\n",
    "# 9. read .msh contents via meshio, generate graph\n",
    "# 10. Note: directly exporting iges from cad via IGESEXPORT command doesn't help; not possible to generate elements\n",
    "\n",
    "\n",
    "# Automated:\n",
    "# 1. run extract_subgraph to retain just one inlet and one outlet in the user-defined graph\n",
    "# 2. update the network topology in input.xlsx\n",
    "# 3. run the untwiched version of the code to find out the right direction of flow (velocity computation) in MATLAB\n",
    "# 4. update the edges with the right direction of flow in input.xlsx\n",
    "# 5. run read_mesh.py to generate the coordinates of the mesh elements\n",
    "# 6. update the pos_mesh coordinates in input.xlsx\n",
    "# 7. run concentration simulations\n",
    "\n",
    "# Note to self: Interpreter : python 3.6\n",
    "# ------------------------------------------------------------------------------\n",
    "\n",
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "import gmsh\n",
    "import pygmsh\n",
    "import meshio\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from vedo import *\n",
    "# from vedo import io\n",
    "from collections import OrderedDict\n",
    "from matpancreas.settings_model import GRAPH_MESH_FILE\n",
    "from matpancreas.utils.graph_utils_py import draw_graph3d_vedo, convert_graph_to_df, get_eucledian_lengths, \\\n",
    "    convert_coordinate_units, create_graph, get_eucledian_lengths_wrt_origin, remove_edge_attribute, get_graph\n",
    "from matpancreas.utils.io_utils_py import write_output\n",
    "\n",
    "\n",
    "def mesh_network(G):\n",
    "\n",
    "    \"\"\"\n",
    "    loads a networkx graph and meshed the edges (/lines) to add intermediate nodes via gmsh library\n",
    "    References:\n",
    "    ----------\n",
    "    https://stackoverflow.com/a/431868/8281509\n",
    "    \"\"\"\n",
    "    # draw_graph3d_vedo([G], label_scale=5) #G=G, point_r=1, lw=1)\n",
    "    H = G.copy()\n",
    "    H = remove_edge_attribute(H, edge_attr=['l'])\n",
    "    pos = nx.get_node_attributes(G, 'pos')\n",
    "\n",
    "    # generate mesh\n",
    "    gmsh.initialize()\n",
    "    p = [gmsh.model.geo.addPoint(*point, tag=tag) for tag, point in pos.items()]\n",
    "    gmsh.model.geo.synchronize()\n",
    "    # print(gmsh.model.getEntities())\n",
    "\n",
    "    for i, e in enumerate(sorted(G.edges())):\n",
    "        t, h = e\n",
    "        gmsh.model.geo.addLine(t, h, tag=i)\n",
    "    gmsh.model.geo.synchronize()\n",
    "\n",
    "    gmsh.option.setNumber(\"Mesh.MeshSizeMin\", 9.5)\n",
    "    gmsh.option.setNumber(\"Mesh.MeshSizeMax\", 13.5)\n",
    "    gmsh.model.mesh.generate(1)  # make a 1d mesh\n",
    "    gmsh.write('test.msh')\n",
    "\n",
    "    for i, e in enumerate(sorted(G.edges(data=True))):\n",
    "        t, h, attr = e\n",
    "        node_tags, node_coords = gmsh.model.mesh.getNodes(dim=1, tag=i)[0:2]\n",
    "        # node_coords = np.reshape(node_coords, len(node_coords))\n",
    "        nx.add_path(H, nodes_for_path=[t, *node_tags, h], d=attr['d'])\n",
    "        nx.set_edge_attributes(G, {(t, h): len(node_tags)}, 'segment_i_node')\n",
    "        # initial edge is removed when meshing creates interior nodes\n",
    "        if len(node_tags) > 0:\n",
    "            H.remove_edge(t, h)\n",
    "    # ------------------------------------------------------------------------------------------------------------------\n",
    "    # set positions of newly added mesh nodes and edge lengths\n",
    "    # ------------------------------------------------------------------------------------------------------------------\n",
    "    mesh_nodes, mesh_coords = gmsh.model.mesh.getNodes(dim=1)[0:2]\n",
    "    mesh_coords = np.reshape(mesh_coords, (len(mesh_nodes), 3))\n",
    "    mesh_pos = dict(zip(mesh_nodes, mesh_coords))\n",
    "\n",
    "    nx.set_node_attributes(H, mesh_pos, 'pos')\n",
    "    H = get_eucledian_lengths(H)\n",
    "    gmsh.finalize()\n",
    "\n",
    "    return G, H\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "\n",
    "    G = create_graph(attributes=True, actual_pos=False, directed=True)\n",
    "    exit()\n",
    "    # draw_graph3d_vedo([G], label_scale=3)\n",
    "    # G.remove_edges_from(\n",
    "    #     ebunch=[(289, 303), (291, 337), (339, 329),\n",
    "    #             (14, 12), (229, 308), (168, 171),\n",
    "    #             (323, 251), (144, 145), (68, 83),\n",
    "    #             (32, 30), (172, 174), (314, 338),\n",
    "    #             (178, 180), (88, 86), (159, 11)]\n",
    "    # )\n",
    "\n",
    "    G = get_eucledian_lengths(G)\n",
    "    G, H = mesh_network(G=G)\n",
    "    graph_df = convert_graph_to_df(H)\n",
    "    write_output(d={'islet': graph_df}, file='data_model')\n",
    "    # draw_graph3d_vedo([H], label_scale=3) #G=G, point_r=1, lw=1)"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "67be7881b41419f4"
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}