{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Checkout www.pygimli.org for more examples"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# The mesh class\n\nThe mesh class holds either geometric definitions (piece-wise linear complex -\nPLC) or discretizations of the subsurface. It contains of nodes,\nboundaries (edges in 2D, faces in 3D) and cells (triangles, quadrangles in 2D,\nhexahedra or tetrahedra in 3D) with associated markers and arbitrary data for\nnodes or cells.\n\nhttps://www.pygimli.org/_tutorials_auto/1_mesh/plot_2-anatomy_of_a_mesh.html\ngives a good overview on the properties of a mesh. Here we demonstrate how to\ncreate and manipulate meshes from the scratch.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We start off with the typical imports\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport pygimli as pg"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We construct a mesh by creating nodes and cells.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh = pg.Mesh(2)  # 2D\nn11 = mesh.createNode((0.0, 0.0))\nn12 = mesh.createNode((1.0, 0.0))\nn13 = mesh.createNode((2.0, 0.0))\nn21 = mesh.createNode((0.0, 1.0))\nn22 = mesh.createNode((1.0, 1.0), marker=1)\nn23 = mesh.createNode((2.0, 1.0))\nn31 = mesh.createNode((0.5, 1.7))\nn32 = mesh.createNode((1.5, 1.7))\nmesh.createQuadrangle(n11, n12, n22, n21, marker=4)\nmesh.createQuadrangle(n12, n13, n23, n22, marker=4)\nmesh.createTriangle(n21, n22, n31, marker=3)\nmesh.createTriangle(n22, n23, n32, marker=3)\nmesh.createTriangle(n31, n32, n22, marker=3)  # leave out\nprint(mesh)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function createNeighborInfos adds boundaries to nodes and cells.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh.createNeighborInfos()\nprint(mesh)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can look only at a given mesh and add matplotlib objects to the plot\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax, _ = pg.show(mesh, showMesh=True)\n\nfor n in mesh.nodes():\n    ax.text(n.x(), n.y(), str(n.id()), color=\"C0\")\nfor b in mesh.boundaries():\n    ax.text(b.center().x(), b.center().y(), str(b.id()), color=\"C1\")\nfor c in mesh.cells():\n    ax.text(c.center().x(), c.center().y(), str(c.id()), color=\"C2\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Or we can change and show all (cell or boundary) markers of a mesh\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh.boundary(2).setMarker(1)\nmesh.boundary(6).setMarker(2)\nax, _ = pg.show(mesh, markers=True, showMesh=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can iterate over all nodes, cells, or boundaries and retrieve or change\ntheir properties like node positions,\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for n in mesh.nodes():\n    print(n.id(), n.pos())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "find all nodes of a given cell or find the neighbors,\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for c in mesh.cells():\n    print(c.id(), len(c.nodes()), c.center())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "or find the neighboring cells for all inner boundaries.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for b in mesh.boundaries():\n    print(b.id(), \": Nodes:\", b.ids(), end=\" \")\n    if not b.outside():\n        left = b.leftCell()\n        right = b.rightCell()\n        print(left.id(), right.id())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We visualize some related property and attribute it the mesh cells.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "voltage = np.arange(mesh.nodeCount()) + 10\nax, cb = pg.show(mesh, voltage)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Node-based data are typically shown in form of contour lines.\n\nWe can add this (node-based) vector to the mesh as a property.\nThe same can be done for cell-based properties.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh[\"voltage\"] = voltage\nmesh[\"velocity\"] = np.arange(mesh.cellCount()) + 2\nax, cb = pg.show(mesh, \"velocity\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Cell-based data are, on the other hand, valid for the whole cell, which is\nwhy they are typically shown by filled patches.\n\nThe VTK format can save these properties along with the mesh, point data like\nthe voltage under POINT_DATA and cell data like velocity under CELL_DATA.\nIt is particularly suited to save inversion results in one file.\n3D vtk files can be nicely visualized by a number of programs, e.g. Paraview.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "if False:\n    mesh.exportVTK(\"mesh.vtk\")\n    with open(\"mesh.vtk\") as f:\n        print(f.read())"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}