{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Import the required packages and open model files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import the required libraries\n",
    "import os, re\n",
    "import numpy as np\n",
    "\n",
    "from workingFunctions import Functions                      #functions from the workingFunctions.py file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#change directory to the script path\n",
    "os.chdir('C:/Users\\Saul\\Documents\\Ih_PlayaroundwithVTK')    #use your own path\n",
    "\n",
    "#open the DIS, BAS and FHD and DRN files\n",
    "disLines= open('Model/Modelo3.dis').readlines()             #discretization data\n",
    "basLines= open('Model/Modelo3.bas').readlines()             #active / inactive data\n",
    "fhdLines= open('Model/Modelo3.fhd').readlines()             #water head data\n",
    "drnLines= open('Model/Modelo3.drn').readlines()             #drain package data\n",
    "\n",
    "#create a empty dictionay to store the model features\n",
    "modDis = {}\n",
    "modBas = {}\n",
    "modFhd = {}\n",
    "modDrn = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with the DIS (Discretization Data) data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General model features as modDis dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#get the extreme coordinates form the dis header\n",
    "for line in disLines[:5]:\n",
    "    line =  re.sub('\\(|\\)|,',' ',line)                      #change the (, )  and , character in the dis file\n",
    "    if \"Lower left corner\" in line:\n",
    "        line = line.split()\n",
    "        modDis[\"vertexXmin\"]=float(line[4])\n",
    "        modDis[\"vertexYmin\"]=float(line[5])\n",
    "    elif \"Upper right corner\" in line:\n",
    "        line = line.split()\n",
    "        modDis[\"vertexXmax\"]=float(line[4])\n",
    "        modDis[\"vertexYmax\"]=float(line[5])\n",
    "\n",
    "#get the number of layers, rows, columns, cell and vertex numbers\n",
    "linelaycolrow = disLines[6].split()\n",
    "modDis[\"cellLays\"] = int(linelaycolrow[0])\n",
    "modDis[\"cellRows\"] = int(linelaycolrow[1])\n",
    "modDis[\"cellCols\"] = int(linelaycolrow[2])\n",
    "modDis[\"vertexLays\"] = modDis[\"cellLays\"] + 1\n",
    "modDis[\"vertexRows\"] = modDis[\"cellRows\"] + 1\n",
    "modDis[\"vertexCols\"] = modDis[\"cellCols\"] + 1\n",
    "modDis[\"vertexperlay\"] = modDis[\"vertexRows\"] * modDis[\"vertexCols\"]\n",
    "modDis[\"cellsperlay\"] = modDis[\"cellRows\"] * modDis[\"cellCols\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Get the DIS Breakers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#get the grid breakers\n",
    "modDis['disBreakers']={}\n",
    "breakerValues = [\"INTERNAL\",\"CONSTANT\"]\n",
    "\n",
    "vertexLay=0\n",
    "for item in breakerValues:\n",
    "    for line in disLines:\n",
    "        if item in line:\n",
    "            if 'DELR' in line:                                     #DELR is cell width along rows\n",
    "                modDis['disBreakers']['DELR']=disLines.index(line) \n",
    "            elif 'DELC' in line:                                   #DELC is cell width along columns\n",
    "                modDis['disBreakers']['DELC']=disLines.index(line)\n",
    "            else: \n",
    "                modDis['disBreakers']['vertexLay'+str(vertexLay)]=disLines.index(line)\n",
    "                vertexLay+=1\n",
    "                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the DEL Info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modDis['DELR'] = Functions.getListFromDEL(modDis['disBreakers']['DELR'],disLines,modDis['cellCols'])\n",
    "modDis['DELC'] = Functions.getListFromDEL(modDis['disBreakers']['DELC'],disLines,modDis['cellRows'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the Cell Centroid Z "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modDis['cellCentroidZList']={}\n",
    "\n",
    "for lay in range(modDis['vertexLays']):\n",
    "    \n",
    "    #add auxiliar variables to identify breakers\n",
    "    lineaBreaker = modDis['disBreakers']['vertexLay'+str(lay)]\n",
    "     \n",
    "    #two cases in breaker line\n",
    "    if 'INTERNAL' in disLines[lineaBreaker]:\n",
    "        lista = Functions.getListFromBreaker(lineaBreaker,modDis,disLines)\n",
    "        modDis['cellCentroidZList']['lay'+str(lay)]=lista        \n",
    "    elif 'CONSTANT' in disLines[lineaBreaker]:\n",
    "        constElevation = float(disLines[lineaBreaker].split()[1])\n",
    "        modDis['cellCentroidZList']['lay'+str(lay)]= [constElevation for x in range(modDis[\"cellsperlay\"])]\n",
    "    else:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List of arrays of cells and vertex coord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modDis['vertexEasting']  = np.array([modDis['vertexXmin']+np.sum(modDis['DELR'][:col]) for col in range(modDis['vertexCols'])])\n",
    "modDis['vertexNorthing'] = np.array([modDis['vertexYmax']-np.sum(modDis['DELC'][:row]) for row in range(modDis['vertexRows'])])\n",
    "\n",
    "modDis['cellEasting']    = np.array([modDis['vertexXmin']+np.sum(modDis['DELR'][:col])+modDis['DELR'][col]/2 for col in range(modDis['cellCols'])])\n",
    "modDis['cellNorthing']   = np.array([modDis['vertexYmax']-np.sum(modDis['DELC'][:row])-modDis['DELC'][row]/2 for row in range(modDis['cellRows'])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation from Z cell centroid to z vertex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the BAS Info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the grid breakers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty dict to store BAS breakers\n",
    "modBas['basBreakers']={}\n",
    "\n",
    "breakerValues = [\"INTERNAL\",\"CONSTANT\"]\n",
    "\n",
    "#store the breakers in the dict\n",
    "lay=0\n",
    "for item in breakerValues:\n",
    "    for line in basLines:\n",
    "        if item in line:\n",
    "            if 'IBOUND' in line: \n",
    "                modBas['basBreakers']['lay'+str(lay)]=basLines.index(line)\n",
    "                lay+=1\n",
    "            else: \n",
    "                pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Store ibound per lay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty dict to store cell ibound per layer\n",
    "modBas['cellIboundList']={}\n",
    "\n",
    "for lay in range(modDis['cellLays']):\n",
    "    \n",
    "    #add auxiliar variables to identify breakers\n",
    "    lineaBreaker = modBas['basBreakers']['lay'+str(lay)]\n",
    "     \n",
    "    #two cases in breaker line\n",
    "    if 'INTERNAL' in basLines[lineaBreaker]:\n",
    "        lista = Functions.getListFromBreaker(lineaBreaker,modDis,basLines)\n",
    "        modBas['cellIboundList']['lay'+str(lay)]=lista        \n",
    "    elif 'CONSTANT' in basLines[lineaBreaker]:\n",
    "        constElevation = float(disLines[lineaBreaker].split()[1]) #todavia no he probado esto\n",
    "        modBas['cellIboundList']['lay'+str(lay)]= [constElevation for x in range(modDis[\"cellsperlay\"])]\n",
    "    else:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Store Cell Centroids as a Numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty list to store cell centroid\n",
    "cellCentroidList = []\n",
    "\n",
    "#numpy array of cell centroid\n",
    "for row in range(modDis['cellRows']):\n",
    "    for col in range(modDis['cellCols']):\n",
    "        cellCentroidList.append([modDis['cellEasting'][col],modDis['cellNorthing'][row]])\n",
    "\n",
    "#store cell centroids as numpy array\n",
    "modDis['cellCentroids']=np.asarray(cellCentroidList)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid of XYZ Vertex Coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modDis['vertexXgrid'] = np.repeat(modDis['vertexEasting'].reshape(modDis['vertexCols'],1),modDis['vertexRows'],axis=1).T\n",
    "modDis['vertexYgrid'] = np.repeat(modDis['vertexNorthing'],modDis['vertexCols']).reshape(modDis['vertexRows'],modDis['vertexCols'])\n",
    "modDis['vertexZGrid'] = Functions.interpolateCelltoVertex(modDis,'cellCentroidZList')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the HDS Info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the grid breakers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty dict to store HDS breakers\n",
    "modFhd['fhdBreakers']={}\n",
    "\n",
    "#store the breakers in the dict\n",
    "lay=0\n",
    "for line in fhdLines:\n",
    "    if 'HEAD' in line:\n",
    "        modFhd['fhdBreakers']['lay'+str(lay)]=fhdLines.index(line)\n",
    "        lay+=1\n",
    "    else: \n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Store heads per lay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty dict to store cell heads per layer\n",
    "modFhd['cellHeadGrid']={}\n",
    "\n",
    "for lay in range(modDis['cellLays']):\n",
    "    #add auxiliar variables to identify breakers\n",
    "    lineaBreaker = modFhd['fhdBreakers']['lay'+str(lay)]     \n",
    "    lista = Functions.getListFromBreaker(lineaBreaker,modDis,fhdLines)\n",
    "    array = np.asarray(lista).reshape(modDis['cellRows'],modDis['cellCols'])\n",
    "    modFhd['cellHeadGrid']['lay'+str(lay)]=array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrange to transform cell centroid heads to vertex heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty temporal dictionary to store transformed heads\n",
    "vertexHeadGridCentroid = {}\n",
    "\n",
    "#arrange to hace positive heads in all vertex of an active cell\n",
    "for lay in range(modDis['cellLays']):\n",
    "    matrix = np.zeros([modDis['vertexRows'],modDis['vertexCols']])\n",
    "    for row in range(modDis['cellRows']):\n",
    "        for col in range(modDis['cellCols']):\n",
    "            headLay = modFhd['cellHeadGrid']['lay'+str(lay)]\n",
    "            neighcartesianlist = [headLay[row-1,col-1],headLay[row-1,col],headLay[row,col-1],headLay[row,col]]\n",
    "            headList = []\n",
    "            for item in neighcartesianlist:\n",
    "                if item > -1e+20:\n",
    "                    headList.append(item)\n",
    "            if len(headList) > 0:\n",
    "                headMean = sum(headList)/len(headList)\n",
    "            else:\n",
    "                headMean = -1e+20\n",
    "            \n",
    "            matrix[row,col]=headMean\n",
    "            \n",
    "    matrix[-1,:-1] = modFhd['cellHeadGrid']['lay'+str(lay)][-1,:]\n",
    "    matrix[:-1,-1] = modFhd['cellHeadGrid']['lay'+str(lay)][:,-1]\n",
    "    matrix[-1,-1] = modFhd['cellHeadGrid']['lay'+str(lay)][-1,-1]\n",
    "\n",
    "    vertexHeadGridCentroid['lay'+str(lay)]=matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modFhd['vertexHeadGrid'] = {}\n",
    "for lay in range(modDis['vertexLays']):\n",
    "    anyGrid = vertexHeadGridCentroid\n",
    "    if lay == modDis['cellLays']:\n",
    "        modFhd['vertexHeadGrid']['lay'+str(lay)] = anyGrid['lay'+str(lay-1)]\n",
    "    elif lay == 0:\n",
    "        modFhd['vertexHeadGrid']['lay0'] = anyGrid['lay0']\n",
    "    else:\n",
    "        \n",
    "        value = np.where(anyGrid['lay'+str(lay)]>-1e+20,\n",
    "                         anyGrid['lay'+str(lay)],\n",
    "                         (anyGrid['lay'+str(lay-1)] + anyGrid['lay'+str(lay)])/2\n",
    "                          )\n",
    "        modFhd['vertexHeadGrid']['lay'+str(lay)] = value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the DRN Info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Get the list of drain cells with attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modDrn['DrainNumber'] = int(drnLines[2].split()[0])\n",
    "modDrn['DrainCells'] = []\n",
    "\n",
    "for item in range(modDrn['DrainNumber']-4):\n",
    "    anyLine = drnLines[item + 4].split()\n",
    "    \n",
    "    anyList = [ int(anyLine[0]), int(anyLine[1]), int(anyLine[2]), float(anyLine[3]), float(anyLine[4])]\n",
    "    \n",
    "    modDrn['DrainCells'].append(anyList)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Tables on Cell and Vextex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Water Table on Cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty numpy array for the water table\n",
    "waterTableCellGrid = np.zeros((modDis['cellRows'],modDis['cellCols']))\n",
    "\n",
    "#obtain the first positive or real head from the head array\n",
    "for row in range(modDis['cellRows']):\n",
    "    for col in range(modDis['cellCols']):\n",
    "        anyList = []\n",
    "        for lay in range(modDis['cellLays']):\n",
    "            anyList.append(modFhd['cellHeadGrid']['lay' + str(lay)][row,col])\n",
    "        a = np.asarray(anyList)\n",
    "        if list(a[a>-1e+10]) != []:                    #just in case there are some inactive zones\n",
    "            waterTableCellGrid[row,col] = a[a>-1e+10][0]\n",
    "        else: waterTableCellGrid[row,col] = -1e+10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Water Table on Vertex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty numpy array for the water table\n",
    "waterTableVertexGrid = np.zeros((modDis['vertexRows'],modDis['vertexCols']))\n",
    "\n",
    "#obtain the first positive or real head from the head array\n",
    "for row in range(modDis['vertexRows']):\n",
    "    for col in range(modDis['vertexCols']):\n",
    "        anyList = []\n",
    "        for lay in range(modDis['cellLays']):\n",
    "            anyList.append(modFhd['vertexHeadGrid']['lay' + str(lay)][row,col])\n",
    "        a = np.asarray(anyList)\n",
    "        if list(a[a>-1e+10]) != []:                    #just in case there are some inactive zones\n",
    "            waterTableVertexGrid[row,col] = a[a>-1e+10][0]\n",
    "        else: waterTableVertexGrid[row,col] = -1e+10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Water Table on Cells and Vertex as List"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "listWaterTableCell = list(waterTableCellGrid.flatten())\n",
    "waterTableListVertex = list(waterTableVertexGrid.flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lists for the VTK file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of xyz points for all vertex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty list to store all vertex XYZ\n",
    "vertexXYZPoints = []\n",
    "\n",
    "#definition of xyz points for all vertex\n",
    "for lay in range(modDis['vertexLays']):\n",
    "    for row in range(modDis['vertexRows']):\n",
    "        for col in range(modDis['vertexCols']):\n",
    "            xyz=[\n",
    "                modDis['vertexEasting'][col], \n",
    "                modDis['vertexNorthing'][row],\n",
    "                modDis['vertexZGrid']['lay'+str(lay)][row, col]\n",
    "                ]\n",
    "            vertexXYZPoints.append(xyz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of xyz points for Water Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty list to store all vertex Water Table XYZ\n",
    "vertexWaterTableXYZPoints = []\n",
    "\n",
    "#definition of xyz points for all vertex\n",
    "for row in range(modDis['vertexRows']):\n",
    "    for col in range(modDis['vertexCols']):\n",
    "        if waterTableVertexGrid[row, col] > -1e+10:\n",
    "            waterTable = waterTableVertexGrid[row, col]\n",
    "        else:\n",
    "            waterTable = 1000\n",
    "        xyz=[\n",
    "            modDis['vertexEasting'][col], \n",
    "            modDis['vertexNorthing'][row],\n",
    "            waterTable\n",
    "            ]\n",
    "        vertexWaterTableXYZPoints.append(xyz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of xyz points for Drain Cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty list to store drain vertex XYZ list\n",
    "listDrainCellQuadXYZPoints = []\n",
    "\n",
    "#definition of xyz points for all drain vertex\n",
    "for row in range(modDis['vertexRows']):\n",
    "    for col in range(modDis['vertexCols']):\n",
    "        xyz=[\n",
    "            modDis['vertexEasting'][col], \n",
    "            modDis['vertexNorthing'][row],\n",
    "            modDis['vertexZGrid']['lay0'][row,col]\n",
    "            ]\n",
    "        listDrainCellQuadXYZPoints.append(xyz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of IBound list for all cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty list to store all ibound\n",
    "listIBound = []\n",
    "\n",
    "#definition of IBOUND\n",
    "for lay in range(modDis['cellLays']):\n",
    "    for item in modBas['cellIboundList']['lay'+str(lay)]:\n",
    "        listIBound.append(item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of all Cell Heads and Vertex Heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty list to store all cell heads\n",
    "listCellHead = []\n",
    "\n",
    "#definition of cellHead\n",
    "for lay in range(modDis['cellLays']):\n",
    "    for item in list(modFhd['cellHeadGrid']['lay'+str(lay)].flatten()):\n",
    "        listCellHead.append(item)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#empty list to store all vertex heads\n",
    "listVertexHead = []\n",
    "\n",
    "#definition of vertexHead\n",
    "for lay in range(modDis['vertexLays']):\n",
    "    for item in list(modFhd['vertexHeadGrid']['lay'+str(lay)].flatten()):\n",
    "        listVertexHead.append(item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of Cell Ibound List"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty list to store drain cells IO\n",
    "listDrainsCellsIO = []    \n",
    "    \n",
    "#definition of drain cells\n",
    "modDrn['DrainCellsGrid'] = np.zeros((modDis['cellRows'],modDis['cellCols']))\n",
    "for item in modDrn['DrainCells']:\n",
    "    modDrn['DrainCellsGrid'][item[1],item[2]] = 1\n",
    "listDrainsCellsIO = list(modDrn['DrainCellsGrid'].flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hexahedrons and Quads sequences for the VTK File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List of Layer Quad Sequences (Works only for a single layer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty list to store cell coordinates\n",
    "listLayerQuadSequence = []\n",
    "\n",
    "#definition of hexahedrons cell coordinates\n",
    "for row in range(modDis['cellRows']):\n",
    "    for col in range(modDis['cellCols']):\n",
    "        pt0 = modDis['vertexCols']*(row+1)+col\n",
    "        pt1 = modDis['vertexCols']*(row+1)+col+1\n",
    "        pt2 = modDis['vertexCols']*(row)+col+1\n",
    "        pt3 = modDis['vertexCols']*(row)+col\n",
    "        anyList = [pt0,pt1,pt2,pt3]\n",
    "        listLayerQuadSequence.append(anyList)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List of Hexa Sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#empty list to store cell coordinates\n",
    "listHexaSequence = []\n",
    "\n",
    "#definition of hexahedrons cell coordinates\n",
    "for lay in range(modDis['cellLays']):\n",
    "    for row in range(modDis['cellRows']):\n",
    "        for col in range(modDis['cellCols']):\n",
    "            pt0 = modDis['vertexperlay']*(lay+1)+modDis['vertexCols']*(row+1)+col\n",
    "            pt1 = modDis['vertexperlay']*(lay+1)+modDis['vertexCols']*(row+1)+col+1\n",
    "            pt2 = modDis['vertexperlay']*(lay+1)+modDis['vertexCols']*(row)+col+1\n",
    "            pt3 = modDis['vertexperlay']*(lay+1)+modDis['vertexCols']*(row)+col\n",
    "            pt4 = modDis['vertexperlay']*(lay)+modDis['vertexCols']*(row+1)+col\n",
    "            pt5 = modDis['vertexperlay']*(lay)+modDis['vertexCols']*(row+1)+col+1\n",
    "            pt6 = modDis['vertexperlay']*(lay)+modDis['vertexCols']*(row)+col+1\n",
    "            pt7 = modDis['vertexperlay']*(lay)+modDis['vertexCols']*(row)+col\n",
    "            anyList = [pt0,pt1,pt2,pt3,pt4,pt5,pt6,pt7]\n",
    "            listHexaSequence.append(anyList)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Active Cells Hexahedrons and Quads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell Heads and Hexa Sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "listHexaSequenceDef = []\n",
    "listCellHeadDef = []\n",
    "\n",
    "for i in range(len(listCellHead)):\n",
    "    if listCellHead[i] > -1e-10:\n",
    "        listHexaSequenceDef.append(listHexaSequence[i])\n",
    "        listCellHeadDef.append(listCellHead[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Active Cells and Hexa Sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "listActiveHexaSequenceDef = []\n",
    "listIBoundDef = []\n",
    "\n",
    "#filter hexahedrons and heads for active cells\n",
    "for i in range(len(listIBound)):\n",
    "    if listIBound[i] > 0:\n",
    "        listActiveHexaSequenceDef.append(listHexaSequence[i])\n",
    "        listIBoundDef.append(listIBound[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Water Table Quads and Cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "listWaterTableQuadSequenceDef = []\n",
    "listWaterTableCellDef = []\n",
    "\n",
    "for item in range(len(listWaterTableCell)):\n",
    "    if listWaterTableCell[item] > -1e10:\n",
    "        listWaterTableQuadSequenceDef.append(listLayerQuadSequence[item])\n",
    "        listWaterTableCellDef.append(listWaterTableCell[item])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drain Cells Quads and Cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the secuence and values for cell with drains\n",
    "listDrainsCellsIODef = []    \n",
    "listDrainsCellsSecuenceDef = []\n",
    "\n",
    "for item in range(len(listDrainsCellsIO)):\n",
    "    if listDrainsCellsIO[item] > 0:\n",
    "        listDrainsCellsIODef.append(listDrainsCellsIO[item])\n",
    "        listDrainsCellsSecuenceDef.append(listLayerQuadSequence[item])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VTK creation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary of lists for the vtk creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Point sets\n",
    "#   vertexXYZPoints                    for XYZ in all cell vertex\n",
    "#   vertexWaterTableXYZPoints          for XYZ in all water table quad vertex\n",
    "#   listDrainCellQuadXYZPoints         for XYZ in all drain cells quad vertex\n",
    "\n",
    "### Quad and Hexa secuences\n",
    "#   listHexaSequenceDef                for Head Hexa Sequence in all active cells\n",
    "#   listActiveHexaSequenceDef          for Active Hexa Sequence in all active cells\n",
    "#   listWaterTableQuadSequenceDef      for Water Table Quad Sequence in all active cells\n",
    "#   listDrainsCellsSecuenceDef         for Drain Cell Quad Sequence in drain cells\n",
    "\n",
    "### Cell data\n",
    "#   listCellHeadDef                    for filtered active cells\n",
    "#   listIBoundDef                      \n",
    "#   listWaterTableCellDef              for filtered water table cells\n",
    "#   listDrainsCellsIODef               for filtered drains cells\n",
    "\n",
    "### Point data\n",
    "#   listVertexHead                     for heads in all cells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Heads on Vertex and Cells VTK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "textoVtk = open('VTUFiles/VTU_Heads.vtu','w')\n",
    "\n",
    "#add header\n",
    "textoVtk.write('<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\\n')\n",
    "textoVtk.write('  <UnstructuredGrid>\\n')\n",
    "textoVtk.write('    <Piece NumberOfPoints=\"'+str(len(vertexXYZPoints))+'\" NumberOfCells=\"'+str(len(listHexaSequenceDef))+'\">\\n')  \n",
    "\n",
    "#point data\n",
    "textoVtk.write('      <PointData Scalars=\"Heads\">\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Heads\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listVertexHead)):\n",
    "    textvalue = str(listVertexHead[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + textvalue + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(textvalue + '\\n          ')\n",
    "    else:\n",
    "        textoVtk.write(textvalue + ' ')\n",
    "textoVtk.write('\\n')\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </PointData>\\n')\n",
    "\n",
    "#cell data\n",
    "textoVtk.write('      <CellData Scalars=\"Model\">\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"CellHead\" format=\"ascii\" RangeMin=\"1\" RangeMax=\"1\">\\n')\n",
    "for item  in range(len(listCellHeadDef)): #cell list\n",
    "    textvalue = str(int(listCellHeadDef[item]))\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + textvalue + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(textvalue + '\\n          ')\n",
    "    else:\n",
    "        textoVtk.write(textvalue + ' ')\n",
    "textoVtk.write('\\n')\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </CellData>\\n')\n",
    "\n",
    "#points definition\n",
    "textoVtk.write('      <Points>\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\\n')\n",
    "for item in range(len(vertexXYZPoints)):\n",
    "    tuplevalue = tuple(vertexXYZPoints[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write(\"          %.2f %.2f %.2f \" % tuplevalue)\n",
    "    elif item % 4 == 0:        \n",
    "        textoVtk.write(\"%.2f %.2f %.2f \\n          \" % tuplevalue)\n",
    "    elif item == len(vertexXYZPoints)-1:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \\n\" % tuplevalue)\n",
    "    else:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \" % tuplevalue)   \n",
    "textoVtk.write('        </DataArray>\\n')        \n",
    "textoVtk.write('      </Points>\\n')         \n",
    "    \n",
    "#cell connectivity\n",
    "textoVtk.write('      <Cells>\\n')               \n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listHexaSequenceDef)):\n",
    "    textoVtk.write('          ')\n",
    "    textoVtk.write('%s %s %s %s %s %s %s %s \\n' % tuple(listHexaSequenceDef[item]))\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell offset\n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\\n')\n",
    "for item in range(len(listHexaSequenceDef)):\n",
    "    offset = str((item+1)*8)\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + offset + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(offset + ' \\n          ')\n",
    "    elif item == len(listHexaSequenceDef)-1:\n",
    "        textoVtk.write(offset + ' \\n')\n",
    "    else:\n",
    "        textoVtk.write(offset + ' ')\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell type\n",
    "textoVtk.write('        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\\n') \n",
    "for item in range(len(listHexaSequenceDef)):\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + '12 ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write('12 \\n          ')\n",
    "    elif item == len(listHexaSequenceDef)-1:\n",
    "        textoVtk.write('12 \\n')\n",
    "    else:\n",
    "        textoVtk.write('12 ')        \n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </Cells>\\n')  \n",
    "        \n",
    "#footer\n",
    "textoVtk.write('    </Piece>\\n')\n",
    "textoVtk.write('  </UnstructuredGrid>\\n')\n",
    "textoVtk.write('</VTKFile>\\n')\n",
    "\n",
    "textoVtk.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Active Cell VTK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "textoVtk = open('VTUFiles/VTU_Active.vtu','w')\n",
    "\n",
    "#add header\n",
    "textoVtk.write('<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\\n')\n",
    "textoVtk.write('  <UnstructuredGrid>\\n')\n",
    "textoVtk.write('    <Piece NumberOfPoints=\"'+str(len(vertexXYZPoints))+'\" NumberOfCells=\"'+str(len(listActiveHexaSequenceDef))+'\">\\n')  \n",
    "\n",
    "#cell data\n",
    "textoVtk.write('      <CellData Scalars=\"Model\">\\n')\n",
    "textoVtk.write('        <DataArray type=\"Int32\" Name=\"Active\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listIBoundDef)): #cell list\n",
    "    textvalue = str(int(listIBoundDef[item]))\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + textvalue + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(textvalue + '\\n          ')\n",
    "    else:\n",
    "        textoVtk.write(textvalue + ' ')\n",
    "textoVtk.write('\\n')\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </CellData>\\n')\n",
    "\n",
    "#points definition\n",
    "textoVtk.write('      <Points>\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\\n')\n",
    "for item in range(len(vertexXYZPoints)):\n",
    "    tuplevalue = tuple(vertexXYZPoints[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write(\"          %.2f %.2f %.2f \" % tuplevalue)\n",
    "    elif item % 4 == 0:        \n",
    "        textoVtk.write('%.2f %.2f %.2f \\n          ' % tuplevalue)\n",
    "    elif item == len(vertexXYZPoints)-1:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \\n\" % tuplevalue)\n",
    "    else:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \" % tuplevalue)\n",
    "textoVtk.write('        </DataArray>\\n')        \n",
    "textoVtk.write('      </Points>\\n')         \n",
    "        \n",
    "#cell connectivity\n",
    "textoVtk.write('      <Cells>\\n')               \n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listActiveHexaSequenceDef)):\n",
    "    textoVtk.write('          ')\n",
    "    textoVtk.write('%s %s %s %s %s %s %s %s \\n' % tuple(listActiveHexaSequenceDef[item]))\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell offsets\n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\\n')\n",
    "for item in range(len(listActiveHexaSequenceDef)):\n",
    "    offset = str((item+1)*8)\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + offset + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(offset + ' \\n          ')\n",
    "    elif item == len(listActiveHexaSequenceDef)-1:\n",
    "        textoVtk.write(offset + ' \\n')\n",
    "    else:\n",
    "        textoVtk.write(offset + ' ')\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell types\n",
    "textoVtk.write('        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\\n') \n",
    "for item in range(len(listActiveHexaSequenceDef)):\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + '12 ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write('12 \\n          ')\n",
    "    elif item == len(listActiveHexaSequenceDef)-1:\n",
    "        textoVtk.write('12 \\n')\n",
    "    else:\n",
    "        textoVtk.write('12 ')        \n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </Cells>\\n')  \n",
    "        \n",
    "#footer\n",
    "textoVtk.write('    </Piece>\\n')\n",
    "textoVtk.write('  </UnstructuredGrid>\\n')\n",
    "textoVtk.write('</VTKFile>\\n')\n",
    "\n",
    "textoVtk.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Water Table VTK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "textoVtk = open('VTUFiles/VTU_WaterTable.vtu','w')\n",
    "\n",
    "\n",
    "#add header\n",
    "textoVtk.write('<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\\n')\n",
    "textoVtk.write('  <UnstructuredGrid>\\n')\n",
    "textoVtk.write('    <Piece NumberOfPoints=\"'+str(len(vertexWaterTableXYZPoints))+'\" NumberOfCells=\"'+\n",
    "               str(len(listWaterTableCellDef))+'\">\\n')  \n",
    "\n",
    "#cell data\n",
    "textoVtk.write('      <CellData Scalars=\"Water Table\">\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Heads\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listWaterTableCellDef)):\n",
    "    textvalue = str(listWaterTableCellDef[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + textvalue + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(textvalue + '\\n          ')\n",
    "    else:\n",
    "        textoVtk.write(textvalue + ' ')\n",
    "textoVtk.write('\\n')\n",
    "\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </CellData>\\n')\n",
    "\n",
    "#points definition\n",
    "textoVtk.write('      <Points>\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\\n')\n",
    "for item in range(len(vertexWaterTableXYZPoints)):\n",
    "    tuplevalue = tuple(vertexWaterTableXYZPoints[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write(\"          %.2f %.2f %.2f \" % tuplevalue)\n",
    "    elif item % 4 == 0:        \n",
    "        textoVtk.write('%.2f %.2f %.2f \\n          ' % tuplevalue)\n",
    "    elif item == len(vertexWaterTableXYZPoints)-1:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \\n\" % tuplevalue)\n",
    "    else:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \" % tuplevalue)\n",
    "textoVtk.write('        </DataArray>\\n')        \n",
    "textoVtk.write('      </Points>\\n')         \n",
    "        \n",
    "#cell connectivity\n",
    "textoVtk.write('      <Cells>\\n')               \n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listWaterTableQuadSequenceDef)):\n",
    "    textoVtk.write('          ')\n",
    "    textoVtk.write('%s %s %s %s \\n' % tuple(listWaterTableQuadSequenceDef[item]))\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell offsets\n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\\n')\n",
    "for item in range(len(listWaterTableQuadSequenceDef)):\n",
    "    offset = str((item+1)*4)\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + offset + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(offset + ' \\n          ')\n",
    "    elif item == len(listWaterTableQuadSequenceDef)-1:\n",
    "        textoVtk.write(offset + ' \\n')\n",
    "    else:\n",
    "        textoVtk.write(offset + ' ')\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "#cell types\n",
    "textoVtk.write('        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\\n') \n",
    "for item in range(len(listWaterTableQuadSequenceDef)):\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + '9 ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write('9 \\n          ')\n",
    "    elif item == len(listWaterTableQuadSequenceDef)-1:\n",
    "        textoVtk.write('9 \\n')\n",
    "    else:\n",
    "        textoVtk.write('9 ')\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </Cells>\\n')  \n",
    "        \n",
    "#footer\n",
    "textoVtk.write('    </Piece>\\n')\n",
    "textoVtk.write('  </UnstructuredGrid>\\n')\n",
    "textoVtk.write('</VTKFile>\\n')\n",
    "\n",
    "textoVtk.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Drain Cell VTK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "textoVtk = open('VTUFiles/VTU_DrainCells.vtu','w')\n",
    "\n",
    "#add header\n",
    "textoVtk.write('<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\\n')\n",
    "textoVtk.write('  <UnstructuredGrid>\\n')\n",
    "textoVtk.write('    <Piece NumberOfPoints=\"'+str(len(listDrainCellQuadXYZPoints))\n",
    "               +'\" NumberOfCells=\"'+str(len(listDrainsCellsSecuenceDef))+'\">\\n')  \n",
    "\n",
    "#cell data\n",
    "textoVtk.write('      <CellData Scalars=\"Water Table\">\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Heads\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listDrainsCellsIODef)):\n",
    "    textvalue = str(listDrainsCellsIODef[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + textvalue + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(textvalue + '\\n          ')\n",
    "    else:\n",
    "        textoVtk.write(textvalue + ' ')\n",
    "textoVtk.write('\\n')\n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </CellData>\\n')\n",
    "\n",
    "\n",
    "#points definition\n",
    "textoVtk.write('      <Points>\\n')\n",
    "textoVtk.write('        <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\\n')\n",
    "for item in range(len(listDrainCellQuadXYZPoints)):\n",
    "    tuplevalue = tuple(listDrainCellQuadXYZPoints[item])\n",
    "    if item == 0:\n",
    "        textoVtk.write(\"          %.2f %.2f %.2f \" % tuplevalue)\n",
    "    elif item % 4 == 0:        \n",
    "        textoVtk.write('%.2f %.2f %.2f \\n          ' % tuplevalue)\n",
    "    elif item == len(listDrainCellQuadXYZPoints)-1:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \\n\" % tuplevalue)\n",
    "    else:\n",
    "        textoVtk.write(\"%.2f %.2f %.2f \" % tuplevalue)\n",
    "textoVtk.write('        </DataArray>\\n')        \n",
    "textoVtk.write('      </Points>\\n')         \n",
    "        \n",
    "#cell definition\n",
    "textoVtk.write('      <Cells>\\n')               \n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\\n')\n",
    "for item  in range(len(listDrainsCellsSecuenceDef)):\n",
    "    textoVtk.write('          ')\n",
    "    textoVtk.write('%s %s %s %s \\n' % tuple(listDrainsCellsSecuenceDef[item]))\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "textoVtk.write('        <DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\\n')\n",
    "for item in range(len(listDrainsCellsSecuenceDef)):\n",
    "    offset = str((item+1)*4)\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + offset + ' ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write(offset + ' \\n          ')\n",
    "    elif item == len(listDrainsCellsSecuenceDef)-1:\n",
    "        textoVtk.write(offset + ' \\n')\n",
    "    else:\n",
    "        textoVtk.write(offset + ' ')\n",
    "textoVtk.write('        </DataArray>\\n') \n",
    "textoVtk.write('        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\\n') \n",
    "for item in range(len(listDrainsCellsSecuenceDef)):\n",
    "    if item == 0:\n",
    "        textoVtk.write('          ' + '9 ')\n",
    "    elif item % 20 == 0:\n",
    "        textoVtk.write('9 \\n          ')\n",
    "    elif item == len(listDrainsCellsSecuenceDef)-1:\n",
    "        textoVtk.write('9 \\n')\n",
    "    else:\n",
    "        textoVtk.write('9 ')       \n",
    "textoVtk.write('        </DataArray>\\n')\n",
    "textoVtk.write('      </Cells>\\n')  \n",
    "        \n",
    "#footer\n",
    "textoVtk.write('    </Piece>\\n')\n",
    "textoVtk.write('  </UnstructuredGrid>\\n')\n",
    "textoVtk.write('</VTKFile>\\n')\n",
    "\n",
    "textoVtk.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
