WrappedPython_ParFlow_parflow.h

Go to the documentation of this file.

// Loadable modules
//
// Generated by /builds/gitlab-kitware-sciviz-ci/build/bin/vtkProcessXML-pv6.0
//
#ifndef WrappedPython_ParFlow_parflow_h
#define WrappedPython_ParFlow_parflow_h

#include <cstring>
#include <cassert>
#include <algorithm>


// From file /builds/gitlab-kitware-sciviz-ci/Plugins/ParFlow/parflow.py
static const char* const module_parflow_parflow_string0 =
"\"\"\"This module provides ParaView equivalents to several ParFlow\n"
"pftools algorithms; these are implemented using VTKPythonAlgorithmBase\n"
"subclasses and demonstrate how to access both surface and subsurface\n"
"computational grids within a single filter.\"\"\"\n"
"\n"
"from paraview.util.vtkAlgorithm import *\n"
"import numpy as np\n"
"\n"
"class pftools:\n"
"    \"\"\"This class defines some utilities to perform computations similar\n"
"    to those implemented by ParFlow's pftools suite.\"\"\"\n"
"\n"
"    @staticmethod\n"
"    def dataPointExtent(dataset):\n"
"        \"\"\"Return the parflow extents in (k,j,i) order of the *point*\n"
"        extents of this dataset, which is assumed to be structured.\n"
"\n"
"        Note that VTK normally assumes data is (i,j,k) ordered.\"\"\"\n"
"        ext = tuple(np.flip(dataset.GetDimensions()))\n"
"        return ext\n"
"\n"
"    @staticmethod\n"
"    def dataCellExtent(dataset):\n"
"        \"\"\"Return the parflow extents in (k,j,i) order of the *cell*\n"
"        extents of this dataset, which is assumed to be structured.\n"
"\n"
"        Note that (1) VTK normally uses point counts along axes rather\n"
"        than cell counts and (2) assumes data is (i,j,k) ordered.\"\"\"\n"
"        # Subtracting 1 from GetDimensions accounts for point-counts\n"
"        # vs cell-counts, while np.flip(...) switches the axis order:\n"
"        ext = tuple(np.flip(dataset.GetDimensions()) - 1)\n"
"        return ext\n"
"\n"
"    @staticmethod\n"
"    def computeTopSurface(ext, mask):\n"
"        # I. Reshape the mask and create zk, a column of k-axis indices:\n"
"        mr = np.reshape(mask > 0, ext)\n"
"        zk = np.reshape(np.arange(ext[0]), (ext[0],1,1))\n"
"        # II. The top surface at each (i,j) is the location\n"
"        #      of the highest zk-value that is not masked:\n"
"        top = np.argmax(mr * zk, axis=0)\n"
"        # III. Mark completely-masked columns of cells with an\n"
"        #      invalid \"top\" index of -1:\n"
"        # Columns in mask that had no valid value for mr * zk result\n"
"        # in entries of \"top\" that are 0. But 0 is also a valid k index\n"
"        # so we mark those which should be marked out with -1:\n"
"        mr = (mr == False) # Invert the mask to catch bad cells\n"
"        top = top - np.all(mr, axis=0) # Subtract 1 from masked column entries.\n"
"        return top\n"
"\n"
"    @staticmethod\n"
"    def computeTopSurfaceIndices(top):\n"
"        \"\"\"Convert the \"top\" surface matrix into an array of indices that\n"
"        only include valid top-surface cells. Any negative entry of \"top\"\n"
"        does not have its index included.\n"
"\n"
"        The result is a matrix; each row corresponds to a top-surface cell\n"
"        and each column holds a cell's k-, j-, and i-indices, respectively.\n"
"        \"\"\"\n"
"        itop = np.array([(top[i,j], j, i) \\\n"
"            for i in range(top.shape[0]) \\\n"
"                for j in range(top.shape[1]) \\\n"
"                    if top[i,j] >= 0])\n"
"        return itop\n"
"\n"
"    @staticmethod\n"
"    def computeSurfaceStorage(ext, itop, xx, pressure):\n"
"        # I. Compute the dx * dy area term using point coordinates from the mesh:\n"
"        xt = xx[itop[:,0], itop[:,1], itop[:,2]]\n"
"        dx = xx[itop[:,0], itop[:,1], itop[:,2] + 1] - xt\n"
"        dy = xx[itop[:,0], itop[:,1] + 1, itop[:,2]] - xt\n"
"        # Now dx and dy are matrices whose rows are vectors.\n"
"        # Compute the norms of the vectors.\n"
"        #  dx = np.sqrt(np.sum(dx * dx,axis=1))\n"
"        #  dy = np.sqrt(np.sum(dy * dy,axis=1))\n"
"        # To more closely match parflow:\n"
"        dx = dx[:,0]\n"
"        dy = dy[:,1]\n"
"\n"
"        pp = np.reshape(pressure, ext)\n"
"        pt = pp[itop[:,0], itop[:,1], itop[:,2]]\n"
"        sus = pt * (pt > 0.0) * dx * dy\n"
"        return sus\n"
"\n"
"    @staticmethod\n"
"    def computeSubsurfaceStorage(saturation, pressure, volume, porosity, specific_storage):\n"
"        \"\"\"Compute the subsurface storage for the entire domain.\"\"\"\n"
"        import numpy as np\n"
"        sbs = saturation * volume * (porosity + pressure * specific_storage)\n"
"        return sbs\n"
"\n"
"    @staticmethod\n"
"    def computeGroundwaterStorage(saturation, pressure, volume, porosity, specific_storage):\n"
"        \"\"\"Compute the groundwater storage.\n"
"        This is the same calculation as subsurface storage, but only\n"
"        sums values over cells with a saturation of 1.0.\n"
"        \"\"\"\n"
"        import numpy as np\n"
"        sbs = (saturation == 1.0) * volume * (porosity + pressure * specific_storage)\n"
"        return sbs\n"
"\n"
"    @staticmethod\n"
"    def computeSurfaceRunoff(top, xx, pressure, slope, mannings):\n"
"        \"\"\"Compute surface runoff (water leaving the domain boundary\n"
"        or flowing into a masked area)\"\"\"\n"
"        def addToRunoff(runoff, sro, top, slope, pressure, mannings, delta):\n"
"            \"\"\"Given a truthy 2-d array (idx) of places where runoff occurs,\n"
"            compute the runoff and add it to the total (sro).\n"
"            \"\"\"\n"
"            # We would like to do this:\n"
"            #   sro += \\\n"
"            #       runoff * np.sqrt(np.abs(slope)) / mannings * \\\n"
"            #       np.power(ptop, 5.0/3.0) * delta\n"
"            # ... but we can't since it might result in NaN values where\n"
"            # runoff is False. Instead, only compute the runoff exactly\n"
"            # at locations where runoff is true:\n"
"            ww = np.where(runoff)\n"
"            sro[ww] += \\\n"
"                np.sqrt(np.abs(slope[ww])) / mannings[ww] * \\\n"
"                np.power(ptop[ww], 5.0/3.0) * delta[ww]\n"
"\n"
"        # Initialize runoff:\n"
"        sro = np.zeros(top.shape)\n"
"        # Prepare indices and tempororary variables:\n"
"        sz = np.prod(top.shape)\n"
"        ext = (int(np.prod(pressure.shape)/sz),) + top.shape\n"
"        tk = np.reshape(top, sz)\n"
"        tj = np.floor(np.arange(sz)/top.shape[1]).astype(int)\n"
"        ti = np.arange(sz) % top.shape[1]\n"
"        tt = np.zeros(top.shape) # temporary array\n"
"        # Subset of pressure at the top surface:\n"
"        ptop = np.reshape(np.reshape(pressure, ext)[tk, tj, ti], top.shape)\n"
"\n"
"        # Determine size of top-surface cells along north-south direction:\n"
"        delta = np.reshape((xx[tk, tj + 1, ti] - xx[tk, tj, ti])[:,1], top.shape)\n"
"\n"
"        # Fill the temporary array with data shifted south by 1:\n"
"        np.concatenate((top[1:,:], -np.ones((1,top.shape[1]))), axis=0, out=tt)\n"
"        # Compute conditions for flow exiting to the north:\n"
"        cnorth = (top >= 0) & (tt < 0) & (slope[:,:,1] < 0) & (ptop > 0)\n"
"        # Now add values to per-cell runoff sro:\n"
"        addToRunoff(cnorth, sro, top, slope[:,:,1], pressure, mannings, delta)\n"
"\n"
"        # Fill the temporary array with data shifted north by 1:\n"
"        np.concatenate((-np.ones((1,top.shape[1])), top[0:-1,:]), axis=0, out=tt)\n"
"        # Compute conditions for flow exiting to the south:\n"
"        csouth = (top >= 0) & (tt < 0) & (slope[:,:,1] > 0) & (ptop > 0)\n"
"        # Now add values to per-cell runoff sro:\n"
"        addToRunoff(csouth, sro, top, slope[:,:,1], pressure, mannings, delta)\n"
"\n"
"        # Determine size of top-surface cells along east-west direction:\n"
"        delta = np.reshape((xx[tk, tj, ti + 1] - xx[tk, tj, ti])[:,0], top.shape)\n"
"\n"
"        # Fill the temporary array with data shifted east by 1:\n"
"        np.concatenate((top[:,1:], -np.ones((top.shape[0], 1))), axis=1, out=tt)\n"
"        # Compute conditions for flow exiting to the west:\n"
"        cwest = (top >= 0) & (tt < 0) & (slope[:,:,0] < 0) & (ptop > 0)\n"
"        # Now add values to per-cell runoff sro:\n"
"        addToRunoff(cwest, sro, top, slope[:,:,0], pressure, mannings, delta)\n"
"\n"
"        # Fill the temporary array with data shifted north by 1:\n"
"        np.concatenate((-np.ones((top.shape[0], 1)), top[:,0:-1]), axis=1, out=tt)\n"
"        # Compute conditions for flow exiting to the south:\n"
"        ceast = (top >= 0) & (tt < 0) & (slope[:,:,0] > 0) & (ptop > 0)\n"
"        # Now add values to per-cell runoff sro:\n"
"        addToRunoff(ceast, sro, top, slope[:,:,0], pressure, mannings, delta)\n"
"\n"
"        return sro\n"
"\n"
"    @staticmethod\n"
"    def computeWaterTableDepth(top, saturation, xx):\n"
"        \"\"\"Compute the depth to the water table.\n"
"\n"
"        The returned array is the distance along the z axis between\n"
"        the water surface and ground level. It should always be\n"
"        negative (i.e., it is a displacement on the z axis starting\n"
"        at the surface).\n"
"        \"\"\"\n"
"        sz = np.prod(top.shape)\n"
"        ext = (int(np.prod(saturation.shape)/sz),) + top.shape\n"
"        sr = np.reshape(saturation >= 1.0, ext)\n"
"        zk = np.reshape(np.arange(ext[0]), (ext[0],1,1))\n"
"        # The top surface at each (i,j) is the location\n"
"        # of the highest zk-value that is not masked:\n"
"        depthIdx = np.argmax(sr * zk, axis=0)\n"
"        tk = np.reshape(top, sz)\n"
"        dk = np.reshape(depthIdx, sz)\n"
"        dj = np.floor(np.arange(sz)/top.shape[1]).astype(int)\n"
"        di = np.arange(sz) % top.shape[1]\n"
"        depth = xx[dk, dj, di, 2] - xx[tk, dj, di, 2]\n"
"\n"
"        # Mark locations where the domain is masked or\n"
"        # no water table appears as NaN (not-a-number):\n"
"        sr = (sr == False) # Invert the saturation mask\n"
"        ww = np.where(np.reshape(np.all(sr, axis=0), sz))\n"
"        depth[ww] = np.nan\n"
"        return depth\n"
"\n"
"class ParFlowAlgorithm(VTKPythonAlgorithmBase):\n"
"    \"\"\"A base class for algorithms that operate on\n"
"    ParFlow simulation data.\n"
"\n"
"    This class provides method implementations that\n"
"    simplify working with data from the vtkParFlowMetaReader.\n"
"    Specifically,\n"
"\n"
"    1. RequestUpdateExtent will ensure valid extents\n"
"    are passed to each of the filter's inputs based on\n"
"    their available whole extents. Without this, passing\n"
"    both the surface and subsurface meshes to a python\n"
"    filter would cause warnings as the default\n"
"    implementation simply mirrors the downstream request\n"
"    to each of its upstream filters.\n"
"\n"
"    2. RequestDataObject will create the same type of\n"
"    output as the input data objects so that either\n"
"    image data or structured grids (both of which may\n"
"    be output by the reader) can be passed through the\n"
"    filter.\n"
"    \"\"\"\n"
"    def FillInputPortInformation(self, port, info):\n"
"        info.Set(self.INPUT_REQUIRED_DATA_TYPE(), \"vtkDataSet\")\n"
"        return 1\n"
"\n"
"    def RequestDataObject(self, request, inInfoVec, outInfoVec):\n"
"        \"\"\"Create a data object of the same type as the input.\"\"\"\n"
"        from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData\n"
"        from vtkmodules import vtkCommonDataModel as dm\n"
"        input0 = vtkDataSet.GetData(inInfoVec[0], 0)\n"
"        opt = dm.vtkDataObject.GetData(outInfoVec)\n"
"\n"
"        if opt and opt.IsA(input0.GetClassName()):\n"
"            return 1\n"
"\n"
"        opt = input0.NewInstance()\n"
"        outInfoVec.GetInformationObject(0).Set(dm.vtkDataObject.DATA_OBJECT(), opt)\n"
"        return 1\n"
"\n"
"    def RequestUpdateExtent(self, request, inInfoVec, outInfoVec):\n"
"        \"\"\"Override the default algorithm for updating extents to handle the\n"
"        surface and subsurface models, which have different dimensionality.\n"
"\n"
"        Take the requested (downstream) extent and intersect it with the\n"
"        available (upstream) extents; use that as the request rather than\n"
"        blindly copying the request upstream.\n"
"        \"\"\"\n"
"        from vtkmodules.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline as sddp\n"
"\n"
"        # Determine the port requesting an update from us:\n"
"        reqPort = request.Get(sddp.FROM_OUTPUT_PORT())\n"
"        # Get that port's information and the extents it is requesting:\n"
"        outInfo = self.GetOutputInformation(reqPort)\n"
"        esrc = outInfo.Get(sddp.UPDATE_EXTENT())\n"
"\n"
"        # Slice the extents into tuples holding the lo and hi indices:\n"
"        losrc = esrc[0:6:2]\n"
"        hisrc = esrc[1:6:2]\n"
"        np = self.GetNumberOfInputPorts()\n"
"        for port in range(np):\n"
"            # Loop over all the input connections and set their update extents:\n"
"            nc = self.GetNumberOfInputConnections(port)\n"
"            for connection in range(nc):\n"
"                inInfo = self.GetInputInformation(port, connection)\n"
"                # Set UPDATE_EXTENT to be the intersection of the input\n"
"                # port's WHOLE_EXTENT with the output port's UPDATE_EXTENT.\n"
"                #\n"
"                # NB/FIXME: Really this is incorrect. If reqPort is 1, then\n"
"                #   someone is asking for an update of the 2-d surface grid\n"
"                #   and they could conceivably need the entire 3-d subsurface\n"
"                #   grid that overlaps the surface grid in x and y...\n"
"                etgt = inInfo.Get(sddp.WHOLE_EXTENT())\n"
"                lotgt = etgt[0:6:2]\n"
"                hitgt = etgt[1:6:2]\n"
"                lodst = [int(max(x)) for x in zip(losrc, lotgt)]\n"
"                hidst = [int(min(x)) for x in zip(hisrc, hitgt)]\n"
"                edst = [val for tup in zip(lodst, hidst) for val in tup]\n"
"                inInfo.Set(sddp.UPDATE_EXTENT(), tuple(edst), 6)\n"
"        return 1\n"
"\n"
"@smproxy.filter(label='Subsurface Storage')\n"
"@smhint.xml(\"\"\"\n"
"    <ShowInMenu category=\"ParFlow\"/>\n"
"    <RepresentationType port=\"0\" view=\"RenderView\" type=\"Surface\" />\n"
"\"\"\")\n"
"@smproperty.input(name=\"Subsurface\", port_index=0)\n"
"@smdomain.datatype(dataTypes=[\"vtkDataSet\"], composite_data_supported=False)\n"
"class SubsurfaceStorage(ParFlowAlgorithm):\n"
"    \"\"\"This algorithm computes the subsurface storage across the entire domain.\n"
"    The result is added as field data and marked as time-varying so it\n"
"    can be plotted as a timeseries.\n"
"    \"\"\"\n"
"    def __init__(self):\n"
"        VTKPythonAlgorithmBase.__init__(self, nInputPorts=1, nOutputPorts=1, outputType=\"vtkDataSet\")\n"
"\n"
"    def RequestData(self, request, inInfoVec, outInfoVec):\n"
"        from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData\n"
"        from vtkmodules.vtkIOExodus import vtkExodusIIReader as e2r\n"
"        from vtkmodules.vtkFiltersVerdict import vtkMeshQuality as mq\n"
"        from vtkmodules.util.numpy_support import vtk_to_numpy as vton\n"
"        from vtkmodules.util.numpy_support import numpy_to_vtk as ntov\n"
"        import numpy as np\n"
"\n"
"        ## Fetch the filter inputs and outputs:\n"
"        input0 = vtkDataSet.GetData(inInfoVec[0], 0)\n"
"        output = vtkDataSet.GetData(outInfoVec, 0)\n"
"        output.ShallowCopy(input0)\n"
"\n"
"        ## Fetch input arrays\n"
"        cd = output.GetCellData()\n"
"        vmask = cd.GetArray('mask')\n"
"        vsaturation = cd.GetArray('saturation')\n"
"        vporosity = cd.GetArray('porosity')\n"
"        vpressure = cd.GetArray('pressure')\n"
"        vspecstor = cd.GetArray('specific storage')\n"
"        if vmask == None or vsaturation == None or \\\n"
"           vporosity == None or vpressure == None or \\\n"
"           vspecstor == None:\n"
"            print('Error: A required array was not present.')\n"
"            return 0\n"
"\n"
"        ## Compute the volume of each cell:\n"
"        mqf = mq()\n"
"        mqf.SetHexQualityMeasureToVolume()\n"
;

// From file /builds/gitlab-kitware-sciviz-ci/Plugins/ParFlow/parflow.py
static const char* const module_parflow_parflow_string1 =
"        mqf.SetInputDataObject(0, output)\n"
"        mqf.Update()\n"
"        volume = vton(mqf.GetOutputDataObject(0).GetCellData().GetArray('Quality'))\n"
"\n"
"        ## Get NumPy versions of each array so we can do arithmetic:\n"
"        mask = vton(vmask)\n"
"        saturation = vton(vsaturation)\n"
"        porosity = vton(vporosity)\n"
"        pressure = vton(vpressure)\n"
"        specstor = vton(vspecstor)\n"
"\n"
"        ## Compute the subsurface storage as sbs\n"
"        ## and store it as field data.\n"
"        sbs = np.sum(pftools.computeSubsurfaceStorage( \\\n"
"                saturation, pressure, volume, porosity, specstor))\n"
"        arr = ntov(sbs)\n"
"        arr.SetName('subsurface storage')\n"
"        arr.GetInformation().Set(e2r.GLOBAL_TEMPORAL_VARIABLE(), 1)\n"
"        output.GetFieldData().AddArray(arr)\n"
"\n"
"        return 1\n"
"\n"
"\n"
"@smproxy.filter(label='Water Table Depth')\n"
"@smhint.xml(\"\"\"\n"
"    <ShowInMenu category=\"ParFlow\"/>\n"
"    <RepresentationType port=\"0\" view=\"RenderView\" type=\"Surface\" />\n"
"\"\"\")\n"
"@smproperty.input(name=\"Surface\", port_index=1)\n"
"@smdomain.datatype(dataTypes=[\"vtkDataSet\"], composite_data_supported=False)\n"
"@smproperty.input(name=\"Subsurface\", port_index=0)\n"
"@smdomain.datatype(dataTypes=[\"vtkDataSet\"], composite_data_supported=False)\n"
"class WaterTableDepth(ParFlowAlgorithm):\n"
"    \"\"\"This algorithm computes the depth to the water table the domain surface.\n"
"    \"\"\"\n"
"    def __init__(self):\n"
"        VTKPythonAlgorithmBase.__init__(self, nInputPorts=2, nOutputPorts=1, outputType=\"vtkDataSet\")\n"
"\n"
"    def RequestInformation(self, request, inInfoVec, outInfoVec):\n"
"        \"\"\"Tell ParaView our output extents match the surface, not the subsurface.\"\"\"\n"
"        from vtkmodules.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline as sddp\n"
"\n"
"        iin = inInfoVec[1].GetInformationObject(0)\n"
"        outInfoVec.GetInformationObject(0).Set(sddp.WHOLE_EXTENT(), iin.Get(sddp.WHOLE_EXTENT()), 6);\n"
"        return 1\n"
"\n"
"    def RequestData(self, request, inInfoVec, outInfoVec):\n"
"        from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData\n"
"        from vtkmodules.vtkIOExodus import vtkExodusIIReader as e2r\n"
"        from vtkmodules.vtkFiltersVerdict import vtkMeshQuality as mq\n"
"        from vtkmodules.util.numpy_support import vtk_to_numpy as vton\n"
"        from vtkmodules.util.numpy_support import numpy_to_vtk as ntov\n"
"        import numpy as np\n"
"\n"
"        ## Fetch the filter inputs and outputs:\n"
"        input0 = vtkDataSet.GetData(inInfoVec[0], 0)\n"
"        input1 = vtkDataSet.GetData(inInfoVec[1], 0)\n"
"        output = vtkDataSet.GetData(outInfoVec, 0)\n"
"        output.ShallowCopy(input1)\n"
"\n"
"        ## Fetch input arrays\n"
"        cd = output.GetCellData()\n"
"        scd = input0.GetCellData()\n"
"        vmask = scd.GetArray('mask')\n"
"        vsaturation = scd.GetArray('saturation')\n"
"        if vmask == None or vsaturation == None:\n"
"            print('Error: A required array was not present.')\n"
"            return 0\n"
"\n"
"        ## Get NumPy versions of each array so we can do arithmetic:\n"
"        mask = vton(vmask)\n"
"        saturation = vton(vsaturation)\n"
"\n"
"        ## Find the top surface\n"
"        ext = pftools.dataCellExtent(input0)\n"
"        top = pftools.computeTopSurface(ext, mask)\n"
"\n"
"        ## Compute the water table depth storage as wtd\n"
"        ## and store it as cell data.\n"
"        ps = pftools.dataPointExtent(input0) + (3,)\n"
"        xx = np.reshape(vton(input0.GetPoints().GetData()), ps)\n"
"        wtd = pftools.computeWaterTableDepth(top, saturation, xx)\n"
"        arr = ntov(wtd)\n"
"        arr.SetName('water table depth')\n"
"        output.GetCellData().AddArray(arr)\n"
"\n"
"        return 1\n"
"\n"
"\n"
"@smproxy.filter(label='Water Balance')\n"
"@smhint.xml(\"\"\"\n"
"    <ShowInMenu category=\"ParFlow\"/>\n"
"    <RepresentationType port=\"0\" view=\"RenderView\" type=\"Surface\" />\n"
"\"\"\")\n"
"@smproperty.input(name=\"Surface\", port_index=1)\n"
"@smdomain.datatype(dataTypes=[\"vtkDataSet\"], composite_data_supported=False)\n"
"@smproperty.input(name=\"Subsurface\", port_index=0)\n"
"@smdomain.datatype(dataTypes=[\"vtkDataSet\"], composite_data_supported=False)\n"
"class WaterBalance(ParFlowAlgorithm):\n"
"    \"\"\"This algorithm computes the subsurface storage across the entire domain.\n"
"    The result is added as field data and marked as time-varying so it\n"
"    can be plotted as a timeseries.\n"
"    \"\"\"\n"
"    def __init__(self):\n"
"        VTKPythonAlgorithmBase.__init__(self, nInputPorts=2, nOutputPorts=1, outputType=\"vtkDataSet\")\n"
"\n"
"    def RequestData(self, request, inInfoVec, outInfoVec):\n"
"        from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData\n"
"        from vtkmodules.vtkIOExodus import vtkExodusIIReader as e2r\n"
"        from vtkmodules.vtkFiltersVerdict import vtkMeshQuality as mq\n"
"        from vtkmodules.util.numpy_support import vtk_to_numpy as vton\n"
"        from vtkmodules.util.numpy_support import numpy_to_vtk as ntov\n"
"        import numpy as np\n"
"\n"
"        ## Fetch the filter inputs and outputs:\n"
"        input0 = vtkDataSet.GetData(inInfoVec[0], 0)\n"
"        input1 = vtkDataSet.GetData(inInfoVec[1], 0)\n"
"        output = vtkDataSet.GetData(outInfoVec, 0)\n"
"        output.ShallowCopy(input0)\n"
"\n"
"        ## Fetch input arrays\n"
"        cd = output.GetCellData()\n"
"        scd = input1.GetCellData()\n"
"        vmask = cd.GetArray('mask')\n"
"        vsaturation = cd.GetArray('saturation')\n"
"        vporosity = cd.GetArray('porosity')\n"
"        vpressure = cd.GetArray('pressure')\n"
"        vspecstor = cd.GetArray('specific storage')\n"
"        vslope = scd.GetArray('slope')\n"
"        vmannings = scd.GetArray('mannings')\n"
"        if vmask == None or vsaturation == None or \\\n"
"           vporosity == None or vpressure == None or \\\n"
"           vspecstor == None or vslope == None or \\\n"
"           vmannings == None:\n"
"            print('Error: A required array was not present.')\n"
"            return 0\n"
"\n"
"        ## Compute the volume of each cell:\n"
"        mqf = mq()\n"
"        mqf.SetHexQualityMeasureToVolume()\n"
"        mqf.SetInputDataObject(0, output)\n"
"        mqf.Update()\n"
"        volume = vton(mqf.GetOutputDataObject(0).GetCellData().GetArray('Quality'))\n"
"\n"
"        ## Get NumPy versions of each array so we can do arithmetic:\n"
"        mask = vton(vmask)\n"
"        saturation = vton(vsaturation)\n"
"        porosity = vton(vporosity)\n"
"        pressure = vton(vpressure)\n"
"        specstor = vton(vspecstor)\n"
"        slope = vton(vslope)\n"
"        mannings = vton(vmannings)\n"
"\n"
"        ## Compute the subsurface storage as sbs\n"
"        ## and store it as field data.\n"
"        sbs = np.sum(pftools.computeSubsurfaceStorage( \\\n"
"                saturation, pressure, volume, porosity, specstor))\n"
"        arr = ntov(sbs)\n"
"        arr.SetName('subsurface storage')\n"
"        arr.GetInformation().Set(e2r.GLOBAL_TEMPORAL_VARIABLE(), 1)\n"
"        output.GetFieldData().AddArray(arr)\n"
"\n"
"        ## Find the top surface\n"
"        ext = pftools.dataCellExtent(input0)\n"
"        top = pftools.computeTopSurface(ext, mask)\n"
"        itop = pftools.computeTopSurfaceIndices(top)\n"
"\n"
"        ## Compute the surface storage\n"
"        ps = pftools.dataPointExtent(input0) + (3,)\n"
"        xx = np.reshape(vton(output.GetPoints().GetData()), ps)\n"
"        sus = np.sum(pftools.computeSurfaceStorage(ext, itop, xx, pressure))\n"
"        arr = ntov(sus)\n"
"        arr.SetName('surface storage')\n"
"        arr.GetInformation().Set(e2r.GLOBAL_TEMPORAL_VARIABLE(), 1)\n"
"        output.GetFieldData().AddArray(arr)\n"
"\n"
"        ## Compute the surface runoff\n"
"        slope = np.reshape(slope, top.shape + (2,))\n"
"        mannings = np.reshape(mannings, top.shape)\n"
"        sro = np.sum(pftools.computeSurfaceRunoff(top, xx, pressure, slope, mannings))\n"
"        arr = ntov(sro)\n"
"        arr.SetName('surface runoff')\n"
"        arr.GetInformation().Set(e2r.GLOBAL_TEMPORAL_VARIABLE(), 1)\n"
"        output.GetFieldData().AddArray(arr)\n"
"\n"
"        return 1\n"
"\n"
"if __name__ == \"__main__\":\n"
"    from paraview.detail.pythonalgorithm import get_plugin_xmls\n"
"    from xml.dom.minidom import parseString\n"
"    for xml in get_plugin_xmls(globals()):\n"
"        dom = parseString(xml)\n"
"        print(dom.toprettyxml(\" \",\"\\n\"))\n"
"\n";
// Get single string
inline char* module_parflow_parflow_source()
{

  const size_t len0 = strlen(module_parflow_parflow_string0);
  const size_t len1 = strlen(module_parflow_parflow_string1);
  size_t len = ( 0
    + len0
    + len1 );
  char* res = new char[ len + 1];
  size_t offset = 0;
  std::copy(module_parflow_parflow_string0, module_parflow_parflow_string0 + len0, res + offset); offset += len0;
  std::copy(module_parflow_parflow_string1, module_parflow_parflow_string1 + len1, res + offset); offset += len1;
  assert(offset == len);
  res[offset] = 0;
  return res;
}



#endif