/*****************************************************************************
*
* 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: avtIsovolumeFilter.C
// ************************************************************************* //

#include <avtIsovolumeFilter.h>

#include <vtkVisItClipper.h>
#include <vtkDataSet.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkPolyData.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellDataToPointData.h>

#include <avtAccessor.h>
#include <avtIntervalTree.h>
#include <avtMetaData.h>

#include <DebugStream.h>
#include <VisItException.h>

#include <float.h>

#include <string>
#include <vector>


// ****************************************************************************
//  Method: avtIsovolumeFilter constructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

avtIsovolumeFilter::avtIsovolumeFilter()
{
}


// ****************************************************************************
//  Method: avtIsovolumeFilter destructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
//  Modifications:
//
// ****************************************************************************

avtIsovolumeFilter::~avtIsovolumeFilter()
{
}


// ****************************************************************************
//  Method:  avtIsovolumeFilter::Create
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

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


// ****************************************************************************
//  Method:      avtIsovolumeFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: Jeremy Meredith
//  Creation:   January 30, 2004
//
// ****************************************************************************

void
avtIsovolumeFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const IsovolumeAttributes*)a;

    // We need to specify that we want a secondary variable as soon as 
    // possible.
    if (strcmp(atts.GetVariable().c_str(), "default") != 0)
    {
        SetActiveVariable(atts.GetVariable().c_str());
    }
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtIsovolumeFilter with the given
//      parameters would result in an equivalent avtIsovolumeFilter.
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

bool
avtIsovolumeFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(IsovolumeAttributes*)a);
}


// ****************************************************************************
//  Method:  avtIsovolumeFilter::ExecuteSingleClip
//
//  Purpose:
//    Clip against either a low or high scalar.  Convert cell to point
//    data as necessary.
//
//  Arguments:
//    in_ds      the data set to clip
//    val        the value to clip against
//    flip       false if val is a max, true if val is a min
//
//  Programmer:  Jeremy Meredith
//  Creation:    May  6, 2004
//
//  Modifications:
//
//    Hank Childs, Sat Sep 29 11:24:12 PDT 2007
//    Pass in vtkDataArrays to the clipper, not "float *".
//
//    Eric Brugger, Wed Jul 30 18:49:46 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

vtkDataSet *
avtIsovolumeFilter::ExecuteSingleClip(vtkDataSet *in_ds, float val, bool flip)
{
    vtkVisItClipper *clipper = vtkVisItClipper::New();
    clipper->SetInsideOut(flip);

    //
    // Get the scalar array we'll use for clipping; it must be nodal
    //
    vtkCellDataToPointData *cd2pd = NULL;
    if (in_ds->GetPointData()->GetScalars())
    {
        vtkDataArray *s = in_ds->GetPointData()->GetScalars();
        clipper->SetClipScalars(s, val);
    }
    else if (in_ds->GetCellData()->GetScalars())
    {
        //
        // Okay, our active variable was cell-centered.  Recenter it....
        //
        vtkDataSet *temp_ds = (vtkDataSet *) in_ds->NewInstance();
        temp_ds->CopyStructure(in_ds);
        temp_ds->GetCellData()->SetScalars(in_ds->GetCellData()->GetScalars());

        cd2pd = vtkCellDataToPointData::New();
        cd2pd->SetInputData(temp_ds);
        cd2pd->Update();

        vtkDataSet *temporary = cd2pd->GetOutput();

        // Now tell the clipper about it....
        vtkDataArray *s = temporary->GetPointData()->GetScalars();
        clipper->SetClipScalars(s, val);

        // Wait until after the clipping is done to delete 'cd2pd' (which
        // will take 'temporary' with it)
        temp_ds->Delete();
    }
    else
    {
        debug1 << "Could not find any data for isovolume operation\n";
        clipper->Delete();
        EXCEPTION1(VisItException, "No variable was present for the Isovolume");
    }

    //
    // Do the clipping!
    //
    clipper->SetInputData(in_ds);
    clipper->Update();
    vtkDataSet *out_ds = clipper->GetOutput();
    out_ds->Register(NULL);
    clipper->Delete();

    //
    // Free the temporary filter used to convert to point data
    //
    if (cd2pd)
        cd2pd->Delete();

    return out_ds;
}

// ****************************************************************************
// Method: IsovolumeMinMax
//
// Purpose: 
//   Determine the min and max for the input data array using an accessor.
//
// Arguments:
//
// Returns:    
//
// Note:       
//
// Programmer: Brad Whitlock
// Creation:   Sun Apr 22 01:21:24 PDT 2012
//
// Modifications:
//   
// ****************************************************************************

template <typename Accessor>
inline void IsovolumeMinMax(double &min, double &max, Accessor access)
{
    min = +FLT_MAX;
    max = -FLT_MAX;
    access.InitTraversal();
    while(access.Iterating())
    {
        double value = access.GetTuple1();
        min = (min < value ? min : value);
        max = (max > value ? max : value);
        ++access;
    }
}

// ****************************************************************************
//  Method: avtIsovolumeFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the Isovolume filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: Jeremy Meredith
//  Creation:   January 30, 2004
//
//  Modifications:
//    Jeremy Meredith, Mon Feb  2 13:13:05 PST 2004
//    Fixed memory leak.
//
//    Jeremy Meredith, Wed May  5 14:56:35 PDT 2004
//    Removed the "3D" from the end of vtkVisItClipper because I made it
//    fully support 2D as well.  I also changed it to only take a single
//    cutoff so the math is more robust, and that required making a "min"
//    pass as well as a "max" pass, so that had to change in this function.
//
//    Jeremy Meredith, Thu May  6 11:38:21 PDT 2004
//    I neglected to create a new scalar array in the event that we do both
//    a min and max pass, so I split most of this routine into a new function
//    and called the new one (ExecuteSingleClip) twice.
//
//    Hank Childs, Wed Oct 20 16:16:15 PDT 2004
//    Basic optimizations -- only do min pass and max pass when the data
//    requires it.  Also pass the dataset through if it is wholly contained
//    within the range and return a NULL dataset if it is outside the range.
//
//    Hank Childs, Wed Nov 17 11:43:53 PST 2004
//    At the end, we try to convert poly-data input back to poly-data output.
//    There was an assumption that after processing the data, it would be
//    in unstructured grid form, which is not true.  Add a check here. ['5640]
//
//    Brad Whitlock, Sun Apr 22 01:21:10 PDT 2012
//    I added double support.
//
//    Gunther H. Weber, Tue May 13 11:51:19 PDT 2014
//    Check whether data set is empty before performing max clip to prevent
//    crash in border case.
//
//    Eric Brugger, Wed Jul 30 18:49:46 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Brad Whitlock, Fri Sep  5 11:09:06 PDT 2014
//    Fix reference counting.
//
//    Eric Brugger, Fri Sep 26 08:49:02 PDT 2014
//    I modified the routine to return a NULL in the case where it previously
//    returned an avtDataRepresentation with a NULL vtkDataSet.
//
//    Eric Brugger, Tue Aug 25 10:13:49 PDT 2015
//    I modified the routine to return NULL if the output data set was NULL.
//
// ****************************************************************************

avtDataRepresentation *
avtIsovolumeFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    //
    // Start off by calculating the range of the dataset.
    //
    vtkDataArray *vals = NULL;
    if (in_ds->GetPointData()->GetScalars() != NULL)
        vals = in_ds->GetPointData()->GetScalars();
    else if (in_ds->GetCellData()->GetScalars() != NULL)
        vals = in_ds->GetCellData()->GetScalars();

    if (vals == NULL)
        return in_dr;

    double min, max;
    if(vals->GetDataType() == VTK_FLOAT)
        IsovolumeMinMax(min, max, avtDirectAccessor<float>(vals));
    else if(vals->GetDataType() == VTK_DOUBLE)
        IsovolumeMinMax(min, max, avtDirectAccessor<double>(vals));
    else
        IsovolumeMinMax(min, max, avtTupleAccessor(vals));

    //
    // Check to see if our range is below the min or above the max.  If so,
    // we will have an empty intersection.
    //
    if (max < atts.GetLbound() || min > atts.GetUbound())
    {
        return NULL;
    }

    //
    // Determine if we need to do the min clip or max clip.  Because of the
    // above logic, we can assume that the dataset's max is bigger than
    // the isovolume's lbound and the min is less than the ubound.
    //
    bool doMinClip = false;
    if ((atts.GetLbound() > -1e37) && (min < atts.GetLbound()))
        doMinClip = true;
    bool doMaxClip = false;
    if ((atts.GetUbound() < 1e37) && (max > atts.GetUbound()))
        doMaxClip = true;

    //
    // Do the clipping!
    //
    vtkDataSet *out_ds = in_ds;
    if(doMinClip && doMaxClip)
    {
        vtkDataSet *intermediate = NULL;
        intermediate = ExecuteSingleClip(in_ds, atts.GetLbound(), true);
        if(intermediate->GetNumberOfCells() > 0)
        {
            out_ds = ExecuteSingleClip(intermediate, atts.GetUbound(), false);
            intermediate->Delete();
        }
        else
            out_ds = intermediate;
    }
    else if(doMinClip)
        out_ds = ExecuteSingleClip(in_ds, atts.GetLbound(), true);
    else if(doMaxClip)
        out_ds = ExecuteSingleClip(in_ds, atts.GetUbound(), false);
    bool own = out_ds != in_ds;

    //
    // Make sure there's something there
    //
    if (out_ds->GetNumberOfCells() <= 0)
    {
        // We can't delete the out_ds unless we made it and own it. 
        // Otherwise, we're deleting the reference out from under the in_ds.
        if(own)
            out_ds->Delete();
        return NULL;
    }

    //
    // If we had poly data input, we want poly data output.  The VTK filter
    // only returns unstructured grids, so convert that now.  Note: we don't
    // necessarily have a ugrid, since it might be that we didn't process the
    // dataset.
    //
    if (in_ds->GetDataObjectType() == VTK_POLY_DATA &&
        out_ds->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
    {
        vtkUnstructuredGrid *ugrid = (vtkUnstructuredGrid *) out_ds;
        vtkPolyData *out_pd = vtkPolyData::New();
        out_pd->SetPoints(ugrid->GetPoints());
        out_pd->GetPointData()->ShallowCopy(ugrid->GetPointData());
        out_pd->GetCellData()->ShallowCopy(ugrid->GetCellData());
        vtkIdType ncells = ugrid->GetNumberOfCells();
        out_pd->Allocate(ncells);
        for (vtkIdType i = 0 ; i < ncells ; i++)
        {
            int celltype = ugrid->GetCellType(i);
            vtkIdType *pts, npts;
            ugrid->GetCellPoints(i, npts, pts);
            out_pd->InsertNextCell(celltype, npts, pts);
        }
        if(own)
            out_ds->Delete();
        out_ds = out_pd;
        own = true;
    }

    avtDataRepresentation *out_dr = NULL;
    if (out_ds != NULL)
        out_dr = new avtDataRepresentation(out_ds,
            in_dr->GetDomain(), in_dr->GetLabel());

    if (own && out_ds != NULL)
        out_ds->Delete();

    return out_dr;
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Indicates the zones no longer correspond to the original problem.
//
//  Programmer: Jeremy Meredith
//  Creation:   February 16, 2004
//
//  Modifications:
//    Kathleen Bonnell, Thu Mar  2 14:26:06 PST 2006 
//    Set ZonesSplit.
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtIsovolumeFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().ZonesSplit();

    char params[100];
    SNPRINTF(params, 100, "var=%s, min=%lg, max=%lg", atts.GetVariable().c_str(),
        atts.GetLbound(), atts.GetUbound());
    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("Isovolume", params);
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::ModifyContract
//
//  Purpose:
//      Tell the database that we will need ghost zones.
//
//  Programmer: Hank Childs
//  Creation:   August 11, 2004
//
//  Modifications:
//
//    Hank Childs, Wed Oct 20 17:04:52 PDT 2004
//    Use interval trees to only read in the domains we need.
//
//    Kathleen Bonnell, Mon Aug 14 16:40:30 PDT 2006
//    API change for avtIntervalTree.
//
// ****************************************************************************

avtContract_p
avtIsovolumeFilter::ModifyContract(avtContract_p in_spec)
{
    avtContract_p spec = new avtContract(in_spec);

    const char *varname = NULL;
    if (atts.GetVariable() != "default")
        varname = atts.GetVariable().c_str();
    else
        varname = in_spec->GetDataRequest()->GetVariable();

    //
    // We will need the ghost zones so that we can interpolate along domain
    // boundaries and get no cracks in our isosurface.
    //
    avtDataAttributes &in_atts = GetInput()->GetInfo().GetAttributes();
    bool skipGhost = false;
    if (in_atts.ValidVariable(varname) &&
        in_atts.GetCentering(varname) == AVT_NODECENT)
        skipGhost = true;
    if (!skipGhost)
        spec->GetDataRequest()->SetDesiredGhostDataType(GHOST_ZONE_DATA);

    std::string iso_var(atts.GetVariable());
    if (iso_var == "default")
        iso_var = in_spec->GetDataRequest()->GetVariable();

    avtIntervalTree *it = GetMetaData()->GetDataExtents(iso_var.c_str());
    if (it != NULL)
    {
        double min = atts.GetLbound();
        double max = atts.GetUbound();
        std::vector<int> dl;
        it->GetElementsListFromRange(&min, &max, dl);
        spec->GetDataRequest()->GetRestriction()->RestrictDomains(dl);
    }

    return spec;
}


