/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkSubdivideTetra.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkSubdivideTetra.h"

#include "vtkCellType.h"
#include "vtkGenericCell.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMergePoints.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkUnstructuredGrid.h"

vtkStandardNewMacro(vtkSubdivideTetra);

//----------------------------------------------------------------------------
// Description:
// Construct with all types of clipping turned off.
vtkSubdivideTetra::vtkSubdivideTetra()
{
}

//----------------------------------------------------------------------------
int vtkSubdivideTetra::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // get the input and output
  vtkUnstructuredGrid *input = vtkUnstructuredGrid::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkIdType numPts = input->GetNumberOfPoints();
  vtkIdType numCells = input->GetNumberOfCells();
  vtkPoints *inPts = input->GetPoints();
  vtkIdType cellId, i;
  vtkIdType pts[4];
  vtkGenericCell *cell;
  vtkPointData *pd = input->GetPointData();
  vtkPointData *outputPD = output->GetPointData();
  vtkPoints *newPts;
  vtkIdType ptId;

  double weights[4], x0[3], x1[3], x2[3], x3[3], x[3];
  int p0, p1, p2, p3;
  vtkIdType center, e01, e02, e03, e12, e13, e23;
  vtkMergePoints *locator;

  vtkDebugMacro(<<"Executing mesh subdivide");

  if (input->IsHomogeneous() == 0 ||
      input->GetCellType(0) != VTK_TETRA)
  {
    vtkErrorMacro(<<"all cells must be tetrahedra.");
    return 1;
  }

  // Copy original points and point data
  newPts = vtkPoints::New();
  newPts->Allocate(5*numPts,numPts);
  outputPD->InterpolateAllocate(pd,5*numPts,numPts);

  output->Allocate(numCells);
  output->SetPoints(newPts);

  locator = vtkMergePoints::New();
  locator->InitPointInsertion (newPts, input->GetBounds());

  for (ptId=0; ptId < numPts; ptId++)
  {
    locator->InsertNextPoint(inPts->GetPoint(ptId));
    outputPD->CopyData(pd,ptId,ptId);
  }

  cell = vtkGenericCell::New();

  // loop over tetrahedra, generating sixteen new ones for each. This is
  // done by introducing mid-edge nodes and a single mid-tetra node.
  for(cellId=0; cellId < numCells; cellId++)
  {
    input->GetCell(cellId, cell);

    // get tetra points
    cell->Points->GetPoint(0,x0);
    cell->Points->GetPoint(1,x1);
    cell->Points->GetPoint(2,x2);
    cell->Points->GetPoint(3,x3);

    p0 = cell->PointIds->GetId(0);
    p1 = cell->PointIds->GetId(1);
    p2 = cell->PointIds->GetId(2);
    p3 = cell->PointIds->GetId(3);

    // compute center point
    weights[0] = weights[1] = weights[2] = weights[3] = 0.25;
    for (i=0; i<3; i++)
    {
      x[i] = 0.25*(x0[i] + x1[i] + x2[i] + x3[i]);
    }
    center = locator->InsertNextPoint(x);
    outputPD->InterpolatePoint(pd, center, cell->PointIds, weights);

    // compute edge points
    // edge 0-1
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x1[i] + x0[i]);
    }
    e01 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e01, p0, p1, 0.5);

    // edge 1-2
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x2[i] + x1[i]);
    }
    e12 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e12, p1, p2, 0.5);

    // edge 2-0
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x2[i] + x0[i]);
    }
    e02 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e02, p2, p0, 0.5);

    // edge 0-3
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x3[i] + x0[i]);
    }
    e03 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e03, p0, p3, 0.5);

    // edge 1-3
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x3[i] + x1[i]);
    }
    e13 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e13, p1, p3, 0.5);

    // edge 2-3
    for (i=0; i<3; i++)
    {
      x[i] = 0.5 * (x3[i] + x2[i]);
    }
    e23 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e23, p2, p3, 0.5);

    // Now create tetrahedra
    // First, four tetra from each vertex
    pts[0] = p0;
    pts[1] = e01;
    pts[2] = e02;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p1;
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p2;
    pts[1] = e02;
    pts[2] = e12;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p3;
    pts[1] = e03;
    pts[2] = e13;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);

    // Now four tetra from cut-off tetra corners
    pts[0] = center;
    pts[1] = e01;
    pts[2] = e02;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e02;
    pts[2] = e12;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e03;
    pts[2] = e13;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);

    // Now four tetra from triangles on tetra faces
    pts[0] = center;
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e02;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e01;
    pts[2] = e13;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e12;
    pts[2] = e23;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e02;
    pts[2] = e23;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);

  } //for all cells
  cell->Delete();

  vtkDebugMacro(<<"Subdivided " << numCells << " cells");

  locator->Delete();
  newPts->Delete();
  output->Squeeze();

  return 1;
}

//----------------------------------------------------------------------------
void vtkSubdivideTetra::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
}