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

#define ENABLE_UNGHOST_CELL_OPTIMIZATION

#include <avtAMRStitchCellFilter.h>
#include "AMRStitchCellTesselations2D.h"
#include "AMRStitchCellTesselations3D.h"

#include <algorithm>
#include <cassert>
#include <limits>
#include <list>
#include <string>
#include <vector>

#include <DebugStream.h>

#include <avtDataTree.h>
#include <avtGhostData.h>
#include <avtIntervalTree.h>
#include <avtMetaData.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtStructuredDomainNesting.h>

#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkFieldData.h>
#include <vtkFloatArray.h>
#include <vtkIdList.h>
#include <vtkIntArray.h>
#include <vtkObjectFactory.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>

#include <vtkDataSetWriter.h>
#include <vtkVisItUtility.h>

// ****************************************************************************
//  Method: avtAMRStitchCellFilter constructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
// ****************************************************************************

avtAMRStitchCellFilter::avtAMRStitchCellFilter()
{
}


// ****************************************************************************
//  Method: avtAMRStitchCellFilter destructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
//  Modifications:
//
// ****************************************************************************

avtAMRStitchCellFilter::~avtAMRStitchCellFilter()
{
}


// ****************************************************************************
//  Method:  avtAMRStitchCellFilter::Create
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
// ****************************************************************************

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


// ****************************************************************************
//  Method:      avtAMRStitchCellFilter::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 Jul 8 15:14:01 PST 2010
//
// ****************************************************************************

void
avtAMRStitchCellFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const AMRStitchCellAttributes*)a;
}


// ****************************************************************************
//  Method: avtAMRStitchCellFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtAMRStitchCellFilter with the given
//      parameters would result in an equivalent avtAMRStitchCellFilter.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
// ****************************************************************************

bool
avtAMRStitchCellFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(AMRStitchCellAttributes*)a);
}


// ****************************************************************************
//  Method: avtAMRStitchCellFilter::ModifyContract
//
//  Purpose:
//      Request creation of ghost zones
//
//  Arguments:
//      in_contract  The input contract.
//
//  Programmer: Hank Childs
//  Creation:   June 6, 2001
//
//  Modifications:
//
// ****************************************************************************

avtContract_p
avtAMRStitchCellFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p contract = new avtContract(in_contract);
    // FIXME: Ensure that data is node centered
    contract->GetDataRequest()->SetDesiredGhostDataType(GHOST_ZONE_DATA);
    return contract;
}

// ****************************************************************************
//  Method: avtAMRStitchCellFilter::PreExecute
//
//  Purpose:
//      Compute extents to obtain origin of all boxes/patches.
//
//  Programmer: Gunther H. Weber
//  Creation:   August 5, 2010
//
//  Modifications:
//
//    Gunther H. Weber, Wed Jul 18 17:29:30 PDT 2012
//    Fixed memory error reported by valgrind
//
// ****************************************************************************

void
avtAMRStitchCellFilter::PreExecute(void)
{
    // Get dimension of data set
    avtDataAttributes &inAtts = GetInput()->GetInfo().GetAttributes();
    topologicalDimension = inAtts.GetTopologicalDimension();
    spatialDimension = inAtts.GetSpatialDimension();

    // Check dimension
    if (topologicalDimension < 2 || topologicalDimension > 3)
        EXCEPTION1(ImproperUseException,
                "Need 2D or 3D data set to generate dual mesh and stitch cells.");
    if (spatialDimension < 2 || spatialDimension > 3)
        EXCEPTION1(ImproperUseException,
                "Need 2D or 3D data set to generate dual mesh and stitch cells.");
    if (spatialDimension != topologicalDimension)
        EXCEPTION1(ImproperUseException,
                "Topological and spatial dimension must be identical to generate "
                "dual mesh and stitch cells.");

    // Obtian data set origin via spatial extents
    avtIntervalTree *iTree = GetMetaData()->GetSpatialExtents();
    if (!iTree)
        EXCEPTION1(ImproperUseException,
                "Cannot determine spatial extents of data set.");
    double domainBoundingBox[6];
    iTree->GetExtents(domainBoundingBox);

    // Set domainOrigin
    for (int i=0; i<spatialDimension; ++i)
        domainOrigin[i] = domainBoundingBox[2*i];
    for (int i=spatialDimension; i<3; ++i)
        domainOrigin[i] = 0.0;

    // We need structured domain nesting to figure out extents of boxes/patches,
    // levels and refinement ratio wrt. to the parent level
    domainNesting =
        dynamic_cast<avtStructuredDomainNesting*>(GetMetaData()->GetDomainNesting());
    if (!domainNesting)
        EXCEPTION1(ImproperUseException,
                "Need structured domain nesting information to generate dual "
                "mesh and stitch cells.");

    // Get logical extents of root level
    size_t nLevels = domainNesting->GetNumberOfLevels();
    if (nLevels == 0)
         EXCEPTION1(ImproperUseException,
                 "Data sets does not contain any levels (according to information "
                 "in structured domain nesting.");

    logicalDomainBoundingBox.resize(nLevels);
    logicalDomainBoundingBox[0].resize(6);
    logicalDomainBoundingBox[0][0] = std::numeric_limits<int>::max();
    logicalDomainBoundingBox[0][1] = std::numeric_limits<int>::max();
    logicalDomainBoundingBox[0][2] = std::numeric_limits<int>::max();
    logicalDomainBoundingBox[0][3] = std::numeric_limits<int>::min();
    logicalDomainBoundingBox[0][4] = std::numeric_limits<int>::min();
    logicalDomainBoundingBox[0][5] = std::numeric_limits<int>::min();
 
    std::vector<int> le;
    for (int dom = 0; dom < (int)domainNesting->GetNumberOfDomains(); ++dom)
    {
        if (domainNesting->GetDomainLevel(dom) == 0)
        {
            le = domainNesting->GetDomainLogicalExtents(dom);
            logicalDomainBoundingBox[0][0] = std::min(logicalDomainBoundingBox[0][0], le[0]);
            logicalDomainBoundingBox[0][1] = std::min(logicalDomainBoundingBox[0][1], le[1]);
            logicalDomainBoundingBox[0][2] = std::min(logicalDomainBoundingBox[0][2], le[2]);
            logicalDomainBoundingBox[0][3] = std::max(logicalDomainBoundingBox[0][3], le[3]);
            logicalDomainBoundingBox[0][4] = std::max(logicalDomainBoundingBox[0][4], le[4]);
            logicalDomainBoundingBox[0][5] = std::max(logicalDomainBoundingBox[0][5], le[5]);
        }
    }

    debug5 << "avtAMRStitchCellFilter::PreExecute(): logicalDomainBoundingBox[0][6] = { ";
    debug5 << logicalDomainBoundingBox[0][0] << ", " << logicalDomainBoundingBox[0][1] << ", ";
    debug5 << logicalDomainBoundingBox[0][2] << ", " << logicalDomainBoundingBox[0][3] << ", ";
    debug5 << logicalDomainBoundingBox[0][4] << ", " << logicalDomainBoundingBox[0][5] << " }";
    debug5 << std::endl;

    // Compute logical bounding box in index space of each level
    for (size_t l=1; l<nLevels; ++l)
    {
        const std::vector<int>& refRatio = domainNesting->GetLevelRefinementRatios((int)l);
        if (refRatio.size() != (size_t)topologicalDimension)
            EXCEPTION1(ImproperUseException,
                    "Refinement ratio provided by database via domain nesting is invalid. "
                    "Expected a vector of length equal to topological dataset dimension.");

        logicalDomainBoundingBox[l].resize(6);
        logicalDomainBoundingBox[l][0] = refRatio[0] * logicalDomainBoundingBox[l-1][0];
        logicalDomainBoundingBox[l][1] = refRatio[1] * logicalDomainBoundingBox[l-1][1];
        logicalDomainBoundingBox[l][3] = refRatio[0] * (logicalDomainBoundingBox[l-1][3] + 1) - 1;
        logicalDomainBoundingBox[l][4] = refRatio[1] * (logicalDomainBoundingBox[l-1][4] + 1) - 1;
        if (topologicalDimension == 2)
        {
            logicalDomainBoundingBox[l][2] = 0;
            logicalDomainBoundingBox[l][5] = 0;
        }
        else
        {
            logicalDomainBoundingBox[l][2] = refRatio[2] * logicalDomainBoundingBox[l-1][2];
            logicalDomainBoundingBox[l][5] = refRatio[2] * (logicalDomainBoundingBox[l-1][5] + 1) - 1;
        }

        debug5 << "avtAMRStitchCellFilter::PreExecute(): logicalDomainBoundingBox[" << l << "][6] = { ";
        debug5 << logicalDomainBoundingBox[l][0] << ", " << logicalDomainBoundingBox[l][1] << ", ";
        debug5 << logicalDomainBoundingBox[l][2] << ", " << logicalDomainBoundingBox[l][3] << ", ";
        debug5 << logicalDomainBoundingBox[l][4] << ", " << logicalDomainBoundingBox[l][5];
        debug5 << " } " << std::endl;
    }

    // Get cell sizes for levels
    cellSize.resize(nLevels);
    for (size_t l=0; l<nLevels; ++l)
    {
        cellSize[l] = domainNesting->GetLevelCellSizes((int)l);
        if (cellSize[l].size() != (size_t)spatialDimension)
        {
            EXCEPTION1(ImproperUseException,
                    "Database plugin did not properly set level cell sizes using "
                    "SetLevelCellSizes() method.");
        }
    }
}

// ****************************************************************************
//  Method: avtAMRStitchCellFilter::ExecuteDataTree
//
//  Purpose:
//      Create dual grid and stitch cell data sets and add them to a data tree
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: Gunther H. Weber
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
//  Modifications:
//    Eric Brugger, Tue Jul 22 17:02:41 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataTree_p
avtAMRStitchCellFilter::ExecuteDataTree(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set, the domain number and the label.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();
    int domain = in_dr->GetDomain();
    std::string label = in_dr->GetLabel();

    // Diagnostic output
    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): Processing domain ";
    debug5 << domain << " (" << label << ")" << std::endl;

    // Obtain level
    const int level = domainNesting->GetDomainLevel(domain);
    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): Patch is in level ";
    debug5 << level << std::endl;

    // Ensure that we are working on rectilinear grid
    vtkRectilinearGrid *rgrid = dynamic_cast<vtkRectilinearGrid*>(in_ds);
    if (!rgrid)
        EXCEPTION1(ImproperUseException,
                "Can only create dual mesh and stitch cells for a rectilinear grid.");

    // Dimensions of data set
    int dims[3];
    rgrid->GetDimensions(dims);

    // We want number of cells as grid dimension, not number of samples
    for (int d=0; d<topologicalDimension; ++d)
        dims[d]--;

    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): Grid is regular and has ";
    debug5 << "extents " << dims[0] << " " << dims[1] << " " << dims[2] << std::endl;

    // Get base_index (measured in number of cells of current level, i.e.,
    // same level as patch, form origin)
    vtkIntArray *baseIdxArray =
        (vtkIntArray*) rgrid->GetFieldData()->GetArray("base_index");
    if (!baseIdxArray)
        EXCEPTION1(ImproperUseException,
                "Need base_index data to generate dual mesh and stitch cells.");
    int baseIdx[3];
    for (int d = 0; d < topologicalDimension; ++d)
        baseIdx[d] = baseIdxArray->GetValue(d);
    for (int d = topologicalDimension; d < 3; ++d)
        baseIdx[d] = 0;

    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): base_index is ";
    debug5 << baseIdx[0] << ", " << baseIdx[1] << ", " << baseIdx[2] << std::endl;

    // Consistency check: layer of ghost cells needs to have a width of one

    // ... Get real dims information
    vtkIntArray *realDimsArray =
        dynamic_cast<vtkIntArray*>(rgrid->GetFieldData()->GetArray("avtRealDims"));
    if (!realDimsArray)
        EXCEPTION1(ImproperUseException, "Need avtRealDims to generate stitch cells.");

    int realIMin = realDimsArray->GetValue(0);
    int realIMax = realDimsArray->GetValue(1);
    int realJMin = realDimsArray->GetValue(2);
    int realJMax = realDimsArray->GetValue(3);
    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): Real dims are ";
    debug5 << realIMin << " " << realIMax << " " << realJMin << " " << realJMax << std::endl;

    if (realIMin > 1 || (baseIdx[0] != logicalDomainBoundingBox[level][0] && realIMin != 1) ||
        (dims[0] - realIMax) > 1 || ((baseIdx[0] + dims[0] - 1 - realIMin) != logicalDomainBoundingBox[level][3] && (dims[0] - realIMax) != 1) ||
        realJMin > 1 || (baseIdx[1] != logicalDomainBoundingBox[level][1] && realJMin != 1) ||
        (dims[1] - realJMax) > 1 || ((baseIdx[1] + dims[1] - 1 - realJMin) != logicalDomainBoundingBox[level][4] && (dims[1] - realJMax) != 1))
    {
        EXCEPTION1(ImproperUseException,
                "Need exactly one layer of ghost cells (except at domain boundary) to create stitch cells.");
    }

    // Offset for accessing data
    int ghostOffset[3] = { realIMin, realJMin, 0 };

    // Real size of data set
    int realDim[3] = {
        dims[0] - realIMin - (dims[0] - realIMax),
        dims[1] - realJMin - (dims[1] - realJMax),
        0
    };

    if (topologicalDimension == 3)
    {
        int realKMin = realDimsArray->GetValue(4);
        int realKMax = realDimsArray->GetValue(5);

        if (realKMin > 1 || (baseIdx[2] != logicalDomainBoundingBox[level][2] && realKMin != 1) ||
            (dims[2] - realKMax) > 1 || ((baseIdx[2] + dims[2] - 1 - realKMin) != logicalDomainBoundingBox[level][5] && (dims[2] - realKMax) != 1)) 
        {
            EXCEPTION1(ImproperUseException,
                    "Need exactly one layer of ghost cells (except at domain boundary) to create stitch cells.");
        }

        ghostOffset[2] = realKMin;
        realDim[2] = dims[2] - realKMin - (dims[2] - realKMax);
    }

    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): ghostOffset[3] = { "
        << ghostOffset[0] << ", " << ghostOffset[1] << ", " << ghostOffset[2]
        << " } " << std::endl;
    debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): realDim[3] = { "
        << realDim[0] << ", " << realDim[1] << ", " << realDim[2]
        << " } " << std::endl;

    // Initialize refineInDameLevelId array -> For all ghost cells surrounding the grid,
    // this array contains -1 if that cell does not have a neighboring grid in the same
    // level. If there is a neighboring grid in the same level, this array contains the 
    // domain number of that grid.
    vtkIdType *refinedInSameLevelDomain;

    avtStructuredDomainBoundaries *db = dynamic_cast<avtStructuredDomainBoundaries*>(
            GetMetaData()->GetDomainBoundaries());
    if (!db)
        EXCEPTION1(ImproperUseException,
                "Need structured domain boundaries to compute stitch cells.");
    std::vector<Neighbor> neighbors = db->GetNeighbors(domain);

    if (topologicalDimension == 2)
    {
        refinedInSameLevelDomain = new vtkIdType[dims[0]*dims[1]]; // FIXME: More memory efficient storage
        for (int i=0; i < dims[0]*dims[1]; ++i)
        {
            // Fill with pre-set value of -1 to mark cells without neighboring grid
            refinedInSameLevelDomain[i] = -1;
        }

        std::vector<int> le;
        debug5 << "Patch " << domain << " has " << neighbors.size() << " neighbors." << std::endl;
        for (std::vector<Neighbor>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
        {
            if (it->refinement_rel != SAME_REFINEMENT_LEVEL)
            {
                debug5 << "Skipping neighbor " << it->domain << " with different refinement level." << std::endl;
                continue;
            }

            switch(it->type)
            {
                case Boundary::IMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]); j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        refinedInSameLevelDomain[j*dims[0]] = it->domain;
                    break;
                case Boundary::IMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]); j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        refinedInSameLevelDomain[j*dims[0]+dims[0]-1] = it->domain;
                    break;
                case Boundary::JMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]); i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        refinedInSameLevelDomain[i] = it->domain;
                    break;
                case Boundary::JMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]); i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        refinedInSameLevelDomain[(dims[1]-1)*dims[0]+i] = it->domain;
                    break;
                case Boundary::IMIN | Boundary::JMIN:
                    refinedInSameLevelDomain[0] = it->domain;
                    break;
                case Boundary::IMIN | Boundary::JMAX:
                    refinedInSameLevelDomain[(dims[1]-1)*dims[0]] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMIN:
                    refinedInSameLevelDomain[dims[0]-1] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMAX:
                    refinedInSameLevelDomain[(dims[1]-1)*dims[0]+dims[0]-1] = it->domain;
                    break;
                default:
                    debug5 << "Encountered invalid boaundary: ";
                    debug5 << (it->type & Boundary::IMIN ? "(IMIN)" : "");
                    debug5 << (it->type & Boundary::IMAX ? "(IMAX)" : "");
                    debug5 << (it->type & Boundary::JMIN ? "(JMIN)" : "");
                    debug5 << (it->type & Boundary::JMAX ? "(JMAX)" : "");
                    EXCEPTION1(ImproperUseException, "Encountered invalid boundary.");
            }
        }
    }
    else
    {
        refinedInSameLevelDomain = new vtkIdType[dims[0]*dims[1]*dims[2]]; // FIXME: More memory efficient storage
        for (int i=0; i < dims[0]*dims[1]*dims[2]; ++i)
        {
            // Fill with pre-set value of -1 to mark cells without neighboring grid
            refinedInSameLevelDomain[i] = -1;
        }

#ifdef LOC
#undef LOC
#endif
#define LOC(i,j,k) (((k)*dims[1])+j)*dims[0]+i
        std::vector<int> le;
        for (std::vector<Neighbor>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
        {
            if (it->refinement_rel != SAME_REFINEMENT_LEVEL)
            {
                debug5 << "Skipping neighbor " << it->domain << " with different refinement level." << std::endl;
                continue;
            }

            switch(it->type)
            {
                case Boundary::IMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        for(int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                            k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                        {
                            refinedInSameLevelDomain[LOC(0, j, k)] = it->domain;
                        }
                    break;
                case Boundary::IMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        for (int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                             k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                        {
                            refinedInSameLevelDomain[LOC(dims[0]-1, j, k)] = it->domain;
                        }
                    break;
                case Boundary::JMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        for (int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                             k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                        {
                            refinedInSameLevelDomain[LOC(i, 0, k)] = it->domain;
                        }
                    break;
                case Boundary::JMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        for (int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                             k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                        {
                            refinedInSameLevelDomain[LOC(i, dims[1]-1, k)] = it->domain;
                        }
                    break;
                case Boundary::KMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                             j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        {
                            refinedInSameLevelDomain[LOC(i, j, 0)] = it->domain;
                        }
                    break;
                case Boundary::KMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                        for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                             j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                        {
                            refinedInSameLevelDomain[LOC(i, j, dims[2]-1)] = it->domain;
                        }
                    break;
                case Boundary::IMIN | Boundary::JMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for(int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                        k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                    {
                        refinedInSameLevelDomain[LOC(0, 0, k)] = it->domain;
                    }
                    break;
                case Boundary::IMIN | Boundary::JMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for(int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                        k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                    {
                        refinedInSameLevelDomain[LOC(0, dims[1]-1, k)] = it->domain;
                    }
                    break;
                case Boundary::IMAX | Boundary::JMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for(int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                        k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                    {
                        refinedInSameLevelDomain[LOC(dims[0]-1, 0, k)] = it->domain;
                    }
                    break;
                case Boundary::IMAX | Boundary::JMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for(int k=std::max(0, le[2]-baseIdx[2]+ghostOffset[2]);
                        k<std::min(dims[2], le[5]-baseIdx[2]+ghostOffset[2]+1); ++k)
                    {
                        refinedInSameLevelDomain[LOC(dims[0]-1, dims[1]-1, k)] = it->domain;
                    }
                    break;
                case Boundary::IMIN | Boundary::KMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                    {
                        refinedInSameLevelDomain[LOC(0, j, 0)] = it->domain;
                    }
                    break;
                case Boundary::IMIN | Boundary::KMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                    {
                        refinedInSameLevelDomain[LOC(0, j, dims[2]-1)] = it->domain;
                    }
                    break;
                case Boundary::IMAX | Boundary::KMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                    {
                        refinedInSameLevelDomain[LOC(dims[0]-1, j, 0)] = it->domain;
                    }
                    break;
                case Boundary::IMAX | Boundary::KMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int j=std::max(0, le[1]-baseIdx[1]+ghostOffset[1]);
                         j<std::min(dims[1], le[4]-baseIdx[1]+ghostOffset[1]+1); ++j)
                    {
                        refinedInSameLevelDomain[LOC(dims[0]-1, j, dims[2]-1)] = it->domain;
                    }
                    break;
                case Boundary::JMIN | Boundary::KMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                    {
                        refinedInSameLevelDomain[LOC(i, 0, 0)] = it->domain;
                    }
                    break;
                case Boundary::JMIN | Boundary::KMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                    {
                        refinedInSameLevelDomain[LOC(i, 0, dims[2]-1)] = it->domain;
                    }
                    break;
                case Boundary::JMAX | Boundary::KMIN:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                    {
                        refinedInSameLevelDomain[LOC(i, dims[1]-1, 0)] = it->domain;
                    }
                    break;
                case Boundary::JMAX | Boundary::KMAX:
                    le = domainNesting->GetDomainLogicalExtents(it->domain);
                    for (int i=std::max(0, le[0]-baseIdx[0]+ghostOffset[0]);
                         i<std::min(dims[0], le[3]-baseIdx[0]+ghostOffset[0]+1); ++i)
                    {
                        refinedInSameLevelDomain[LOC(i, dims[1]-1, dims[2]-1)] = it->domain;
                    }
                    break;
                case Boundary::IMIN | Boundary::JMIN | Boundary::KMIN:
                    refinedInSameLevelDomain[LOC(0, 0, 0)] = it->domain;
                    break;
                case Boundary::IMIN | Boundary::JMIN | Boundary::KMAX:
                    refinedInSameLevelDomain[LOC(0, 0, dims[2]-1)] = it->domain;
                    break;
                case Boundary::IMIN | Boundary::JMAX | Boundary::KMIN:
                    refinedInSameLevelDomain[LOC(0, dims[1]-1, 0)] = it->domain;
                    break;
                case Boundary::IMIN | Boundary::JMAX | Boundary::KMAX:
                    refinedInSameLevelDomain[LOC(0, dims[1]-1, dims[2]-1)] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMIN | Boundary::KMIN:
                    refinedInSameLevelDomain[LOC(dims[0]-1, 0, 0)] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMIN | Boundary::KMAX:
                    refinedInSameLevelDomain[LOC(dims[0]-1, 0, dims[2]-1)] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMAX | Boundary::KMIN:
                    refinedInSameLevelDomain[LOC(dims[0]-1, dims[1]-1, 0)] = it->domain;
                    break;
                case Boundary::IMAX | Boundary::JMAX | Boundary::KMAX:
                    refinedInSameLevelDomain[LOC(dims[0]-1, dims[1]-1, dims[2]-1)] = it->domain;
                    break;
                default:
                    EXCEPTION1(ImproperUseException, "Encountered invalid boundary.");
            }
        }

#undef LOC
    }

    // Output data sets: Dual grid and stitch cells
    vtkDataSet *out_ds[2];

    if (atts.GetCreateCellsOfType() == AMRStitchCellAttributes::DualGridAndStitchCells ||
        atts.GetCreateCellsOfType() == AMRStitchCellAttributes::DualGrid)
        out_ds[0] = CreateDualGrid(rgrid, domain, level, dims, baseIdx, ghostOffset, refinedInSameLevelDomain);
    else
        out_ds[0] = 0;

    if ((atts.GetCreateCellsOfType() == AMRStitchCellAttributes::DualGridAndStitchCells ||
         atts.GetCreateCellsOfType() == AMRStitchCellAttributes::StitchCells) &&
        level != 0)
    {
        // We are not in the root level and thus need to generate stich cells.
        // Obtain refinement ratio wrt. to parent level needed to translate indices
        // between levels.
        std::vector<int> refinementRatio = domainNesting->GetRatiosForLevel(level-1, domain);
        debug5 << "avtAMRStitchCellFilter::ExecuteDataTree(): Refinement ratio ";
        debug5 << " wrt. parent level is " << refinementRatio[0] << " ";
        debug5 << refinementRatio[1] << " " << refinementRatio[2] << std::endl;

        // Consitency check, patch must start at the boundaries of a coarse level cell
        if ((baseIdx[0] % refinementRatio[0]) || (baseIdx[1] % refinementRatio[1]) ||
            ((topologicalDimension == 3) && (baseIdx[2] % refinementRatio[2])))
            EXCEPTION1(ImproperUseException,
                    "base_index must be divisible by refinement_ratio, i.e. "
                    "fine box/patch must start at coarse cell boundaries.");

        out_ds[1] = CreateStitchCells(rgrid, domain, level, dims, realDim, baseIdx,
               ghostOffset, refinementRatio, refinedInSameLevelDomain);
    }
    else
    {
        // We are in the root level. There is no lower level to connect to.
        // Connection to boxes/patches in the same level is handled by CreateDualGrid
        out_ds[1] = 0;
        debug5 << "Skipping stitch cell generation for patch in root level.";
        debug5 << std::endl;
    }

    // Clean-up no longer needed information
    delete[] refinedInSameLevelDomain;

    // Create data tree
    avtDataTree_p rv = new avtDataTree(2, out_ds, domain, label);

    // Clean-up intermediate results
    for (int dsNo=0; dsNo<2; ++dsNo)
        if (out_ds[dsNo]) out_ds[dsNo]->Delete();

    // Return result
    return rv;
}

// ****************************************************************************
//  Method: avtAMRStitchCellFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Indicates that zone information has changed.
//
//  Programmer: Gunther H. Weber
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Wed Jul 10 14:54:09 PDT 2013
//    Add updating centering
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtAMRStitchCellFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateDataMetaData();
    avtDataAttributes &out_atts = GetOutput()->GetInfo().GetAttributes();
    // Change centering for all variables from zone to node centered
    for (int varNo = 0; varNo < out_atts.GetNumberOfVariables(); ++varNo)
    {
        std::string varname = out_atts.GetVariableName(varNo);
        if (out_atts.GetCentering(varname.c_str()) == AVT_ZONECENT)
            out_atts.SetCentering(AVT_NODECENT, varname.c_str());
    }
    //GetOutput()->GetInfo().GetAttributes().SetCentering(AVT_NODECENT);

    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("AMRStitchCell");
}
 
// ****************************************************************************
//  Method: avtAMRStitchCellFilter::CreateStitchCells
//
//  Purpose:
//      Create stitch cells connecting the dual grid to its surroundings
//
//  Arguments:
//      in_ds      The input dataset.
//      domain     The domain number.
//
//  Returns:       The output dataset.
//
//  Programmer: Gunther H. Weber
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
//  Modifications:
//    Kathleen Biagas, Tue Oct 16 15:28:09 MST 2012
//    Preserve Coordinate type.
//
// ****************************************************************************

vtkDataSet *
avtAMRStitchCellFilter::CreateStitchCells(vtkRectilinearGrid *rgrid,
        int domain, int level, int dims[3], int realDim[3], int baseIdx[3],
        int ghostOffset[3], const std::vector<int> &refinementRatio,
        const vtkIdType* refinedInSameLevelDomain)
{
    // Create unstructured grid for stitch cells
    vtkPoints *stitchCellPts = vtkVisItUtility::NewPoints(rgrid);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints(stitchCellPts);
    stitchCellPts->Delete();

    // Get pointer to ghost information for input data
    vtkUnsignedCharArray *inputGhostZones =
        (vtkUnsignedCharArray *) rgrid->GetCellData()->GetArray("avtGhostZones");
 
    // We store a list of IDs for data attributes we need to copy from
    // the original mesh to the stitch mesh.
    std::list<vtkIdType> originalDataIds;
    std::list<vtkIdList*> interpolatedDataIds;

    vtkIdType *pointIds = new vtkIdType[dims[0]*dims[1]*dims[2]]; // FIXME: More memory efficient storage
    for (int i=0; i < dims[0]*dims[1]*dims[2]; ++i)
    {
        // Fill with specified invalid value for debugging purposes
        pointIds[i] = -2;
    }

    // Needed to check whether coarse vertex is out of bounds
    int maxGlobalI = logicalDomainBoundingBox[level][3];
    int maxGlobalJ = logicalDomainBoundingBox[level][4];

    if (topologicalDimension == 2)
    {
#ifdef LOC
#undef LOC
#endif
#define LOC(i,j) ((j+ghostOffset[1])*dims[0]+(i+ghostOffset[0]))
        // Create stitch cells
        for (int j=-1; j<realDim[1]; ++j)
        {
            if (baseIdx[1]+j<0 || baseIdx[1]+j+1 > maxGlobalJ)
                continue; // Skip cells connecting to outside the domain

            for (int i=-1; i<realDim[0]; i+=((j==-1 || j==realDim[1]-1) ? 1 : realDim[0]))
            {
                if (baseIdx[0]+i<0 || baseIdx[0]+i+1 > maxGlobalI)
                    continue; // Skip cells connecting to outside the domain

                // Indices of the four vertices defining a potential stitch cell
                int vIdx[4][2] = {
                    { i, j },
                    { i, j + 1 },
                    { i + 1, j + 1},
                    { i + 1, j }
                };

                // First, compute case number and check whether cell can be added.
                // A cell cannot be added if (i) it would extent outside the domain,
                // (ii) coincides partially or entirely with a refining grid, or (iii)
                // connects to an adjacent grid in the same level with a smaller domain
                // number (to avoid duplicate cells, each stitch cell belongs to the 
                // patch in a level with the smallest number).
                bool skipCell = false;
                int caseNo = 0;

                for (int vtxNo=0; vtxNo < 4; ++vtxNo)
                {
                    if (// Vertex is in interior, completely refined region:
                            ((vIdx[vtxNo][0] >= 0) && (vIdx[vtxNo][0] < realDim[0]) &&
                             (vIdx[vtxNo][1] >= 0) && (vIdx[vtxNo][1] < realDim[1])) || 
                            // Vertex belongs to a neighboring refining grid:
                            refinedInSameLevelDomain[LOC(vIdx[vtxNo][0], vIdx[vtxNo][1])] != -1)
                    {
                        // Vertex is in same level as the patch for which we are
                        // creating stitch cells.

                        // Set bit in case number indicating that this point is refined
                        caseNo |= 1 << vtxNo;

                        // If cell is in same level, it cannot be outside the domain
                        // since it belongs to a fine patch

                        // Check whether it coincides with a fine grid
                        if (avtGhostData::IsGhostZoneType(
                                    inputGhostZones->GetValue(
                                        LOC(vIdx[vtxNo][0], vIdx[vtxNo][1])),
                                    REFINED_ZONE_IN_AMR_GRID))
                        {
                            skipCell = true;
                            break;
                        }

                        // Check whether it connects to a patch in the same level
                        // with a smaller domain number
                        if (// Has to lie on boundary to connect to another patch
                                ((vIdx[vtxNo][0] == -1) || (vIdx[vtxNo][0] == realDim[0]) ||
                                 (vIdx[vtxNo][1] == -1) || (vIdx[vtxNo][1] == realDim[1])) && 
                                // And have an domain number smaller than us
                                refinedInSameLevelDomain[LOC(vIdx[vtxNo][0], vIdx[vtxNo][1])] < domain)
                        {
                            skipCell = true;
                            break;
                        }
                    }
                    else
                    {
                        // Vertex is in coarse level
                        int globalI = baseIdx[0] + vIdx[vtxNo][0];
                        int globalJ = baseIdx[1] + vIdx[vtxNo][1];

                        if(globalI<0 || globalI > maxGlobalI ||
                           globalJ<0 || globalJ > maxGlobalJ)
                            EXCEPTION1(VisItException, "Failed assert (internal error).");
                        /*
                        if (globalI < 0 || globalI > maxGlobalI ||
                                globalJ < 0 || globalJ > maxGlobalJ)
                        {
                            skipCell = true;
                            break;
                        }
                        */
                    }
                }

                if (!skipCell)
                {
                    // Compute point locations and store them
                    vtkIdType vPtId[4];
                    for (int vtxNo = 0; vtxNo < 4; ++ vtxNo)
                    {
                        int vI = vIdx[vtxNo][0];
                        int vJ = vIdx[vtxNo][1];

                        if (pointIds[LOC(vI, vJ)] < 0)
                        {
                            // Vertex was not computed for previous stitch cell
                            if (caseNo & (1 << vtxNo))
                            {
                                // Fine cell
                                pointIds[LOC(vI, vJ)] = stitchCellPts->InsertNextPoint(
                                        domainOrigin[0] + (baseIdx[0] + vI  + 0.5) * cellSize[level][0],
                                        domainOrigin[1] + (baseIdx[1] + vJ  + 0.5) * cellSize[level][1],
                                        0
                                        );
                                // FIXME: For mapped grids (not supported) topological and spatial dimension
                                // could differ. If implemented, we need to check spatial dimension here.
                                originalDataIds.push_back(vtkIdType(LOC(vI,vJ)));
                            }
                            else
                            {
                                // Vertex i and j in parent level indices
                                int pI = (baseIdx[0] + vI) / refinementRatio[0];
                                int pJ = (baseIdx[1] + vJ) / refinementRatio[1];

                                pointIds[LOC(vI, vJ)] = stitchCellPts->InsertNextPoint(
                                        domainOrigin[0] + (pI + 0.5) * cellSize[level-1][0],
                                        domainOrigin[1] + (pJ + 0.5) * cellSize[level-1][1],
                                        0
                                        );
                                // FIXME: For mapped grids (not supported) topological and spatial dimension
                                // could differ. If implemented, we need to check spatial dimension here.
                                originalDataIds.push_back(vtkIdType(LOC(vI, vJ)));

                                if ((vI == -1) || (vI == realDim[0]))
                                {
                                    if ((vJ != -1) && (vJ != realDim[1]))
                                    {
                                        int baseJ = vJ - (vJ % refinementRatio[1]);
                                        for (int refVtxNo = 0; refVtxNo < refinementRatio[1];
                                                ++refVtxNo)
                                        {
                                            pointIds[LOC(vI,baseJ+refVtxNo)]=pointIds[LOC(vI,vJ)];
                                        }
                                    }
                                }
                                else
                                {
                                    assert ((vJ == -1) || (vJ == realDim[1]));

                                    if ((vI != -1) && (vI != realDim[0]))
                                    {
                                        int baseI = vI - (vI % refinementRatio[0]);
                                        for (int refVtxNo = 0; refVtxNo < refinementRatio[0];
                                                ++refVtxNo)
                                        {
                                            pointIds[LOC(baseI+refVtxNo,vJ)]=pointIds[LOC(vI,vJ)];
                                        }
                                    }
                                }
                            }
                        }

                        // Copy previously or newly computed point
                        vPtId[vtxNo] = pointIds[LOC(vI, vJ)];
                    }

#ifdef ENABLE_UNGHOST_CELL_OPTIMIZATION
                    if (caseNo != 15)
#endif
                    {
                        unsigned char subCaseNo = 0;
                        if (tesselationSubCaseDir2D[caseNo] == 1)
                        {
                            subCaseNo = ((i + 1) % refinementRatio[0] == 0) ? 0 : 1;
                        }
                        else if (tesselationSubCaseDir2D[caseNo] == 2)
                        {
                            subCaseNo = ((j + 1) % refinementRatio[1] == 0) ? 0 : 1;
                        }

                        int tessStart = tesselationCaseStart2D[caseNo][subCaseNo];

                        if (tessStart == -1)
                        {
                            EXCEPTION1(VisItException, "Invalid tesselation case (internal error).");
                        }
                        else
                        {
                            if (tesselationArray2D[tessStart] == 4)
                            {
                                ugrid->InsertNextCell(VTK_QUAD, 4, vPtId);
                            }
                            else if (tesselationArray2D[tessStart] == 3)
                            {
                                vtkIdType cPtId[3];
                                for (int i=0; i<3; ++i)
                                    cPtId[i]=vPtId[tesselationArray2D[++tessStart]];
                                ugrid->InsertNextCell(VTK_TRIANGLE, 3, cPtId);
                            }
                            else
                                EXCEPTION1(VisItException, "Invalid cell type (internal error).");
                        }
                    }
                }
            }
        }
    }
    else
    {
        // Create 3D stitch cells
        int maxGlobalK = logicalDomainBoundingBox[level][5];
#ifdef LOC
#undef LOC
#endif
#define LOC(i,j,k) (((k+ghostOffset[2])*dims[1]+j+ghostOffset[1])*dims[0]+i+ghostOffset[0])
        // Create stitch cells
        for (int k=-1; k<realDim[2]; ++k)
        {
            if (baseIdx[2]+k<0 || baseIdx[2]+k+1 > maxGlobalK) continue; // Skip cells connecting to outside the domain

            for (int j=-1; j<realDim[1]; ++j)
            {
                if (baseIdx[1]+j<0 || baseIdx[1]+j+1 > maxGlobalJ) continue; // Skip cells connecting to outside the domain

                for (int i=-1; i<realDim[0]; i+=((k == -1 || k == realDim[2]-1 ||
                               j==-1 || j==realDim[1]-1) ? 1 : realDim[0]))
                {
                    if (baseIdx[0]+i<0 || baseIdx[0]+i+1 > maxGlobalI) continue; // Skip cells connecting to outside the domain

                    // Indices of the four vertices defining a potential stitch cell
                    int vIdx[8][3] = {
                        { i, j, k },
                        { i + 1, j, k },
                        { i + 1, j + 1, k},
                        { i, j + 1, k },
                        { i, j, k +1 },
                        { i + 1, j, k + 1 },
                        { i + 1, j + 1, k + 1 },
                        { i, j + 1, k + 1 }
                    };

                    // First, compute case number determine for each vertex whether the sub-
                    // cell that connects to it can be added. and check whether cell can be added.
                    // A cell cannot be added if (i) it would extent outside the domain,
                    // (ii) coincides partially or entirely with a refining grid, or (iii)
                    // connects to an adjacent grid in the same level with a smaller domain
                    // number (to avoid duplicate cells, each stitch cell belongs to the 
                    // patch in a level with the smallest number).
                    int caseNo = 0;

                    bool skipCellVtx[8];
                    for (int vtxNo=0; vtxNo < 8; ++vtxNo)
                        skipCellVtx[vtxNo] = false;

                    for (int vtxNo=0; vtxNo < 8; ++vtxNo)
                    {
                        if (// Vertex is in interior, completely refined region:
                                ((vIdx[vtxNo][0] >= 0) && (vIdx[vtxNo][0] < realDim[0]) &&
                                 (vIdx[vtxNo][1] >= 0) && (vIdx[vtxNo][1] < realDim[1]) && 
                                 (vIdx[vtxNo][2] >= 0) && (vIdx[vtxNo][2] < realDim[2])) || 
                                // Vertex belongs to a neighboring refining grid:
                                refinedInSameLevelDomain[LOC(vIdx[vtxNo][0], vIdx[vtxNo][1], vIdx[vtxNo][2])] != -1)
                        {
                            // Vertex is in same level as the patch for which we are
                            // creating stitch cells.

                            // Set bit in case number indicating that this point is refined
                            caseNo |= 1 << vtxNo;

                            // If cell is in same level, it cannot be outside the domain
                            // since it belongs to a fine patch

                            // Check whether it coincides with a fine grid
                            if (avtGhostData::IsGhostZoneType(
                                        inputGhostZones->GetValue(
                                            LOC(vIdx[vtxNo][0], vIdx[vtxNo][1], vIdx[vtxNo][2])),
                                        REFINED_ZONE_IN_AMR_GRID))
                            {
                                skipCellVtx[vtxNo] = true;
                            }

                            // Check whether it connects to a patch in the same level
                            // with a smaller domain number
                            if (// Has to lie on boundary to connect to another patch
                                    ((vIdx[vtxNo][0] == -1) || (vIdx[vtxNo][0] == realDim[0]) ||
                                     (vIdx[vtxNo][1] == -1) || (vIdx[vtxNo][1] == realDim[1]) ||
                                     (vIdx[vtxNo][2] == -1) || (vIdx[vtxNo][2] == realDim[2])) && 
                                    // And have an domain number smaller than us
                                    refinedInSameLevelDomain[LOC(vIdx[vtxNo][0], vIdx[vtxNo][1], vIdx[vtxNo][2])] < domain)
                            {
                                skipCellVtx[vtxNo] = true;
                            }
                        }
                        else
                        {
                            // Vertex is in coarse level
                            int globalI = baseIdx[0] + vIdx[vtxNo][0];
                            int globalJ = baseIdx[1] + vIdx[vtxNo][1];
                            int globalK = baseIdx[2] + vIdx[vtxNo][2];

                            if (globalI < 0 || globalI > maxGlobalI ||
                                globalJ < 0 || globalJ > maxGlobalJ ||
                                globalK < 0 || globalK > maxGlobalK)
                            {
                                skipCellVtx[vtxNo] = true;
                            }
                        }
                    }

#ifdef ENABLE_UNGHOST_CELL_OPTIMIZATION
                    if (caseNo != 255)
#endif
                    {
                        // Compute point locations and store them
                        int subCaseNo = 0;

                        int currIdx[3] = { i, j, k };
                        for (int subDir = 2; subDir >= 0; --subDir)
                            if (tesselationSubCaseDir3D[caseNo] & ( 1 << subDir))
                            {
                                subCaseNo <<= 1;
                                if ((currIdx[subDir] + 1) % refinementRatio[subDir])
                                {
                                    subCaseNo |=  1;
                                }
                            }

                        if (subCaseNo < 0 || subCaseNo > 4)
                        {
                            EXCEPTION1(VisItException,
                                    "Internal error: Invalid subcase.");
                        }

                        int tessPos = tesselationCaseStart3D[caseNo][subCaseNo];
                        if (tessPos == -1)
                        {
                            EXCEPTION1(VisItException,
                                    "Internal error: Invlid case/subcase combintation");
                        }

                        while (tesselationArray3D[tessPos] != -1)
                        {
                            int numVtcs = tesselationArray3D[tessPos++];
                            std::vector<vtkIdType> vPtId(numVtcs);
                            std::vector<vtkIdType> vSrcPtId(numVtcs);

                            bool skipCell = false;
                            for (int vtxNo = 0; vtxNo < numVtcs; ++ vtxNo)
                                if (skipCellVtx[tesselationArray3D[tessPos+vtxNo]])
                                {
                                    skipCell = true;
                                    break;
                                }

                            if (skipCell)
                            {
                                tessPos += numVtcs;
                                continue; // Skip on to next cell
                            }

                            for (int vtxNo =0; vtxNo < numVtcs; ++vtxNo)
                            {
                                int cubeVtxNo = tesselationArray3D[tessPos++];
                                int vI = vIdx[cubeVtxNo][0];
                                int vJ = vIdx[cubeVtxNo][1];
                                int vK = vIdx[cubeVtxNo][2];
                                vSrcPtId[vtxNo] = LOC(vI, vJ, vK);

                                if (pointIds[LOC(vI, vJ, vK)] < 0)
                                {
                                    // Vertex was not computed for a previous stitch cell
                                    if (caseNo & (1 << cubeVtxNo))
                                    {
                                        // Fine cell
                                        pointIds[LOC(vI, vJ, vK)] = stitchCellPts->InsertNextPoint(
                                                domainOrigin[0] + (baseIdx[0] + vI  + 0.5) * cellSize[level][0],
                                                domainOrigin[1] + (baseIdx[1] + vJ  + 0.5) * cellSize[level][1],
                                                domainOrigin[2] + (baseIdx[2] + vK  + 0.5) * cellSize[level][2]
                                                );
                                        originalDataIds.push_back(vtkIdType(LOC(vI, vJ, vK)));
                                    }
                                    else
                                    {
                                        // Vertex i and j in parent level indices
                                        int pI = (baseIdx[0] + vI) / refinementRatio[0];
                                        int pJ = (baseIdx[1] + vJ) / refinementRatio[1];
                                        int pK = (baseIdx[2] + vK) / refinementRatio[2];

                                        pointIds[LOC(vI,vJ,vK)] = stitchCellPts->InsertNextPoint(
                                                domainOrigin[0] + (pI + 0.5) * cellSize[level-1][0],
                                                domainOrigin[1] + (pJ + 0.5) * cellSize[level-1][1],
                                                domainOrigin[2] + (pK + 0.5) * cellSize[level-1][2]
                                                );
                                        originalDataIds.push_back(vtkIdType(LOC(vI, vJ, vK)));

                                        if (((vI == -1) || (vI == realDim[0])) &&
                                            (vJ != -1) && (vJ != realDim[1]) &&
                                            (vK != -1) && (vK != realDim[2]))
                                        {
                                            int baseJ = vJ - (vJ % refinementRatio[1]);
                                            int baseK = vK - (vK % refinementRatio[2]);

                                            for (int refJVtxNo = 0; refJVtxNo < refinementRatio[1];
                                                    ++refJVtxNo)
                                                for (int refKVtxNo = 0; refKVtxNo < refinementRatio[2];
                                                        ++refKVtxNo)
                                                {
                                                    pointIds[LOC(vI, baseJ+refJVtxNo, baseK+refKVtxNo)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                                }
                                        }
                                        else if (((vJ == -1) || (vJ == realDim[1])) &&
                                                 (vI != -1) && (vI != realDim[0]) &&
                                                 (vK != -1) && (vK != realDim[2]))
                                        {
                                            int baseI = vI - (vI % refinementRatio[0]);
                                            int baseK = vK - (vK % refinementRatio[2]);

                                            for (int refIVtxNo = 0; refIVtxNo < refinementRatio[0];
                                                    ++refIVtxNo)
                                                for (int refKVtxNo = 0; refKVtxNo < refinementRatio[2];
                                                        ++refKVtxNo)
                                                {
                                                    pointIds[LOC(baseI+refIVtxNo, vJ, baseK+refKVtxNo)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                                }
                                        }
                                        else if (((vK == -1) || (vK == realDim[2])) &&
                                                 (vI != -1) && (vI != realDim[0]) &&
                                                 (vJ != -1) && (vJ != realDim[1]))
                                        {
                                            int baseI = vI - (vI % refinementRatio[0]);
                                            int baseJ = vJ - (vJ % refinementRatio[1]);

                                            for (int refIVtxNo = 0; refIVtxNo < refinementRatio[0];
                                                    ++refIVtxNo)
                                                for (int refJVtxNo = 0; refJVtxNo < refinementRatio[1];
                                                        ++refJVtxNo)
                                                {
                                                    pointIds[LOC(baseI+refIVtxNo, baseJ+refJVtxNo, vK)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                                }
                                        }

                                        else if (((vI == -1) || (vI == realDim[0])) &&
                                                 ((vJ == -1) || (vJ == realDim[1])) &&
                                                 (vK != -1) && (vK != realDim[2]))
                                        {
                                            int baseK = vK - (vK % refinementRatio[2]);
                                            for (int refKVtxNo = 0; refKVtxNo < refinementRatio[2];
                                                    ++refKVtxNo)
                                                    pointIds[LOC(vI, vJ, baseK+refKVtxNo)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                        } 
                                        else if (((vI == -1) || (vI == realDim[0]) ||
                                                  (vK == -1) || (vK == realDim[2])) &&
                                                 (vJ != -1) && (vJ != realDim[1]))
                                        {
                                            int baseJ = vJ - (vJ % refinementRatio[1]);
                                            for (int refJVtxNo = 0; refJVtxNo < refinementRatio[1];
                                                    ++refJVtxNo)
                                                    pointIds[LOC(vI, baseJ+refJVtxNo, vK)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                        }
                                        else if (((vJ == -1) || (vJ == realDim[1]) ||
                                                  (vK == -1) || (vK == realDim[2])) &&
                                                 (vI != -1) && (vI != realDim[0]))
                                        {
                                            int baseI = vI - (vI % refinementRatio[0]);
                                            for (int refIVtxNo = 0; refIVtxNo < refinementRatio[0];
                                                    ++refIVtxNo)
                                                    pointIds[LOC(baseI+refIVtxNo, vJ, vK)] = 
                                                        pointIds[LOC(vI,vJ,vK)];
                                        }
                                    }
                                }

                                // Copy previously or newly computed point
                                vPtId[vtxNo] = pointIds[LOC(vI, vJ, vK)];
                            }

                            if (numVtcs == 4)
                            {
                                ugrid->InsertNextCell(VTK_TETRA, 4, &vPtId[0]);
                            }
                            else if (numVtcs == 5)
                            {
                                ugrid->InsertNextCell(VTK_PYRAMID, 5, &vPtId[0]);
                            }
                            else if (numVtcs == 6)
                            {
                                ugrid->InsertNextCell(VTK_WEDGE, 6, &vPtId[0]);
                            }
                            else if (numVtcs == 7)
                            {
                                double centroid[3] = { 0, 0, 0 };
                                for (int cPtNo=0; cPtNo<7; ++cPtNo)
                                {
                                    double pt[3];
                                    stitchCellPts->GetPoint(vPtId[cPtNo], pt);
                                    for (int dim=0; dim<3; ++dim)
                                        centroid[dim]+=pt[dim];
                                }
                                for (int dim=0; dim<3; ++dim)
                                    centroid[dim]/=7.;
                                vtkIdType centroidId = stitchCellPts->InsertNextPoint(centroid);
                                //std::cout << "Centroid ID: " << centroidId << std::endl;
                                originalDataIds.push_back(-1);
                                interpolatedDataIds.push_back(vtkIdList::New());
                                interpolatedDataIds.back()->SetNumberOfIds(7);
                                //std::cout << "Interpolated points: ";
                                for (vtkIdType cPtNo=0; cPtNo<7; ++cPtNo)
                                {
                                    //std::cout << vSrcPtId[cPtNo] << " (ID="<<vPtId[cPtNo]<<") ";
                                    interpolatedDataIds.back()->SetId(cPtNo, vSrcPtId[cPtNo]);
                                }
                                //std::cout << std::endl;

                                vtkIdType c0[5] = { vPtId[0], vPtId[1], vPtId[2], vPtId[3], centroidId };
                                vtkIdType c1[5] = { vPtId[0], vPtId[4], vPtId[5], vPtId[1], centroidId };
                                vtkIdType c2[5] = { vPtId[1], vPtId[5], vPtId[6], vPtId[2], centroidId };
                                vtkIdType c3[4] = { vPtId[2], vPtId[6], vPtId[3], centroidId };
                                vtkIdType c4[4] = { vPtId[3], vPtId[6], vPtId[4], centroidId };
                                vtkIdType c5[4] = { vPtId[0], vPtId[3], vPtId[4], centroidId };
                                vtkIdType c6[4] = { vPtId[4], vPtId[6], vPtId[5], centroidId };
                                ugrid->InsertNextCell(VTK_PYRAMID, 5, c0);
                                ugrid->InsertNextCell(VTK_PYRAMID, 5, c1);
                                ugrid->InsertNextCell(VTK_PYRAMID, 5, c2);
                                ugrid->InsertNextCell(VTK_TETRA, 4, c3);
                                ugrid->InsertNextCell(VTK_TETRA, 4, c4);
                                ugrid->InsertNextCell(VTK_TETRA, 4, c5);
                                ugrid->InsertNextCell(VTK_TETRA, 4, c6);
                            }
                            else if (numVtcs == 8)
                            {
                                ugrid->InsertNextCell(VTK_HEXAHEDRON, 8, &vPtId[0]);
                            }
                            else
                                EXCEPTION1(VisItException, "Invalid cell type (internal error).");
                        }
                    }
                }
            }
        }
    }


    // Copy attributes
    ugrid->GetPointData()->InterpolateAllocate(rgrid->GetCellData(), originalDataIds.size());
    int idx = 0;
    std::list<vtkIdList*>::const_iterator iVLIt = interpolatedDataIds.begin();
    vtkIdList *oneIdList = vtkIdList::New();
    oneIdList->SetNumberOfIds(1);
    double oneWeight[1] = { 1.0 };
    double sevenWeights[7] = { 1./7., 1./7., 1./7., 1./7., 1./7., 1./7., 1./7. };

    // FIXME: Possibly use combination of Copy and Interpolate point to improve performance.
    // However, VTK does not seem to support mixing these calls.
    for (std::list<vtkIdType>::const_iterator it = originalDataIds.begin();
            it != originalDataIds.end(); ++it)
    {
        if (*it != -1)
        {
            //std::cout << idx << ": Copying point with ID: " << *it << std::endl;
            oneIdList->SetId(0, *it);
            ugrid->GetPointData()->InterpolatePoint(rgrid->GetCellData(), idx, oneIdList, oneWeight);
        }
        else
        {
            ugrid->GetPointData()->InterpolatePoint(rgrid->GetCellData(), idx, *iVLIt, sevenWeights);
            (*iVLIt)->Delete();
            ++iVLIt;
        }
        ++idx;
    }

    // Clean-up
    oneIdList->Delete(); // FIXME: Move to regular clean-up after debugging is done

    // Clean-up
    delete[] pointIds;
    return ugrid;
}

// ****************************************************************************
//  Method: avtAMRStitchCellFilter::CreateDualGrid
//
//  Purpose:
//      Create dual grid from input data set
//
//  Arguments:
//      in_ds      The input dataset.
//      domain     The domain number.
//
//  Returns:       The output dataset.
//
//  Programmer: Gunther H. Weber (based on Cyrus Harrison's DualMesh filter)
//  Creation:   Thu Jul 8 15:14:01 PST 2010
//
//  Modifications:
//    Kathleen Biagas, Tue Oct 16 13:27:01 MST 2012
//    Preserve coordinate type.
//
// ****************************************************************************

vtkDataSet *
avtAMRStitchCellFilter::CreateDualGrid(vtkRectilinearGrid  *rgrid, int domain,
        int level, int dims[3], int baseIdx[3], int ghostOffset[3],
        const vtkIdType* refinedInSameLevelDomain)
{

    // Create reuslt mesh
    vtkRectilinearGrid *result = vtkRectilinearGrid::New();
    // VTK expects number of samples not cells. The number of cells stored in our
    // dims[] corresponds to the number of samples of the dual grid.
    result->SetDimensions(dims);

    vtkDataArray *xCoords = rgrid->GetXCoordinates()->NewInstance();
    for (int i = 0; i<dims[0]; ++i)
        xCoords->InsertNextTuple1(
                domainOrigin[0] + (i + baseIdx[0] - ghostOffset[0] +  0.5) * cellSize[level][0]
                );
    result->SetXCoordinates(xCoords);
    xCoords->Delete();

    vtkDataArray *yCoords = rgrid->GetYCoordinates()->NewInstance();
    for (int j = 0; j<dims[1]; ++j)
        yCoords->InsertNextTuple1(
                domainOrigin[1] + (j + baseIdx[1] - ghostOffset[1] +  0.5) * cellSize[level][1]
                );
    result->SetYCoordinates(yCoords);
    yCoords->Delete();

    if (topologicalDimension == 3)
    {
        vtkDataArray *zCoords = rgrid->GetZCoordinates()->NewInstance();
        for (int k = 0; k<dims[2]; ++k)
            zCoords->InsertNextTuple1(
                    domainOrigin[2] + (k + baseIdx[2] - ghostOffset[2] +  0.5) * cellSize[level][2]
                    );
        result->SetZCoordinates(zCoords);
        zCoords->Delete();
    }

    // Shallow copy the field data
    result->GetFieldData()->ShallowCopy(rgrid->GetFieldData());

    // Computing the dual grid invalidates field data
    // -> Remove that information
    result->GetFieldData()->RemoveArray("avtRealDims");

    // Cell data of the original grid becomse point data of the 
    // dual grid
    result->GetPointData()->ShallowCopy(rgrid->GetCellData());

    // Update ghost zones
    vtkUnsignedCharArray *inputGhostZones =
        dynamic_cast<vtkUnsignedCharArray *>(rgrid->GetCellData()->GetArray("avtGhostZones"));

    if (!inputGhostZones)
        EXCEPTION1(ImproperUseException,
                "Need ghost zone information to compute dual mesh.");

    vtkUnsignedCharArray *outputGhostZones = vtkUnsignedCharArray::New();
    outputGhostZones->SetName("avtGhostZones");

    if (topologicalDimension == 2)
    {
        outputGhostZones->SetNumberOfTuples((dims[0]-1)*(dims[1]-1));

        // 2D
        for (int i=0; i < dims[0]-1; ++i)
            for (int j=0; j < dims[1]-1; ++j)
            {
#ifdef IGV
#undef IGV
#endif
#define IGV(i, j) inputGhostZones->GetValue((j)*dims[0]+i)
                unsigned char ghostInfo = IGV(i, j);
                ghostInfo |= IGV(i+1, j);
                ghostInfo |= IGV(i,   j+1);
                ghostInfo |= IGV(i+1, j+1);
#undef IGV

#ifdef IS_REMOVE_GHOST_SAMPLE
#undef IS_REMOVE_GHOST_SAMPLE
#endif
#define IS_REMOVE_GHOST_SAMPLE(i,j) \
                (!avtGhostData::IsGhostZone(inputGhostZones->GetValue((j)*dims[0]+(i))) || \
                 (refinedInSameLevelDomain[(j)*dims[0]+(i)] > domain))
                if (
#ifndef ENABLE_UNGHOST_CELL_OPTIMIZATION
                    (level == 0) &&
#endif
                    ((i == 0) || (i == dims[0] - 2) ||
                     (j == 0) || (j == dims[1] - 2)) &&
                    IS_REMOVE_GHOST_SAMPLE(i  , j) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j) &&
                    IS_REMOVE_GHOST_SAMPLE(i,   j+1) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j+1) &&
                    !(avtGhostData::IsGhostZoneType(ghostInfo,
                            REFINED_ZONE_IN_AMR_GRID)))
                {
                    ghostInfo = 0;
                }
#undef IS_REMOVE_GHOST_SAMPLE

                outputGhostZones->SetValue(j*(dims[0]-1)+i, ghostInfo);
            }
    }
    else
    {
        // 3D
       outputGhostZones->SetNumberOfTuples((dims[0]-1)*(dims[1]-1)*(dims[2]-1));

        for (int i=0; i < dims[0]-1; ++i)
            for (int j=0; j < dims[1]-1; ++j)
                for (int k=0; k < dims[2]-1; ++k)
                {
#ifdef IGV
#undef IGV
#endif
#define IGV(i, j, k) inputGhostZones->GetValue(((k)*dims[1]+j)*dims[0]+i)
                    unsigned char ghostInfo = IGV(i, j, k);
                    ghostInfo |= IGV(i+1, j  , k);
                    ghostInfo |= IGV(i,   j+1, k);
                    ghostInfo |= IGV(i+1, j+1, k);
                    ghostInfo |= IGV(i,   j,   k+1);
                    ghostInfo |= IGV(i+1, j,   k+1);
                    ghostInfo |= IGV(i,   j+1, k+1);
                    ghostInfo |= IGV(i+1, j+1, k+1);
#undef IGV

#ifdef IS_REMOVE_GHOST_SAMPLE
#undef IS_REMOVE_GHOST_SAMPLE
#endif
#define IS_REMOVE_GHOST_SAMPLE(i,j,k) \
                (!avtGhostData::IsGhostZone(inputGhostZones->GetValue(((k)*dims[1]+j)*dims[0]+i)) || \
                 (refinedInSameLevelDomain[((k)*dims[1]+j)*dims[0]+i] > domain))
                if (
#ifndef ENABLE_UNGHOST_CELL_OPTIMIZATION
                    (level == 0) &&
#endif
                    ((i == 0) || (i == dims[0] - 2) ||
                     (j == 0) || (j == dims[1] - 2) ||
                     (k == 0) || (k == dims[2] - 2)) &&
                    IS_REMOVE_GHOST_SAMPLE(i,   j,   k) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j,   k) &&
                    IS_REMOVE_GHOST_SAMPLE(i,   j+1, k) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j+1, k) &&
                    IS_REMOVE_GHOST_SAMPLE(i,   j,   k+1) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j,   k+1) &&
                    IS_REMOVE_GHOST_SAMPLE(i,   j+1, k+1) &&
                    IS_REMOVE_GHOST_SAMPLE(i+1, j+1, k+1) &&
                    !(avtGhostData::IsGhostZoneType(ghostInfo,
                            REFINED_ZONE_IN_AMR_GRID)))
                {
                    ghostInfo = 0;
                }
#undef IS_REMOVE_GHOST_SAMPLE

                outputGhostZones->SetValue((k*(dims[1]-1)+j)*(dims[0]-1)+i, ghostInfo);
            }
    }
 
    result->GetCellData()->AddArray(outputGhostZones);

    return result;
}
