examples/example10/DFN_meshing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Meshing of Discrete Fracture Networks\n",
    "\n",
    "This tutorial shows how to create a GridBucket that represents a discrete fracture matrix model. Two meshing strategies are available; they differ in how they treat interseciton lines: \n",
    "1. A fully conforming mesh, where the fracture grid that meet along the intersection have common nodes. This is the simplest option for numerical methods, with no hanging nodes. However, it requires joint meshing of all fractures in the network.\n",
    "2. A partially conforming mesh, where the intersection lines are resolved in all fracture grids, but the grid nodes are not necessarily matching. The grid will thus have hanging nodes, and should be treated by numerical methods that can deal with this. In PorePy, try the virtual element method, and probably TPFA. On the positive side, the fracture planes can be meshed independently, which can be a great benefit for complex geometries.\n",
    "\n",
    "Importantly, both strategies result in a GridBucket which can be used in further numerical computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports needed\n",
    "import numpy as np\n",
    "import porepy as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create three fractures with a common intersection point at the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_1 = pp.Fracture(np.array([[-2, -1, 0], [1, -1, 0], [1, 1, 0], [-2, 1, 0]]).T)\n",
    "f_2 = pp.Fracture(np.array([[-1, 0, -1], [1, 0, -1], [1, 0, 1], [-1, 0, 1]]).T)\n",
    "f_3 = pp.Fracture(np.array([[0, -1, -1], [0, 1, -1], [0, 1, 1], [0, -1, 1]]).T)\n",
    "\n",
    "network = pp.FractureNetwork([f_1, f_2, f_3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fully conforming mesh\n",
    "\n",
    "First create a fully conforming network. This will mesh the fracture networks simultaneously, and ensure that the nodes on the intersections are coinciding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found no information on mesh sizes. Returning\n",
      "\n",
      "\n",
      "Created 3 2-d grids with 359 cells\n",
      "Created 6 1-d grids with 18 cells\n",
      "Created 1 0-d grids with 1 cells\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gb_conforming = pp.meshing.dfn(network, conforming=True, verbose=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partially conforming mesh\n",
    "Next, create a non-conforming grid. In this approach, every fracture plane is meshed independently, and the grids are then merged together with the use of hanging nodes. The merging is somewhat complex, but from the user side, it is simple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mixed dimensional grid. \n",
      "Maximum dimension 2\n",
      "Minimum dimension 1\n",
      "Size of highest dimensional grid: Cells: 478. Nodes: 381\n",
      "In lower dimensions: \n",
      "3 grids of dimension 1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gb_non_conforming = pp.meshing.dfn([f_1, f_2, f_3], conforming=False)\n",
    "# Also print out some information on the mesh\n",
    "print(gb_non_conforming)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the number of cells in 2d is different (actually higher, controlled by default mesh size parameters). \n",
    "\n",
    "Also, observe that the non-conforming mesh has no 0d objects. Including this should not be too much work, but it has not been required yet, thus not prioritized."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualized meshes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First export\n",
    "gb_conforming.assign_node_ordering()\n",
    "gb_non_conforming.assign_node_ordering()\n",
    "export_conforming = pp.Exporter(gb_conforming, 'dfn_conrforming')\n",
    "export_conforming.write_vtk()\n",
    "export_non_conforming = pp.Exporter(gb_non_conforming, 'dfn_non_conforming')\n",
    "export_non_conforming.write_vtk()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some manipulation in paraview then produce the following plots\n",
    "\n",
    "<img src=\"fig/dfn_conforming.png\" alt=\"Drawing\" style=\"width: 400px;\"/>\n",
    "Conforming mesh\n",
    "\n",
    "<img src=\"fig/dfn_non_conforming.png\" alt=\"Drawing\" style=\"width: 400px;\"/>\n",
    "Non-conforming mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is there a fully non-conforming approach?\n",
    "The fully non-conforming approach, in the terminology adapted here, would entail meshes that do not conform to the fracture intersections at all. To create such a mesh is very simple, it is a plane discretization of a 2d surface with no constraints. However, the burden of coupling the equation in each fracture plane is then fully transferred to the numerical method. No such approach is implemented in PorePy - to see how to deal with this, see for instance the work of Berrone, Scialo and coworkers in Turin."
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}