/*****************************************************************************
*
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

// ************************************************************************* //
//  File: avtAxisAlignedSlice4DFilter.C
// ************************************************************************* //

#include <avtAxisAlignedSlice4DFilter.h>

#include <algorithm>
#include <float.h>

#ifdef PARALLEL
#include <avtParallel.h>
#include <mpi.h>
#endif

#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDoubleArray.h>
#include <vtkFieldData.h>
#include <vtkIntArray.h>
#include <vtkRectilinearGrid.h>

#include <vtkDataSetWriter.h>

// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter constructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

avtAxisAlignedSlice4DFilter::avtAxisAlignedSlice4DFilter()
{
    range[0] = 0;
    range[1] = 0;
    spatialExtents[0] = FLT_MAX;
    spatialExtents[1] = -FLT_MAX;
    spatialExtents[2] = FLT_MAX;
    spatialExtents[3] = -FLT_MAX;
    spatialExtents[4] = 0;
    spatialExtents[5] = 0;
}


// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter destructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
//  Modifications:
//
// ****************************************************************************

avtAxisAlignedSlice4DFilter::~avtAxisAlignedSlice4DFilter()
{
}


// ****************************************************************************
//  Method:  avtAxisAlignedSlice4DFilter::Create
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

avtFilter *
avtAxisAlignedSlice4DFilter::Create()
{
    return new avtAxisAlignedSlice4DFilter();
}


// ****************************************************************************
//  Method:      avtAxisAlignedSlice4DFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

void
avtAxisAlignedSlice4DFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const AxisAlignedSlice4DAttributes*)a;
}


// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtAxisAlignedSlice4DFilter with the given
//      parameters would result in an equivalent avtAxisAlignedSlice4DFilter.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

bool
avtAxisAlignedSlice4DFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(AxisAlignedSlice4DAttributes*)a);
}


// ****************************************************************************
//  Method: avAxisAlignedSlice4DFilter::PreExecute
//
//  Purpose:
//      Called before "Execute". Initialize variable we use to keep track of
//      value range.
//
//  Programmer: Gunther H. Weber 
//  Creation:   June 6, 2013
//
//  Modifications:
//
// ****************************************************************************

void
avtAxisAlignedSlice4DFilter::PreExecute(void)
{
    avtPluginDataTreeIterator::PreExecute();

    range[0] =  FLT_MAX;
    range[1] = -FLT_MAX;
    spatialExtents[0] = FLT_MAX;
    spatialExtents[1] = -FLT_MAX;
    spatialExtents[2] = FLT_MAX;
    spatialExtents[3] = -FLT_MAX;
    spatialExtents[4] = 0;
    spatialExtents[5] = 0;

    intVector *sliceIndices[4] = { &atts.GetI(), &atts.GetJ(), &atts.GetK(), &atts.GetL() };
    int nFixed = 0;
    int nVarying = 0;
    for (int i = 0; i < 4; ++i)
    {
        if (sliceIndices[i]->size())
        {
            if (nFixed < 2)
                fixedAxis[nFixed++] = i;
            else
                break;
        }
        else
        {
            if (nVarying < 2)
                varyingAxis[nVarying++] = i;
            else
                break;
        }
    }

    if (nFixed != 2 || nVarying != 2)
    {
        EXCEPTION1(ImproperUseException,
                "Need to specify indices for exactly two dimensions");
    }

    if (sliceIndices[fixedAxis[0]]->size() != sliceIndices[fixedAxis[1]]->size())
    {
        EXCEPTION1(ImproperUseException,
                "Both coordinate arrays need to have the same number of elements.");
    }

    if (sliceIndices[fixedAxis[0]]->size() != 1)
    {
        EXCEPTION1(ImproperUseException, "Currently only one slice supported.");
    }

    for (int i = 0; i < 2; ++i)
    {
        sliceIndex[i] = sliceIndices[fixedAxis[i]]->front();
    }
}

// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter::Execute
//
//  Purpose:
//      Extracts a 2D slice from a 4D data set in a location (or a list of locations)
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
//  Modifications:
//    Kathleen Biagas, Tue Jul  8 14:52:03 PDT 2014
//    Moved avtRealDims code into FIXME comment, to suppress warnings until
//    such a time as they are used.
//
//    Eric Brugger, Thu Jul 24 13:27:16 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtAxisAlignedSlice4DFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    // Ensure that we are working on rectilinear grid
    vtkRectilinearGrid *rgrid = dynamic_cast<vtkRectilinearGrid*>(in_dr->GetDataVTK());
    if (!rgrid)
        EXCEPTION1(ImproperUseException,
                "Can only extract point function for a rectilinear grid.");

    double dx[4] = { 1.0, 1.0, 1.0, 1.0 };
    if (vtkDoubleArray *dx_arr = dynamic_cast<vtkDoubleArray*>(rgrid->GetFieldData()->GetArray("dx_array")))
    {
        dx[2] = dx_arr->GetValue(0);
        dx[3] = dx_arr->GetValue(1);
    }

    // Dims (number of cells along each direction)
    int dims[4];
    rgrid->GetDimensions(dims);
    // Translate from VTK dims (number of vertices) to number of cells
    for (int d=0; d<2; ++d)
        dims[d]--;
    if (vtkIntArray *v_dims_arr = dynamic_cast<vtkIntArray*>(rgrid->GetFieldData()->GetArray("v_dims")))
    {
        for (int i = 0; i < 2; ++i)
            dims[i+2] = v_dims_arr->GetValue(i);
    }
    else
    {
        EXCEPTION1(ImproperUseException,
                "Internal error: Velocity dimensions not set by database plugin.");
    }
    //std::cout << "DIMS: " << dims[0] << " " << dims[1] << " " << dims[2] << " " << dims[3] << std::endl;

    // Base index
    int base_idx[4] = { 0, 0, 0, 0 };
    if (vtkIntArray *base_idx_arr = dynamic_cast<vtkIntArray*>(rgrid->GetFieldData()->GetArray("base_index")))
        for (int d = 0; d < 2; ++d)
            base_idx[d] = base_idx_arr->GetValue(d);
    if (vtkIntArray *v_base_idx_arr = dynamic_cast<vtkIntArray*>(rgrid->GetFieldData()->GetArray("v_base_index")))
        for (int d = 0; d < 2; ++d)
            base_idx[d+2] = v_base_idx_arr->GetValue(d);
    //std::cout << "BASE IDX: " << base_idx[0] << " " << base_idx[1] << " " << base_idx[2] << " " << base_idx[3] << std::endl;

    // Max idx
    int max_idx[4];
    for (int d=0; d<4; ++d)
    {
        max_idx[d] = base_idx[d] + dims[d] - 1;
    }
    //std::cout << "MAX IDX: " << max_idx[0] << " " << max_idx[1] << " " << max_idx[2] << " " << max_idx[3] << std::endl;

    // FIXME: Take avtRealDims into account
    // arr = dynamic_cast<vtkIntArray*>(rgrid->GetFieldData()->GetArray("avtRealDims"));
    // int avtRealDims[6] = { 0, 0, 0, 0, 0, 0 };
    // if (arr)
    //     for (int d = 0; d < 6; ++d)
    //         avtRealDims[d] = arr->GetValue(d);


    if (base_idx[fixedAxis[0]] <= sliceIndex[0] && sliceIndex[0] <= max_idx[fixedAxis[0]] &&
        base_idx[fixedAxis[1]] <= sliceIndex[1] && sliceIndex[1] <= max_idx[fixedAxis[1]])
    {
        const char *justTheVar = pipelineVariable + strlen("operators/AxisAlignedSlice4D/");
        vtkDataArray *data = rgrid->GetCellData()->GetArray(justTheVar);
        if (!data)
        {
            EXCEPTION1(VisItException,
                    "Internal error: Could not get data for array variable (maybe due to an operator requesting ghost zones).");
        }

        vtkRectilinearGrid *ogrid = vtkRectilinearGrid::New();
        ogrid->SetDimensions(dims[varyingAxis[0]]+1, dims[varyingAxis[1]]+1, 1);

        if (varyingAxis[0] == 0)
            ogrid->SetXCoordinates(rgrid->GetXCoordinates());
        else if(varyingAxis[0] == 1)
            ogrid->SetXCoordinates(rgrid->GetYCoordinates());
        else
        {
            vtkDataArray *xCoords = rgrid->GetXCoordinates()->NewInstance();
            xCoords->SetNumberOfTuples(dims[varyingAxis[0]]+1);
            //std::cout << "X: ";
            for (int i = 0; i <= dims[varyingAxis[0]]; ++i)
            {
                //std::cout << (base_idx[varyingAxis[0]] + i) * dx[varyingAxis[0]] << " ";
                xCoords->SetTuple1(i, (base_idx[varyingAxis[0]] + i) * dx[varyingAxis[0]]);
            }
            //std::cout << std::endl;
            ogrid->SetXCoordinates(xCoords);
            xCoords->Delete();
        }

        if (varyingAxis[1] == 0)
            ogrid->SetYCoordinates(rgrid->GetXCoordinates());
        else if(varyingAxis[1] == 1)
            ogrid->SetYCoordinates(rgrid->GetYCoordinates());
        else
        {
            vtkDataArray *yCoords = rgrid->GetYCoordinates()->NewInstance();
            yCoords->SetNumberOfTuples(dims[varyingAxis[1]]+1);
            //std::cout << "Y: ";
            for (int j = 0; j <= dims[varyingAxis[1]]; ++j)
            {
                //std::cout << (base_idx[varyingAxis[1]] + j) * dx[varyingAxis[1]] << " ";
                yCoords->SetTuple1(j, (base_idx[varyingAxis[1]] + j) * dx[varyingAxis[1]]);
            }
            //std::cout << std::endl;
            ogrid->SetYCoordinates(yCoords);
            yCoords->Delete();
        }

        vtkDataArray *odata = data->NewInstance();
        odata->SetName(outVarName.c_str());
        odata->SetNumberOfComponents(1);
        odata->SetNumberOfTuples(dims[varyingAxis[0]]*dims[varyingAxis[1]]);
        //std::cout << dims[varyingAxis[0]]*dims[varyingAxis[1]] << std::endl;
        ogrid->GetCellData()->SetScalars(odata);

        int ijkl[4];
        ijkl[fixedAxis[0]] = sliceIndex[0];
        ijkl[fixedAxis[1]] = sliceIndex[1];
        for (ijkl[varyingAxis[0]] = base_idx[varyingAxis[0]]; ijkl[varyingAxis[0]] <= max_idx[varyingAxis[0]]; ++ijkl[varyingAxis[0]])
            for (ijkl[varyingAxis[1]] = base_idx[varyingAxis[1]]; ijkl[varyingAxis[1]] <= max_idx[varyingAxis[1]]; ++ijkl[varyingAxis[1]])
            {
                //std::cout << ijkl[0] << " " << ijkl[1] << " " << ijkl[2] << " " << ijkl[3] << ": ";
                int in_ijk[3] = { ijkl[0] - base_idx[0], ijkl[1] - base_idx[1], 0 };
                int in_component = (ijkl[2] - base_idx[2])*dims[3] + ijkl[3] - base_idx[3];
                //int in_component = (ijkl[3] - base_idx[3])*dims[2] + ijkl[2] - base_idx[2];
                int out_ijk[3] = { ijkl[varyingAxis[0]] -  base_idx[varyingAxis[0]], ijkl[varyingAxis[1]] -  base_idx[varyingAxis[1]], 0 };
                //std::cout << in_ijk[0] << " " << in_ijk[1] << " " << in_ijk[2] << " comp: " << in_component << " -> " << out_ijk[0] << " " << out_ijk[1] << " " << out_ijk[2] << std::endl;
                double val = data->GetComponent(rgrid->ComputeCellId(in_ijk), in_component);
                odata->SetTuple1(ogrid->ComputeCellId(out_ijk), val);
                if (val < range[0]) range[0] = val;
                if (val > range[1]) range[1] = val;
            }

        if (varyingAxis[0] == 0)
        {
            vtkDataArray *xCoords = rgrid->GetXCoordinates();
            spatialExtents[0] = std::min(spatialExtents[0], xCoords->GetTuple1(0));
            spatialExtents[1] = std::max(spatialExtents[1], xCoords->GetTuple1(xCoords->GetNumberOfTuples()-1));
        }
        else if (varyingAxis[0] == 1)
        {
            vtkDataArray *yCoords = rgrid->GetYCoordinates();
            spatialExtents[0] = std::min(spatialExtents[0], yCoords->GetTuple1(0));
            spatialExtents[1] = std::max(spatialExtents[1], yCoords->GetTuple1(yCoords->GetNumberOfTuples()-1));
        }
        else
        {
            spatialExtents[0] = std::min(spatialExtents[0], base_idx[varyingAxis[0]] * dx[varyingAxis[0]]);
            spatialExtents[1] = std::max(spatialExtents[1], (base_idx[varyingAxis[0]] + dims[varyingAxis[0]]) * dx[varyingAxis[0]]);
        }

        if (varyingAxis[1] == 0)
        {
            vtkDataArray *xCoords = rgrid->GetXCoordinates();
            spatialExtents[2] = std::min(spatialExtents[2], xCoords->GetTuple1(0));
            spatialExtents[3] = std::max(spatialExtents[3], xCoords->GetTuple1(xCoords->GetNumberOfTuples()-1));
        }
        else if (varyingAxis[1] == 1)
        {
            vtkDataArray *yCoords = rgrid->GetYCoordinates();
            spatialExtents[2] = std::min(spatialExtents[2], yCoords->GetTuple1(0));
            spatialExtents[3] = std::max(spatialExtents[3], yCoords->GetTuple1(yCoords->GetNumberOfTuples()-1));
        }
        else
        {
            spatialExtents[2] = std::min(spatialExtents[2], base_idx[varyingAxis[1]] * dx[varyingAxis[1]]);
            spatialExtents[3] = std::max(spatialExtents[3], (base_idx[varyingAxis[1]] + dims[varyingAxis[1]]) * dx[varyingAxis[1]]);
        }
        spatialExtents[4] = 0;
        spatialExtents[5] = 0;
        //std::cout << " " << spatialExtents[0] << " " << spatialExtents[1] << " " << spatialExtents[2] << " " << spatialExtents[3] << " " << spatialExtents[4] << " " << spatialExtents[5] << std::endl;
        odata->Delete();

        avtDataRepresentation *out_dr = new avtDataRepresentation(ogrid, in_dr->GetDomain(), in_dr->GetLabel());
        return out_dr;
    }
    else
    {
        return 0;
    }
}


// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter::PostExecute
//
//  Purpose:
//      Called after "Execute". Update value range.
//
//  Programmer: Gunther H. Weber 
//  Creation:   June 6, 2013
//
//  Modifications:
//
// ****************************************************************************

    void
avtAxisAlignedSlice4DFilter::PostExecute(void)
{
    avtPluginDataTreeIterator::PostExecute();
    if (outVarName != "")
    {
        avtDataAttributes &atts = GetOutput()->GetInfo().GetAttributes();
        atts.GetThisProcsOriginalDataExtents(outVarName.c_str())->Set(range);
        atts.GetThisProcsActualDataExtents(outVarName.c_str())->Set(range);
        atts.GetOriginalSpatialExtents()->Clear();
        atts.GetThisProcsOriginalSpatialExtents()->Set(spatialExtents);
    }
}


// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter::ModifyContract
//
//  Purpose:
//      Creates a contract the removes the operator-created-expression.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

    avtContract_p
avtAxisAlignedSlice4DFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p rv;
    if (strncmp(pipelineVariable, "operators/AxisAlignedSlice4D/", strlen("operators/AxisAlignedSlice4D/")) == 0)
    {
        outVarName = pipelineVariable;
        const char *justTheVar = pipelineVariable + strlen("operators/AxisAlignedSlice4D/");
        avtDataRequest_p dr = new avtDataRequest(in_contract->GetDataRequest(), justTheVar);
        rv = new avtContract(in_contract, dr);
    }
    else
    {
        rv = new avtContract(in_contract);
    }
    return rv;
}


// ****************************************************************************
//  Method: avtAxisAlignedSlice4DFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Tells output that we have a new variable.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jan 7 16:44:23 PST 2016
//
// ****************************************************************************

    void
avtAxisAlignedSlice4DFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &data_atts = GetOutput()->GetInfo().GetAttributes();
    data_atts.SetSpatialDimension(2);
    data_atts.SetTopologicalDimension(2);
    data_atts.AddFilterMetaData("AxisAlignedSlice4D");

    GetOutput()->GetInfo().GetValidity().SetPointsWereTransformed(true);
    GetOutput()->GetInfo().GetAttributes().SetCanUseTransform(false);
    GetOutput()->GetInfo().GetAttributes().SetCanUseInvTransform(false);
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateSpatialMetaData();
    GetOutput()->GetInfo().GetValidity().InvalidateDataMetaData();

    if (outVarName != "")
    {
        avtDataAttributes &atts = GetOutput()->GetInfo().GetAttributes();
        if (! atts.ValidVariable(outVarName))
        {
            atts.AddVariable(outVarName.c_str());
            atts.SetActiveVariable(outVarName.c_str());
            atts.SetVariableDimension(1);
            atts.SetVariableType(AVT_SCALAR_VAR);
            atts.SetCentering(AVT_ZONECENT);
            // you might want to set up extents too
        }
    }
}
