//------------------------------------------------------------------------------
//
// Project:	Crack-free Isosurfaces for AMR Data
// Module:	$RCSfile$
// Language:	C++
// Date:	$Date$
// Author:	$Author$
// Version:	$Revision$
//
//------------------------------------------------------------------------------

#include "AMRGridStitcher.hh"

// Own Helper Includes
#include <assume.hh>
#include <Vertex.hh>

// AMR Specific Includes
#include <AMRGrid.hh>
#include <AMRLevel.hh>

namespace AMRGridStitcherAuxFct
{
  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism045(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      stitchCellHandler.handleTetrahedronCell(p3, p2, p5, p0);
      stitchCellHandler.handlePyramidCell(p1, p4, p5, p2, p0);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism145(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      stitchCellHandler.handleTetrahedronCell(p3, p2, p5, p1);
      stitchCellHandler.handlePyramidCell(p0, p3, p5, p4, p1);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism245(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      generateCellsForDeformedTrianglePrism045(stitchCellHandler, p2, p3, p0, p1, p5, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism345(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      generateCellsForDeformedTrianglePrism145(stitchCellHandler, p2, p3, p0, p1, p5, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism0145(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      stitchCellHandler.handleTetrahedronCell(p0, p4, p1, p5);
      stitchCellHandler.handlePyramidCell(p0, p1, p2, p3, p5);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism2345(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5) 
    {
      generateCellsForDeformedTrianglePrism0145(stitchCellHandler, p2, p3, p0, p1, p5, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism12345(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p4, p1, p3);
      stitchCellHandler.handlePyramidCell(p1, p4, p5, p2, p3);
    }
  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism02345(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p4, p1, p2);
      stitchCellHandler.handlePyramidCell(p0, p3, p5, p4, p2);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism01345(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5)
    {
      generateCellsForDeformedTrianglePrism12345(stitchCellHandler, p2, p3, p0, p1, p5, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedTrianglePrism01245(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5)
    {
      generateCellsForDeformedTrianglePrism02345(stitchCellHandler, p2, p3, p0, p1, p5, p4);
    }
  template <class StitchCellHandler>
    inline void generateCellsForDeformedPyramid04(StitchCellHandler& stitchCellHandler,  const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p1, p2, p4);
      stitchCellHandler.handleTetrahedronCell(p0, p2, p3, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedPyramid023(StitchCellHandler& stitchCellHandler,  const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p1, p2, p4);
      stitchCellHandler.handleTetrahedronCell(p0, p2, p3, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedPyramid024(StitchCellHandler& stitchCellHandler,  const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p1, p2, p4);
      stitchCellHandler.handleTetrahedronCell(p0, p2, p3, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedPyramid1234(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4)
    {
      stitchCellHandler.handleTetrahedronCell(p0, p1, p3, p4);
      stitchCellHandler.handleTetrahedronCell(p1, p2, p3, p4);
    }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedCube024567(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5, const AMRGridCell& p6, const AMRGridCell& p7)
    {
      generateCellsForDeformedPyramid1234(stitchCellHandler, p3, p7, p4, p0, p2);
      generateCellsForDeformedPyramid1234(stitchCellHandler, p1, p0, p4, p5, p2);
      stitchCellHandler.handlePyramidCell(p4, p7, p6, p5, p2);
    }

  Vertex centroid(const Vertex& v0, const Vertex& v1, const Vertex& v2, const Vertex& v3, const Vertex& v4, const Vertex& v5, const Vertex& v6)
  {
    return affineCombination(7, &v0, 1.0/7.0, &v1, 1.0/7.0, &v2, 1.0/7.0, &v3, 1.0/7.0, &v4, 1.0/7.0, &v5, 1.0/7.0, &v6, 1.0/7.0);
  }

  template <class StitchCellHandler>
    inline void generateCellsForDeformedCube1234567(StitchCellHandler& stitchCellHandler, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3, const AMRGridCell& p4, const AMRGridCell& p5, const AMRGridCell& p6, const AMRGridCell& p7)
    {
      stitchCellHandler.handleTetrahedronCell(p3, p1, p4, p0);
      std::cout << "Inserting vertex to avoid illegal tesselation." << std::endl;
      Vertex newPos = AMRGridStitcherAuxFct::centroid(p1.pos(), p2.pos(), p3.pos(), p4.pos(), p5.pos(), p6.pos(), p7.pos());
      double newVal = p1.val() + p2.val() + p3.val() + p4.val() + p5.val() + p6.val() + p7.val();
      newVal /= 7;
      stitchCellHandler.handleTetrahedronCell(p3, p4, p1, newPos, newVal);
      stitchCellHandler.handleTetrahedronCell(p4, p3, p7, newPos, newVal);
      stitchCellHandler.handlePyramidCell(p4, p7, p6, p5, newPos, newVal);
      stitchCellHandler.handlePyramidCell(p3, p2, p6, p7, newPos, newVal);
      stitchCellHandler.handleTetrahedronCell(p1, p2, p3, newPos, newVal);
      stitchCellHandler.handlePyramidCell(p1, p5, p6, p2, newPos, newVal);
      stitchCellHandler.handleTetrahedronCell(p1, p4, p5, newPos, newVal);
    }
}

template <class StitchCellHandler>
AMRGridStitcher<StitchCellHandler>::AMRGridStitcher(StitchCellHandler& stitchCellHandler, AMRHierarchy& amrHier, AMRGrid *stitchGrid, bool considerNextLevel)
: mStitchCellHandler(stitchCellHandler), mAMRDataSet(amrHier), mStitchGrid(stitchGrid), mConsiderNextLevel(considerNextLevel), mParentLevel(mAMRDataSet.getLevel(stitchGrid->level-1))
{
  assert(mParentLevel != 0);
  // <ATTENTION> Assume spatialrefinement == gridplacementrefinement </ATTENTION>
  // <ATTENTION> Assume spatialrefinement % parentLevel->getSpatialRefinemet == 0 </ATTENTION>
  int refinementWRTParentLevel[3];
  refinementWRTParentLevel[xAxis] = mStitchGrid->spatialrefinement[xAxis] / mParentLevel->getSpatialRefinement(xAxis);
  refinementWRTParentLevel[yAxis] = mStitchGrid->spatialrefinement[yAxis] / mParentLevel->getSpatialRefinement(yAxis);
  refinementWRTParentLevel[zAxis] = mStitchGrid->spatialrefinement[zAxis] / mParentLevel->getSpatialRefinement(zAxis);

  mOriginInParentLevel[xAxis] = mStitchGrid->iorigin[xAxis] / refinementWRTParentLevel[xAxis];
  mOriginInParentLevel[yAxis] = mStitchGrid->iorigin[yAxis] / refinementWRTParentLevel[yAxis];
  mOriginInParentLevel[zAxis] = mStitchGrid->iorigin[zAxis] / refinementWRTParentLevel[zAxis];
  // <ATTENTION> Assume refined grid starts at cell boundaries of parent level </ATTENTION>
  assume(mStitchGrid->iorigin[xAxis] %  refinementWRTParentLevel[xAxis] == 0);
  assume(mStitchGrid->iorigin[yAxis] %  refinementWRTParentLevel[yAxis] == 0);
  assume(mStitchGrid->iorigin[zAxis] %  refinementWRTParentLevel[zAxis] == 0);

  mMaxExtInParentLevel[xAxis] = mOriginInParentLevel[xAxis] + (mStitchGrid->dims[xAxis] / refinementWRTParentLevel[xAxis]) -1;
  mMaxExtInParentLevel[yAxis] = mOriginInParentLevel[yAxis] + (mStitchGrid->dims[yAxis] / refinementWRTParentLevel[yAxis]) -1;
  mMaxExtInParentLevel[zAxis] = mOriginInParentLevel[zAxis] + (mStitchGrid->dims[zAxis] / refinementWRTParentLevel[zAxis]) -1;
  // <ATTENTION> Assume refined grid ends at cell boundaries of parent level </ATTENTION>
  assume(mStitchGrid->dims[xAxis] %  refinementWRTParentLevel[xAxis] == 0);
  assume(mStitchGrid->dims[yAxis] %  refinementWRTParentLevel[yAxis] == 0);
  assume(mStitchGrid->dims[zAxis] %  refinementWRTParentLevel[zAxis] == 0);

  // Load all grids necessary
  AMRLevel *thisLevel = mAMRDataSet.getLevel(mStitchGrid->level);
  for (AMRGrid *currGridPtr = thisLevel->getFirstGrid(); currGridPtr != 0; currGridPtr = currGridPtr->getNextGridInLevel())
    currGridPtr->loadCells();

  for (AMRGrid *currGridPtr = mParentLevel->getFirstGrid(); currGridPtr != 0; currGridPtr = currGridPtr->getNextGridInLevel())
    currGridPtr->loadCells();
}

template <class StitchCellHandler>
AMRGridStitcher<StitchCellHandler>::~AMRGridStitcher()
{
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::generateStitchCells(bool xyMin, bool xyMax, bool yzMin, bool yzMax, bool xzMin, bool xzMax)
{
  if (xyMin) generateCellsForFace(xAxis, yAxis, zAxis, min);
  if (xzMin) generateCellsForFace(zAxis, xAxis, yAxis, min);
  if (yzMin) generateCellsForFace(yAxis, zAxis, xAxis, min);
  if (xyMax) generateCellsForFace(yAxis, xAxis, zAxis, max);
  if (xzMax) generateCellsForFace(xAxis, zAxis, yAxis, max);
  if (yzMax) generateCellsForFace(zAxis, yAxis, xAxis, max);
  if (xyMin || xzMin) generateCellsForEdge(xAxis, yAxis, min, zAxis, min);
  if (xyMax || xzMin) generateCellsForEdge(xAxis, yAxis, min, zAxis, max);
  if (xyMin || xzMax) generateCellsForEdge(xAxis, yAxis, max, zAxis, min);
  if (xyMax || xzMax) generateCellsForEdge(xAxis, yAxis, max, zAxis, max);
  if (xyMin || yzMin) generateCellsForEdge(yAxis, zAxis, min, xAxis, min);
  if (xyMin || yzMax) generateCellsForEdge(yAxis, zAxis, min, xAxis, max);
  if (xyMax || yzMin) generateCellsForEdge(yAxis, zAxis, max, xAxis, min);
  if (xyMax || yzMax) generateCellsForEdge(yAxis, zAxis, max, xAxis, max);
  if (xzMin || yzMin) generateCellsForEdge(zAxis, xAxis, min, yAxis, min);
  if (xzMin || yzMax) generateCellsForEdge(zAxis, xAxis, min, yAxis, max);
  if (xzMax || yzMin) generateCellsForEdge(zAxis, xAxis, max, yAxis, min);
  if (xzMax || yzMax) generateCellsForEdge(zAxis, xAxis, max, yAxis, max);
  generateCellsForVertex(xAxis, yAxis, zAxis, min, min, min);
  generateCellsForVertex(zAxis, yAxis, xAxis, min, min, max);
  generateCellsForVertex(xAxis, yAxis, zAxis, max, min, max);
  generateCellsForVertex(zAxis, yAxis, xAxis, max, min, min);
  generateCellsForVertex(zAxis, yAxis, xAxis, min, max, min);
  generateCellsForVertex(xAxis, yAxis, zAxis, max, max, min);
  generateCellsForVertex(zAxis, yAxis, xAxis, max, max, max);
  generateCellsForVertex(xAxis, yAxis, zAxis, min, max, max);
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectQuadToVertex(const AMRGridCell& qp0, const AMRGridCell& qp1, const AMRGridCell& qp2, const AMRGridCell& qp3, const AMRGridCell& vp, AxisType axis0, AxisType axis1, AxisType axis2, MinMaxType faceType)
{
  // If the coarse grid points are from different grids, i.e. the quad results from several 
  // coarse grids rather than a single one, only generate stitch cells if we are stitching
  // the grid with the smallest index.
  if (qp0.containingId() < mStitchGrid->index() || qp1.containingId() < mStitchGrid->index() ||
      qp2.containingId() < mStitchGrid->index() || qp3.containingId() < mStitchGrid->index()) return;

  if (vp.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= vp.refiningId()) {
      // No cell bounadaries are crossed, check only one cell for refinement
      AMRGrid *otherGrid = mAMRDataSet.getGridByIndex(vp.refiningId());
      AxisAlignedBox<int> refiningCells = otherGrid->getRefiningCells(vp.grid(), vp.idx(xAxis), vp.idx(yAxis), vp.idx(zAxis));
      int oIdx00[3];
      oIdx00[axis0] = refiningCells.min(axis0);
      oIdx00[axis1] = refiningCells.min(axis1);
      oIdx00[axis2] = faceType == min ? refiningCells.max(axis2) : refiningCells.min(axis2);

      int oIdx01[3];
      oIdx01[axis0] = oIdx00[axis0] + 1;
      oIdx01[axis1] = oIdx00[axis1];
      oIdx01[axis2] = oIdx00[axis2];

      int oIdx10[3];
      oIdx10[axis0] = oIdx00[axis0];
      oIdx10[axis1] = oIdx00[axis1] + 1;
      oIdx10[axis2] = oIdx00[axis2];

      int oIdx11[3];
      oIdx11[axis0] = oIdx00[axis0] + 1;
      oIdx11[axis1] = oIdx00[axis1] + 1;
      oIdx11[axis2] = oIdx00[axis2];

      AMRGridCell vpRef0(otherGrid, oIdx00);
      AMRGridCell vpRef1(otherGrid, oIdx10);
      AMRGridCell vpRef2(otherGrid, oIdx11);
      AMRGridCell vpRef3(otherGrid, oIdx01);

      if (!mConsiderNextLevel ||
	  vpRef0.refiningId() == AMRGrid::NoGrid &&
	  vpRef1.refiningId() == AMRGrid::NoGrid &&
	  vpRef2.refiningId() == AMRGrid::NoGrid &&
	  vpRef3.refiningId() == AMRGrid::NoGrid)
	// Only add stitch cell when no stitch cells for next level are generated (i.e.
	// one of the refining cells is refined itself by an even finer grid)
	mStitchCellHandler.handleHexahedronCell(qp0, qp1, qp2, qp3, vpRef0, vpRef1, vpRef2, vpRef3);
    }
    // else Do nothing!
  }
  else {
    mStitchCellHandler.handlePyramidCell(qp0, qp1, qp2, qp3, vp);
  }
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectQuadToEdge(const AMRGridCell& qp0, const AMRGridCell& qp1, const AMRGridCell& qp2, const AMRGridCell& qp3, const AMRGridCell& ep0, const AMRGridCell& ep1, AxisType axis01, AxisType axis03, AxisType axisPerp0123, MinMaxType qp0Axis01Type, MinMaxType qp0Axis03Type, MinMaxType quadFaceType)
{
  // If the coarse grid points are from different grids, i.e. the quad results from several 
  // coarse grids rather than a single one, only generate stitch cells if we are stitching
  // the grid with the smallest index.
  if (qp0.containingId() < mStitchGrid->index() || qp1.containingId() < mStitchGrid->index() ||
      qp2.containingId() < mStitchGrid->index() || qp3.containingId() < mStitchGrid->index()) return;

  int connectCase = 0;

  AMRGridCell *epRef0 = 0;
  AMRGridCell *epRef1 = 0;
  if (ep0.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 1;

    AMRGrid *otherGrid0 =  mAMRDataSet.getGridByIndex(ep0.refiningId());
    AxisAlignedBox<int> refiningCells0 = otherGrid0->getRefiningCells(ep0.grid(), ep0.idx(xAxis), ep0.idx(yAxis), ep0.idx(zAxis));

    int oIdx0[3];
    oIdx0[axis01] = qp0Axis01Type == min ? refiningCells0.min(axis01) : refiningCells0.max(axis01);
    oIdx0[axis03] = qp0Axis03Type == min ? refiningCells0.min(axis03) : refiningCells0.max(axis03);
    oIdx0[axisPerp0123] = quadFaceType == min ? refiningCells0.max(axisPerp0123) : refiningCells0.min(axisPerp0123);
    epRef0 = new AMRGridCell(otherGrid0, oIdx0);

    int oIdx1[3];
    oIdx1[axis01] = qp0Axis01Type == max ? refiningCells0.min(axis01) : refiningCells0.max(axis01);
    oIdx1[axis03] = qp0Axis03Type == min ? refiningCells0.min(axis03) : refiningCells0.max(axis03);
    oIdx1[axisPerp0123] = quadFaceType == min ? refiningCells0.max(axisPerp0123) : refiningCells0.min(axisPerp0123);
    epRef1 = new AMRGridCell(otherGrid0, oIdx1);
  }

  AMRGridCell *epRef2 = 0;
  AMRGridCell *epRef3 = 0;
  if (ep1.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 2;

    AMRGrid *otherGrid1 = mAMRDataSet.getGridByIndex(ep1.refiningId()); 
    AxisAlignedBox<int> refiningCells1 = otherGrid1->getRefiningCells(ep1.grid(), ep1.idx(xAxis), ep1.idx(yAxis), ep1.idx(zAxis));

    int oIdx2[3];
    oIdx2[axis01] = qp0Axis01Type == max ? refiningCells1.min(axis01) : refiningCells1.max(axis01);
    oIdx2[axis03] = qp0Axis03Type == max ? refiningCells1.min(axis03) : refiningCells1.max(axis03);
    oIdx2[axisPerp0123] = quadFaceType == min ? refiningCells1.max(axisPerp0123) : refiningCells1.min(axisPerp0123);
    epRef2 = new AMRGridCell(otherGrid1, oIdx2);

    int oIdx3[3];
    oIdx3[axis01] = qp0Axis01Type == min ? refiningCells1.min(axis01) : refiningCells1.max(axis01);
    oIdx3[axis03] = qp0Axis03Type == max ? refiningCells1.min(axis03) : refiningCells1.max(axis03);
    oIdx3[axisPerp0123] = quadFaceType == min ? refiningCells1.max(axisPerp0123) : refiningCells1.min(axisPerp0123);
    epRef3 = new AMRGridCell(otherGrid1, oIdx3);
  }

  switch (connectCase) {
    case 0:
      mStitchCellHandler.handleTrianglePrismCell(qp0, qp1, qp2, qp3, ep0, ep1);
      break;
    case 1:
      assert(epRef0->refiningId() == AMRGrid::NoGrid &&
	  epRef1->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= ep0.refiningId()) {
	mStitchCellHandler.handlePyramidCell(*epRef0, *epRef1, qp2, qp3, ep1);
	mStitchCellHandler.handleTrianglePrismCell(qp3, qp0, qp1, qp2, *epRef0, *epRef1);
      }
      break;
    case 2:
      assert(epRef2->refiningId() == AMRGrid::NoGrid &&
	  epRef3->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= ep1.refiningId()) {
	mStitchCellHandler.handlePyramidCell(qp0, qp1, *epRef2, *epRef3, ep0);
	mStitchCellHandler.handleTrianglePrismCell(qp1, qp2, qp3, qp0, *epRef2, *epRef3);
      }
      break;
    case 3:
      if (!mConsiderNextLevel ||
	  epRef0->refiningId() == AMRGrid::NoGrid &&
	  epRef1->refiningId() == AMRGrid::NoGrid &&
	  epRef2->refiningId() == AMRGrid::NoGrid &&
	  epRef3->refiningId() == AMRGrid::NoGrid)
	if (mStitchGrid->index() <= ep0.refiningId() && mStitchGrid->index() <= ep1.refiningId()) {
	  mStitchCellHandler.handleHexahedronCell(qp0, qp1, qp2, qp3, *epRef0, *epRef1, *epRef2, *epRef3);
	}
      break;
    default:
      assert(false);
      break;
  }
  delete epRef0;
  delete epRef1;
  delete epRef2;
  delete epRef3;
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectQuadToQuad(const AMRGridCell& qp0, const AMRGridCell& qp1, const AMRGridCell& qp2, const AMRGridCell& qp3, const AMRGridCell& oqp0, const AMRGridCell& oqp1, const AMRGridCell& oqp2, const AMRGridCell& oqp3, AxisType axis01, AxisType axis03, AxisType axisPerp0123, MinMaxType qp0Axis01Type, MinMaxType qp0Axis03Type, MinMaxType quadFaceType)
{
  // If the coarse grid points are from different grids, i.e. the quad results from several 
  // coarse grids rather than a single one, only generate stitch cells if we are stitching
  // the grid with the smallest index.
  if (qp0.containingId() < mStitchGrid->index() || qp1.containingId() < mStitchGrid->index() ||
      qp2.containingId() < mStitchGrid->index() || qp3.containingId() < mStitchGrid->index()) return;


  int connectCase = 0;

  AMRGridCell *oqp0Ref = 0;
  if (oqp0.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 1;
    AMRGrid *otherGrid0 = mAMRDataSet.getGridByIndex(oqp0.refiningId());
    AxisAlignedBox<int> refiningCells0 = otherGrid0->getRefiningCells(oqp0.grid(), oqp0.idx(xAxis), oqp0.idx(yAxis), oqp0.idx(zAxis));
    int oIdx0[3];
    oIdx0[axis01] = qp0Axis01Type == max ? refiningCells0.max(axis01) : refiningCells0.min(axis01);
    oIdx0[axis03] = qp0Axis03Type == max ? refiningCells0.max(axis03) : refiningCells0.min(axis03);
    oIdx0[axisPerp0123] = quadFaceType == min ? refiningCells0.max(axisPerp0123) : refiningCells0.min(axisPerp0123);
    oqp0Ref = new AMRGridCell(otherGrid0, oIdx0);
  }

  AMRGridCell *oqp1Ref = 0;
  if (oqp1.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 2;
    AMRGrid *otherGrid1 = mAMRDataSet.getGridByIndex(oqp1.refiningId());
    AxisAlignedBox<int> refiningCells1 = otherGrid1->getRefiningCells(oqp1.grid(), oqp1.idx(xAxis), oqp1.idx(yAxis), oqp1.idx(zAxis));
    int oIdx1[3];
    oIdx1[axis01] = qp0Axis01Type == max ? refiningCells1.min(axis01) : refiningCells1.max(axis01);
    oIdx1[axis03] = qp0Axis03Type == max ? refiningCells1.max(axis03) : refiningCells1.min(axis03);
    oIdx1[axisPerp0123] = quadFaceType == min ? refiningCells1.max(axisPerp0123) : refiningCells1.min(axisPerp0123);
    oqp1Ref = new AMRGridCell(otherGrid1, oIdx1);
  }

  AMRGridCell *oqp2Ref = 0;
  if (oqp2.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 4;
    AMRGrid *otherGrid2 = mAMRDataSet.getGridByIndex(oqp2.refiningId());
    AxisAlignedBox<int> refiningCells2 = otherGrid2->getRefiningCells(oqp2.grid(), oqp2.idx(xAxis), oqp2.idx(yAxis), oqp2.idx(zAxis));
    int oIdx2[3];
    oIdx2[axis01] = qp0Axis01Type == max ? refiningCells2.min(axis01) : refiningCells2.max(axis01);
    oIdx2[axis03] = qp0Axis03Type == max ? refiningCells2.min(axis03) : refiningCells2.max(axis03);
    oIdx2[axisPerp0123] = quadFaceType == min ? refiningCells2.max(axisPerp0123) : refiningCells2.min(axisPerp0123);
    oqp2Ref = new AMRGridCell(otherGrid2, oIdx2);
  }

  AMRGridCell *oqp3Ref = 0;
  if (oqp3.refiningId() != AMRGrid::NoGrid) {
    connectCase |= 8;
    AMRGrid *otherGrid3 = mAMRDataSet.getGridByIndex(oqp3.refiningId());
    AxisAlignedBox<int> refiningCells3 = otherGrid3->getRefiningCells(oqp3.grid(), oqp3.idx(xAxis), oqp3.idx(yAxis), oqp3.idx(zAxis));
    int oIdx3[3];
    oIdx3[axis01] = qp0Axis01Type == max ? refiningCells3.max(axis01) : refiningCells3.min(axis01);
    oIdx3[axis03] = qp0Axis03Type == max ? refiningCells3.min(axis03) : refiningCells3.max(axis03);
    oIdx3[axisPerp0123] = quadFaceType == min ? refiningCells3.max(axisPerp0123) : refiningCells3.min(axisPerp0123);
    oqp3Ref = new AMRGridCell(otherGrid3, oIdx3);
  }

  switch (connectCase) {
    case 0:
      mStitchCellHandler.handleHexahedronCell(qp0, qp1, qp2, qp3, oqp0, oqp1, oqp2, oqp3);
      break;
    case 1:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp0.refiningId()) {
	mStitchCellHandler.handlePyramidCell(qp0, qp1, qp2, qp3, *oqp0Ref);
	mStitchCellHandler.handlePyramidCell(qp3, qp2, oqp2, oqp3, *oqp0Ref);
	mStitchCellHandler.handlePyramidCell(qp2, qp1, oqp1, oqp2, *oqp0Ref);
      }
      break;
    case 2:
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp1.refiningId()) {
	mStitchCellHandler.handlePyramidCell(qp0, qp1, qp2, qp3, *oqp1Ref);
	mStitchCellHandler.handlePyramidCell(qp0, qp3, oqp3, oqp0, *oqp1Ref);
	mStitchCellHandler.handlePyramidCell(qp3, qp2, oqp2, oqp3, *oqp1Ref);
      }
      break;
    case 3:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp0.refiningId() && mStitchGrid->index() <= oqp1.refiningId()) {
	mStitchCellHandler.handleTrianglePrismCell(qp1, qp2, qp3, qp0, *oqp1Ref, *oqp0Ref);
	mStitchCellHandler.handleTrianglePrismCell(qp2, oqp2, oqp3, qp3,*oqp1Ref, *oqp0Ref);
      }
      break;
    case 4:
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp2.refiningId()) {
	mStitchCellHandler.handlePyramidCell(qp0, qp1, qp2, qp3, *oqp2Ref);
	mStitchCellHandler.handlePyramidCell(qp1, qp0, oqp0, oqp1, *oqp2Ref);
	mStitchCellHandler.handlePyramidCell(qp0, qp3, oqp3, oqp0, *oqp2Ref);
      }
      break;
    case 5:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp0.refiningId() && mStitchGrid->index() <= oqp2.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube024567(mStitchCellHandler, *oqp0Ref, oqp3, *oqp2Ref, oqp1, qp0, qp3, qp2, qp1);
      }
      break;
    case 6:
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp1.refiningId() && mStitchGrid->index() <= oqp2.refiningId()) {
	mStitchCellHandler.handleTrianglePrismCell(qp0, qp1, qp2, qp3, *oqp1Ref, *oqp2Ref);
	mStitchCellHandler.handleTrianglePrismCell(oqp0, qp0, qp3, oqp3, *oqp1Ref, *oqp2Ref);
      }
      break;
    case 7:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp2.refiningId() && mStitchGrid->index() <= oqp1.refiningId() && mStitchGrid->index() <= oqp0.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube1234567(mStitchCellHandler, oqp3, *oqp2Ref, *oqp1Ref, *oqp0Ref, qp3, qp2, qp1, qp0);
      }
      break;
    case 8:
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp3.refiningId()) {
	mStitchCellHandler.handlePyramidCell(qp0, qp1, qp2, qp3, *oqp3Ref);
	mStitchCellHandler.handlePyramidCell(qp1, qp0, oqp0, oqp1, *oqp3Ref);
	mStitchCellHandler.handlePyramidCell(qp2, qp1, oqp1, oqp2, *oqp3Ref);
      }
      break;
    case 9:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp0.refiningId() && mStitchGrid->index() <= oqp3.refiningId()) {
	mStitchCellHandler.handleTrianglePrismCell(qp0, qp1, qp2, qp3, *oqp0Ref, *oqp3Ref);
	mStitchCellHandler.handleTrianglePrismCell(qp1, oqp1, oqp2, qp2, *oqp0Ref, *oqp3Ref);
      }
      break;
    case 10:
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp1.refiningId() && mStitchGrid->index() <= oqp3.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube024567(mStitchCellHandler, *oqp1Ref, oqp0, *oqp3Ref, oqp2, qp1, qp0, qp3, qp2);
      }
      break;
    case 11:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp3.refiningId() && mStitchGrid->index() <= oqp1.refiningId() && mStitchGrid->index() <= oqp0.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube1234567(mStitchCellHandler, oqp2, *oqp1Ref, *oqp0Ref, *oqp3Ref, qp2, qp1, qp0, qp3);
      }
      break;
    case 12:
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp2.refiningId() && mStitchGrid->index() <= oqp3.refiningId()) {
	mStitchCellHandler.handleTrianglePrismCell(qp1, qp2, qp3, qp0, *oqp2Ref, *oqp3Ref);
	mStitchCellHandler.handleTrianglePrismCell(oqp1, qp1, qp0, oqp0, *oqp2Ref, *oqp3Ref);
      }
      break;
    case 13:
      assert(oqp0Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp3.refiningId() && mStitchGrid->index() <= oqp2.refiningId() && mStitchGrid->index() <= oqp0.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube1234567(mStitchCellHandler, oqp1, *oqp0Ref, *oqp3Ref, *oqp2Ref, qp1, qp0, qp3, qp2);
      }
      break;
    case 14:
      assert(oqp1Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp2Ref->refiningId() == AMRGrid::NoGrid);
      assert(oqp3Ref->refiningId() == AMRGrid::NoGrid);
      if (mStitchGrid->index() <= oqp3.refiningId() && mStitchGrid->index() <= oqp2.refiningId() && mStitchGrid->index() <= oqp1.refiningId()) {
	AMRGridStitcherAuxFct::generateCellsForDeformedCube1234567(mStitchCellHandler, oqp0, *oqp3Ref, *oqp2Ref, *oqp1Ref, qp0, qp3, qp2, qp1);
      }
      break;
    case 15:
      if (!mConsiderNextLevel ||
	  oqp0Ref->refiningId() == AMRGrid::NoGrid &&
	  oqp1Ref->refiningId() == AMRGrid::NoGrid &&
	  oqp2Ref->refiningId() == AMRGrid::NoGrid &&
	  oqp3Ref->refiningId() == AMRGrid::NoGrid) 
	if (mStitchGrid->index() <= oqp0.refiningId() && mStitchGrid->index() <= oqp1.refiningId() && mStitchGrid->index() <= oqp2.refiningId() && mStitchGrid->index() <= oqp3.refiningId()) {
	  mStitchCellHandler.handleHexahedronCell(qp0, qp1, qp2, qp3, *oqp0Ref, *oqp1Ref, *oqp2Ref, *oqp3Ref);
	}
      break;
    default:
      assert(false);
  }
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectEdgeSegmentToEdges(const AMRGridCell& ep0, const AMRGridCell& ep1, const AMRGridCell& cp0, const AMRGridCell& cp1, const AMRGridCell& cp2, AxisType edgeAxis, AxisType edgeAxis01, MinMaxType edgeAxis01Type, AxisType edgeAxis02, MinMaxType edgeAxis02Type)
{
  if (ep0.containingId() < mStitchGrid->index() || ep1.containingId() < mStitchGrid->index()) return;

  if (cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    mStitchCellHandler.handleTetrahedronCell(ep0, ep1, cp0, cp1);
    mStitchCellHandler.handleTetrahedronCell(ep0, cp0, ep1, cp2);
  }
  else if (cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0Ref0Idx[3];
      cp0Ref0Idx[edgeAxis] = cp0RefiningCells.min(edgeAxis);
      cp0Ref0Idx[edgeAxis01] = edgeAxis01Type == min ? cp0RefiningCells.max(edgeAxis01) : cp0RefiningCells.min(edgeAxis01);
      cp0Ref0Idx[edgeAxis02] = edgeAxis02Type == min ? cp0RefiningCells.max(edgeAxis01) : cp0RefiningCells.min(edgeAxis01);
      AMRGridCell cp0Ref0(cp0RefiningGrid, cp0Ref0Idx[0], cp0Ref0Idx[1], cp0Ref0Idx[2]);
      int cp0Ref1Idx[3];
      cp0Ref1Idx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0Ref1Idx[edgeAxis01] = cp0Ref0Idx[edgeAxis01];
      cp0Ref1Idx[edgeAxis02] = cp0Ref0Idx[edgeAxis02];
      AMRGridCell cp0Ref1(cp0RefiningGrid, cp0Ref1Idx[0], cp0Ref1Idx[1], cp0Ref1Idx[2]);
      mStitchCellHandler.handlePyramidCell(ep0, ep1, cp0Ref1, cp0Ref0, cp1);
      mStitchCellHandler.handlePyramidCell(ep0, cp0Ref0, cp0Ref1, ep1, cp2);
    }
  }
  else if (cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid) {
    AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
    AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
    int cp1Ref0Idx[3];
    cp1Ref0Idx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
    cp1Ref0Idx[edgeAxis01] = edgeAxis01Type == min ? cp1RefiningCells.min(edgeAxis01) : cp1RefiningCells.max(edgeAxis01);
    cp1Ref0Idx[edgeAxis02] = edgeAxis02Type == min ? cp1RefiningCells.max(edgeAxis02) : cp1RefiningCells.min(edgeAxis02);
    AMRGridCell cp1Ref0(cp1RefiningGrid, cp1Ref0Idx[0], cp1Ref0Idx[1], cp1Ref0Idx[2]);
    int cp1Ref1Idx[3];
    cp1Ref1Idx[edgeAxis] = cp1RefiningCells.max(edgeAxis);
    cp1Ref1Idx[edgeAxis01] = cp1Ref0Idx[edgeAxis01];
    cp1Ref1Idx[edgeAxis02] = cp1Ref0Idx[edgeAxis02];
    AMRGridCell cp1Ref1(cp1RefiningGrid, cp1Ref1Idx[0], cp1Ref1Idx[1], cp1Ref1Idx[2]);
    connectQuadToEdge(cp1Ref0, cp1Ref1, ep1, ep0, cp0, cp2, edgeAxis, edgeAxis02, edgeAxis01, min, edgeAxis02Type == min ? max : min, edgeAxis01Type);
  }
  else {
    assert(cp2.refiningId() != AMRGrid::NoGrid);
    AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
    AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
    int cp2Ref0Idx[3];
    cp2Ref0Idx[edgeAxis] = cp2RefiningCells.min(edgeAxis);
    cp2Ref0Idx[edgeAxis01] = edgeAxis01Type == min ? cp2RefiningCells.max(edgeAxis01) : cp2RefiningCells.min(edgeAxis01);
    cp2Ref0Idx[edgeAxis02] = edgeAxis02Type == min ? cp2RefiningCells.min(edgeAxis02) : cp2RefiningCells.max(edgeAxis02);
    AMRGridCell cp2Ref0(cp2RefiningGrid, cp2Ref0Idx[0], cp2Ref0Idx[1], cp2Ref0Idx[2]);
    int cp2Ref1Idx[3];
    cp2Ref1Idx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
    cp2Ref1Idx[edgeAxis01] = cp2Ref0Idx[edgeAxis01];
    cp2Ref1Idx[edgeAxis02] = cp2Ref0Idx[edgeAxis02];
    AMRGridCell cp2Ref1(cp2RefiningGrid, cp2Ref1Idx[0], cp2Ref1Idx[1], cp2Ref1Idx[2]);
    connectQuadToEdge(cp2Ref1, cp2Ref0, ep0, ep1, cp0, cp1, edgeAxis, edgeAxis01, edgeAxis02, max, edgeAxis01Type == min ? max : min, edgeAxis02Type);
  }
}

template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectEdgeSegmentToQuads(const AMRGridCell& ep0, const AMRGridCell& ep1, const AMRGridCell& cp0, const AMRGridCell& cp1, const AMRGridCell& cp2, const AMRGridCell& cp3, const AMRGridCell& cp4, const AMRGridCell& cp5, AxisType edgeAxis, AxisType addAxis0132, MinMaxType edgeAxis0132Type, AxisType addAxis0451, MinMaxType edgeAxis0451Type)
{
  if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 0
    mStitchCellHandler.handleTrianglePrismCell(cp2, cp0, cp1, cp3, ep0, ep1);
    mStitchCellHandler.handleTrianglePrismCell(cp0, cp4, cp5, cp1, ep0, ep1);
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 1
    if (mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism145(mStitchCellHandler, cp2, cp0Ref, cp1, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism045(mStitchCellHandler, cp0Ref, cp4, cp5, cp1, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 2
    if (mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism245(mStitchCellHandler, cp2, cp0, cp1Ref, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism345(mStitchCellHandler, cp0, cp4, cp5, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 3
    if (mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      mStitchCellHandler.handleTrianglePrismCell(cp0Ref, ep0, ep1, cp1Ref, cp2, cp3);
      mStitchCellHandler.handleTrianglePrismCell(ep0, cp0Ref, cp1Ref, ep1, cp4, cp5);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 4
    if (mStitchGrid->index() <= cp2.refiningId()) {
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp1, cp2Ref, cp3);
      std::cout << "Inserting vertex to avoid illegal tesselation (Case 4)." << std::endl;
      Vertex newPos = AMRGridStitcherAuxFct::centroid(ep0.pos(), ep1.pos(), cp0.pos(), cp1.pos(), cp2Ref.pos(), cp4.pos(), cp5.pos());
      double newVal = ep0.val() + ep1.val() + cp0.val() + cp1.val() + cp2Ref.val() + cp4.val() + cp5.val();
      newVal /= 7;
      mStitchCellHandler.handlePyramidCell(ep0, cp4, cp0, cp2Ref, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(ep0, ep1, cp5, cp4, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(cp0, cp4, cp5, cp1, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp1, cp5, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp2Ref, cp1, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(cp0, cp1, cp2Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp2Ref, ep1, newPos, newVal);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 5
    if (mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism0145(mStitchCellHandler, cp2Ref, cp0Ref, cp1, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism045(mStitchCellHandler, cp0Ref, cp4, cp5, cp1, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 6
    if (mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      mStitchCellHandler.handlePyramidCell(ep0, cp4, cp0, cp2Ref, cp1Ref);
      mStitchCellHandler.handlePyramidCell(ep0, ep1, cp5, cp4, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp3, ep1, ep0, cp2Ref, cp1Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 7
    if (mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01245(mStitchCellHandler, cp2Ref, cp0Ref, cp1Ref, cp3, ep0, ep1);
      mStitchCellHandler.handleTrianglePrismCell(cp0Ref, cp4, cp5, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 8
    if (mStitchGrid->index() <= cp3.refiningId()) {
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp3Ref, cp0, cp2);
      std::cout << "Inserting vertex to avoid illegal tesselation (Case 8)." << std::endl;
      Vertex newPos = AMRGridStitcherAuxFct::centroid(ep0.pos(), ep1.pos(), cp0.pos(), cp1.pos(), cp3Ref.pos(), cp4.pos(), cp5.pos());
      double newVal = ep0.val() + ep1.val() + cp0.val() + cp1.val() + cp3Ref.val() + cp4.val() + cp5.val();
      newVal /= 7;
      mStitchCellHandler.handlePyramidCell(cp0, cp4, cp5, cp1, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(ep0, ep1, cp5, cp4, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(ep1, cp3Ref, cp1, cp5, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp4, cp0, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp0, cp3Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(cp0, cp1, cp3Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp3Ref, ep1, newPos, newVal);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 9
    if (mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);

      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);

      mStitchCellHandler.handlePyramidCell(cp1, cp5, ep1, cp3Ref, cp0Ref);
      mStitchCellHandler.handlePyramidCell(ep0, ep1, cp5, cp4, cp0Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp2, cp3Ref, ep1, ep0, cp0Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 10
    if (mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);

      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism2345(mStitchCellHandler, cp2, cp0, cp1Ref, cp3Ref, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism345(mStitchCellHandler, cp0, cp4, cp5, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 11
    if (mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism12345(mStitchCellHandler, cp2, cp0Ref, cp1Ref, cp3Ref, ep0, ep1);
      mStitchCellHandler.handleTrianglePrismCell(cp0Ref, cp4, cp5, cp1Ref, ep0, ep1);
    }
  }
  else if (!(cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid) && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid) {
    // Case 12 ... 15, Case 28 ... 31, Case 44 ... 47
    AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
    AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
    int cp2RefIdx[3];
    cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
    cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
    cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
    AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
    AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
    AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
    int cp3RefIdx[3];
    cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
    cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
    cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
    AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
    connectQuadToQuad(ep0, cp2Ref, cp3Ref, ep1, cp4, cp0, cp1, cp5, addAxis0451, edgeAxis, addAxis0132, edgeAxis0451Type, max, edgeAxis0132Type);
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 16
    if (mStitchGrid->index() <= cp4.refiningId()) {
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp4Ref, cp1, cp5);
      std::cout << "Inserting vertex to avoid illegal tesselation (Case 16)." << std::endl;
      Vertex newPos = AMRGridStitcherAuxFct::centroid(ep0.pos(), ep1.pos(), cp0.pos(), cp1.pos(), cp2.pos(), cp3.pos(), cp4Ref.pos());
      double newVal = ep0.val() + ep1.val() + cp0.val() + cp1.val() + cp2.val() + cp3.val() + cp4Ref.val();
      newVal /= 7;
      mStitchCellHandler.handlePyramidCell(ep0, cp4Ref, cp0, cp2, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(cp0, cp1, cp3, cp2, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(ep0, cp2, cp3, ep1, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, ep1, cp4Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp1, cp4Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep1, cp3, cp1, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(cp0, cp4Ref, cp1, newPos, newVal);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 17
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism145(mStitchCellHandler, cp2, cp0Ref, cp1, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism0145(mStitchCellHandler, cp0Ref, cp4Ref, cp5, cp1, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 18
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);

      mStitchCellHandler.handlePyramidCell(ep0, cp2, cp3, ep1, cp1Ref);
      mStitchCellHandler.handlePyramidCell(ep0, cp4Ref, cp0, cp2, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp5, cp4Ref, ep0, ep1, cp1Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 19
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <=cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handleTrianglePrismCell(cp2, cp0Ref, cp1Ref, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01345(mStitchCellHandler, cp0Ref, cp4Ref, cp5, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 20
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp2.refiningId()) {
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp3, ep1, ep0, cp2Ref, cp4Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, cp2Ref, cp0, cp1, cp3, cp4Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, ep1, cp3, cp1, cp5, cp4Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 21
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handlePyramidCell(ep0, cp4Ref, cp0Ref, cp2Ref, ep1);
      mStitchCellHandler.handlePyramidCell(cp0Ref, cp4Ref, cp5, cp1, ep1);
      mStitchCellHandler.handlePyramidCell(cp0Ref, cp1, cp3, cp2Ref, ep1);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 22
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid024(mStitchCellHandler, ep1, cp3, cp1Ref, cp5, cp4Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid024(mStitchCellHandler, cp1Ref, cp3, cp2Ref, cp0, cp4Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp3, ep1, ep0, cp2Ref, cp4Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 23
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handlePyramidCell(cp0Ref, cp2Ref, ep0, cp4Ref, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp3, ep1, ep0, cp2Ref, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp5, cp4Ref, ep0, ep1, cp1Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 24
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index()) {
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handlePyramidCell(cp0, cp2, ep0, cp4Ref, cp3Ref);
      mStitchCellHandler.handleTetrahedronCell(cp0, cp4Ref, cp1, cp3Ref);
      mStitchCellHandler.handlePyramidCell(ep1, cp3Ref, cp1, cp5, cp4Ref);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp3Ref, ep1, cp4Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 25
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() && cp0.refiningId() <= mStitchGrid->index()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
      mStitchCellHandler.handleTrianglePrismCell(cp3Ref, cp1, cp5, ep1, cp0Ref, cp4Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism12345(mStitchCellHandler, cp2, cp3Ref, ep1, ep0, cp0Ref, cp4Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 26
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() && cp1.refiningId() <= mStitchGrid->index()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);

      mStitchCellHandler.handleTrianglePrismCell(cp4Ref, cp0, cp2, ep0, cp1Ref, cp3Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism12345(mStitchCellHandler, cp5, cp4Ref, ep0, ep1, cp1Ref, cp3Ref);
    }
  }
  else if (cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 27
    if (mStitchGrid->index() <= cp4.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() && cp1.refiningId() <= mStitchGrid->index() && cp0.refiningId() <= mStitchGrid->index()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
      AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
      int cp4RefIdx[3];
      cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
      cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
      cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
      AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism12345(mStitchCellHandler, cp2, cp0Ref, cp1Ref, cp3Ref, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01345(mStitchCellHandler, cp0Ref, cp4Ref, cp5, cp1Ref, ep0, ep1);
    }
  }
  // Case 28 ... 31 see Case 12 ... 15
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 32
    if (mStitchGrid->index() <= cp5.refiningId()) {
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp0, cp5Ref, cp4);
      std::cout << "Inserting vertex to avoid illegal tesselation (Case 32)." << std::endl;
      Vertex newPos = AMRGridStitcherAuxFct::centroid(ep0.pos(), ep1.pos(), cp0.pos(), cp1.pos(), cp2.pos(), cp3.pos(), cp5Ref.pos());
      double newVal = ep0.val() + ep1.val() + cp0.val() + cp1.val() + cp2.val() + cp3.val() + cp5Ref.val();
      newVal /= 7;
      mStitchCellHandler.handlePyramidCell(ep0, cp2, cp3, ep1, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(cp0, cp1, cp3, cp2, newPos, newVal);
      mStitchCellHandler.handlePyramidCell(ep1, cp3, cp1, cp5Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, ep1, cp5Ref, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp5Ref, cp0, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(ep0, cp0, cp2, newPos, newVal);
      mStitchCellHandler.handleTetrahedronCell(cp0, cp5Ref, cp1, newPos, newVal);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 33
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      mStitchCellHandler.handlePyramidCell(cp1, cp5Ref, ep1, cp3, cp0Ref);
      mStitchCellHandler.handlePyramidCell(ep0, cp2, cp3, ep1, cp0Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp4, ep0, ep1, cp5Ref, cp0Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 34
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism245(mStitchCellHandler, cp2, cp0, cp1Ref, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism2345(mStitchCellHandler, cp0, cp4, cp5Ref, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 35
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp1.refiningId() &&  mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      mStitchCellHandler.handleTrianglePrismCell(cp2, cp0Ref, cp1Ref, cp3, ep0, ep1);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism02345(mStitchCellHandler, cp0Ref, cp4, cp5Ref, cp1Ref, ep0, ep1);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 36
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp2.refiningId()) {
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism145(mStitchCellHandler, cp3, cp2Ref, cp0, cp1, ep1, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism0145(mStitchCellHandler, cp2Ref, ep0, cp4, cp0, ep1, cp5Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 37
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      mStitchCellHandler.handleTrianglePrismCell(cp3, cp2Ref, cp0Ref, cp1, ep1, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01345(mStitchCellHandler, cp2Ref, ep0, cp4, cp0Ref, ep1, cp5Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 38
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      mStitchCellHandler.handleTrianglePrismCell(ep0, cp4, cp0, cp2Ref, cp5Ref, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01245(mStitchCellHandler, ep1, ep0, cp2Ref, cp3, cp5Ref, cp1Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() != AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 39
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp2.refiningId() && mStitchGrid->index() <= cp1.refiningId() && cp0.refiningId() <= mStitchGrid->index()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp2RefiningGrid = mAMRDataSet.getGridByIndex(cp2.refiningId());
      AxisAlignedBox<int> cp2RefiningCells = cp2RefiningGrid->getRefiningCells(cp2.grid(), cp2.idx(xAxis), cp2.idx(yAxis), cp2.idx(zAxis));
      int cp2RefIdx[3];
      cp2RefIdx[edgeAxis] = cp2RefiningCells.max(edgeAxis);
      cp2RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp2RefiningCells.min(addAxis0132) : cp2RefiningCells.max(addAxis0132);
      cp2RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp2RefiningCells.max(addAxis0451) : cp2RefiningCells.min(addAxis0451);
      AMRGridCell cp2Ref(cp2RefiningGrid, cp2RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism02345(mStitchCellHandler, ep0, cp4, cp0Ref, cp2Ref, cp5Ref, cp1Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedTrianglePrism01245(mStitchCellHandler, ep1, ep0, cp2Ref, cp3, cp5Ref, cp1Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 40
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp3.refiningId()) {
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, ep0, cp4, cp0, cp2, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, cp3Ref, cp2, cp0, cp1, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp2, cp3Ref, ep1, ep0, cp5Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 41
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp3.refiningId() && cp0.refiningId() <= mStitchGrid->index()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);

      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid024(mStitchCellHandler, ep0, cp4, cp0Ref, cp2, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid024(mStitchCellHandler, cp3Ref, cp2, cp0Ref, cp1, cp5Ref);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp2, cp3Ref, ep1, ep0, cp5Ref);
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
    // Case 42
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() <= cp1.refiningId()) {
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      mStitchCellHandler.handlePyramidCell(ep1, cp3Ref, cp1Ref, cp5Ref, ep0);
      mStitchCellHandler.handlePyramidCell(cp0, cp4, cp5Ref, cp1Ref, ep0);
      mStitchCellHandler.handlePyramidCell(cp0, cp1Ref, cp3Ref, cp2, ep0); 
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() != AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() != AMRGrid::NoGrid && cp0.refiningId() != AMRGrid::NoGrid) {
    // Case 43
    if (mStitchGrid->index() <= cp5.refiningId() && mStitchGrid->index() <= cp3.refiningId() && mStitchGrid->index() <= cp1.refiningId() && mStitchGrid->index() <= cp0.refiningId()) {
      AMRGrid *cp0RefiningGrid = mAMRDataSet.getGridByIndex(cp0.refiningId());
      AxisAlignedBox<int> cp0RefiningCells = cp0RefiningGrid->getRefiningCells(cp0.grid(), cp0.idx(xAxis), cp0.idx(yAxis), cp0.idx(zAxis));
      int cp0RefIdx[3];
      cp0RefIdx[edgeAxis] = cp0RefiningCells.max(edgeAxis);
      cp0RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp0RefiningCells.max(addAxis0132) : cp0RefiningCells.min(addAxis0132);
      cp0RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp0RefiningCells.max(addAxis0451) : cp0RefiningCells.min(addAxis0451);
      AMRGridCell cp0Ref(cp0RefiningGrid, cp0RefIdx);
      AMRGrid *cp1RefiningGrid = mAMRDataSet.getGridByIndex(cp1.refiningId());
      AxisAlignedBox<int> cp1RefiningCells = cp1RefiningGrid->getRefiningCells(cp1.grid(), cp1.idx(xAxis), cp1.idx(yAxis), cp1.idx(zAxis));
      int cp1RefIdx[3];
      cp1RefIdx[edgeAxis] = cp1RefiningCells.min(edgeAxis);
      cp1RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp1RefiningCells.max(addAxis0132) : cp1RefiningCells.min(addAxis0132);
      cp1RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp1RefiningCells.max(addAxis0451) : cp1RefiningCells.min(addAxis0451);
      AMRGridCell cp1Ref(cp1RefiningGrid, cp1RefIdx);
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[edgeAxis] = cp3RefiningCells.min(edgeAxis);
      cp3RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp3RefiningCells.min(addAxis0132) : cp3RefiningCells.max(addAxis0132);
      cp3RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp3RefiningCells.max(addAxis0451) : cp3RefiningCells.min(addAxis0451);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
      cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
      cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      mStitchCellHandler.handlePyramidCell(ep1, cp3Ref, cp1Ref, cp5Ref, ep0);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp4, cp5Ref, cp1Ref, cp0Ref, ep0);
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, cp2, cp0Ref, cp1Ref, cp3Ref, ep0); 
    }
  }
  else if (cp5.refiningId() != AMRGrid::NoGrid && cp4.refiningId() != AMRGrid::NoGrid) {
    // Case 48 ... 63
    AMRGrid *cp4RefiningGrid = mAMRDataSet.getGridByIndex(cp4.refiningId());
    AxisAlignedBox<int> cp4RefiningCells = cp4RefiningGrid->getRefiningCells(cp4.grid(), cp4.idx(xAxis), cp4.idx(yAxis), cp4.idx(zAxis));
    int cp4RefIdx[3];
    cp4RefIdx[edgeAxis] = cp4RefiningCells.max(edgeAxis);
    cp4RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp4RefiningCells.max(addAxis0132) : cp4RefiningCells.min(addAxis0132);
    cp4RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp4RefiningCells.min(addAxis0451) : cp4RefiningCells.max(addAxis0451);
    AMRGridCell cp4Ref(cp4RefiningGrid, cp4RefIdx);
    AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
    AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
    int cp5RefIdx[3];
    cp5RefIdx[edgeAxis] = cp5RefiningCells.min(edgeAxis);
    cp5RefIdx[addAxis0132] = edgeAxis0132Type == min ? cp5RefiningCells.max(addAxis0132) : cp5RefiningCells.min(addAxis0132);
    cp5RefIdx[addAxis0451] = edgeAxis0451Type == min ? cp5RefiningCells.min(addAxis0451) : cp5RefiningCells.max(addAxis0451);
    AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
    connectQuadToQuad(ep0, ep1, cp5Ref, cp4Ref, cp2, cp3, cp1, cp0, edgeAxis, addAxis0132, addAxis0451, max, edgeAxis0132Type, edgeAxis0451Type);
  }
  else {
    assert(false);
  }
}

  template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::connectVertexToQuad(const AMRGridCell& vtx, const AMRGridCell& p0, const AMRGridCell& p1, const AMRGridCell& p2, const AMRGridCell& p3)
{
  assert(p0.refiningId() == AMRGrid::NoGrid);
  if (p1.refiningId() == AMRGrid::NoGrid && p2.refiningId() == AMRGrid::NoGrid && p3.refiningId() == AMRGrid::NoGrid) {
    mStitchCellHandler.handlePyramidCell(p0, p1, p2, p3, vtx);
  }
  else if (p1.refiningId() != AMRGrid::NoGrid && p2.refiningId() == AMRGrid::NoGrid && p3.refiningId() == AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p1.refiningId()) {
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, p1, p2, p3, p0, vtx);
    }
  }
  else if (p1.refiningId() == AMRGrid::NoGrid && p2.refiningId() != AMRGrid::NoGrid && p3.refiningId() == AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p2.refiningId()) {
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, p2, p3, p0, p1, vtx);
    }
  }
  else if (p1.refiningId() == AMRGrid::NoGrid && p2.refiningId() == AMRGrid::NoGrid && p3.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p3.refiningId()) {
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid04(mStitchCellHandler, p3, p0, p1, p2, vtx);
    }
  }
  else if (p1.refiningId() != AMRGrid::NoGrid && p2.refiningId() != AMRGrid::NoGrid && p3.refiningId() == AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p1.refiningId() && mStitchGrid->index() <= p2.refiningId()) {
      mStitchCellHandler.handlePyramidCell(p0, p1, p2, p3, vtx);
    }
  }
  else if (p1.refiningId() != AMRGrid::NoGrid && p2.refiningId() == AMRGrid::NoGrid && p3.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p1.refiningId() && mStitchGrid->index() <= p3.refiningId()) {
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid024(mStitchCellHandler, p1, p2, p3, p0, vtx);
    }
  }
  else if (p1.refiningId() == AMRGrid::NoGrid && p2.refiningId() != AMRGrid::NoGrid && p3.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p1.refiningId() && mStitchGrid->index() <= p2.refiningId()) {
      mStitchCellHandler.handlePyramidCell(p0, p1, p2, p3, vtx);
    }
  }
  else if (p1.refiningId() != AMRGrid::NoGrid && p2.refiningId() != AMRGrid::NoGrid && p3.refiningId() != AMRGrid::NoGrid) {
    if (mStitchGrid->index() <= p1.refiningId() && mStitchGrid->index() <= p2.refiningId() && mStitchGrid->index() <= p3.refiningId()) {
      AMRGridStitcherAuxFct::generateCellsForDeformedPyramid1234(mStitchCellHandler, p0, p1, p2, p3, vtx);
    }
  }
  else {
    assert(false);
  }

}

  template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::generateCellsForFace(AxisType axis0, AxisType axis1, AxisType axis2, MinMaxType faceType)
{
  assert(axis0 != axis1 && axis0 != axis2 && axis1 != axis2);

  try {
    // Here for k=0. k=refiningGrid->dims[2]-1 analogous
    // In this step we connect the faces of the refining grid with points in the
    // refined grid.
    for (int fFaceIdx0=0; fFaceIdx0<mStitchGrid->dims[axis0]-1; ++fFaceIdx0) 
      for (int fFaceIdx1=0; fFaceIdx1<mStitchGrid->dims[axis1]-1; ++fFaceIdx1) {
	int cIdx00[3];
	cIdx00[axis0] = mOriginInParentLevel[axis0] + fFaceIdx0/2;
	cIdx00[axis1] = mOriginInParentLevel[axis1] + fFaceIdx1/2;
	if (faceType == min)
	  cIdx00[axis2] = mOriginInParentLevel[axis2] - 1;
	else
	  cIdx00[axis2] = mMaxExtInParentLevel[axis2] + 1;
	AMRGridCell cp0 = mParentLevel->getGridCell(cIdx00);

	int cIdx01[3];
	cIdx01[axis0] = cIdx00[axis0] + 1;
	cIdx01[axis1] = cIdx00[axis1];
	cIdx01[axis2] = cIdx00[axis2];
	AMRGridCell cp3 = mParentLevel->getGridCell(cIdx01);

	int cIdx10[3];
	cIdx10[axis0] = cIdx00[axis0];
	cIdx10[axis1] = cIdx00[axis1] + 1;
	cIdx10[axis2] = cIdx00[axis2];
	AMRGridCell cp1 = mParentLevel->getGridCell(cIdx10);

	int cIdx11[3];
	cIdx11[axis0] = cIdx00[axis0] + 1;
	cIdx11[axis1] = cIdx00[axis1] + 1;
	cIdx11[axis2] = cIdx00[axis2];
	AMRGridCell cp2 = mParentLevel->getGridCell(cIdx11);;

	int fIdx00[3];
	fIdx00[axis0] = fFaceIdx0;
	fIdx00[axis1] = fFaceIdx1;
	if (faceType == min)
	  fIdx00[axis2] = 0;
	else
	  fIdx00[axis2] = mStitchGrid->dims[axis2] - 1;
	AMRGridCell fp0(mStitchGrid, fIdx00[0], fIdx00[1], fIdx00[2]);

	int fIdx01[3];
	fIdx01[axis0] = fIdx00[axis0] + 1;
	fIdx01[axis1] = fIdx00[axis1];
	fIdx01[axis2] = fIdx00[axis2];
	AMRGridCell fp3(mStitchGrid, fIdx01[0], fIdx01[1], fIdx01[2]);

	int fIdx10[3];
	fIdx10[axis0] = fIdx00[axis0];
	fIdx10[axis1] = fIdx00[axis1] + 1;
	fIdx10[axis2] = fIdx00[axis2];
	AMRGridCell fp1(mStitchGrid, fIdx10[0], fIdx10[1], fIdx10[2]);

	int fIdx11[3];
	fIdx11[axis0] = fIdx00[axis0] + 1;
	fIdx11[axis1] = fIdx00[axis1] + 1;
	fIdx11[axis2] = fIdx00[axis2];
	AMRGridCell fp2(mStitchGrid, fIdx11[0], fIdx11[1], fIdx11[2]);

	// Compute combination from 2D case. If i is even, we connect with the corresponding
	// point, if it is odd, we connect with the corresponding edge. The same for j.
	// The combination yields:
	// Axis 1  Axis 0  Result
	// ======================
	// point   point   point
	// point   edge    edge 0 (parallel axis 0)
	// edge    point   edge 1 (parallel to axis 1)
	// edge    edge 	 quad
	enum ConnectionTypeType { Point = 0, Edge0 = 1, Edge1 = 2, Quad = 3 };
	ConnectionTypeType connectionType = ConnectionTypeType((fFaceIdx1%2) << 1 | fFaceIdx0%2);
	switch (connectionType) {
	  case Point:
	    connectQuadToVertex(fp0, fp1, fp2, fp3, cp0, axis0, axis1, axis2, faceType);
	    break;
	  case Edge0:
	    //connectQuadToAxis0Edge(fp0, fp1, fp2, fp3, cp0, cp3, axis0, axis1, axis2, faceType);
	    connectQuadToEdge(fp0, fp1, fp2, fp3, cp0, cp3, axis1, axis0, axis2, min, max, faceType);
	    break;
	  case Edge1:
	    //connectQuadToAxis1Edge(fp0, fp1, fp2, fp3, cp0, cp1, axis0, axis1, axis2, faceType);
	    connectQuadToEdge(fp3, fp0, fp1, fp2, cp0, cp1, axis0, axis1, axis2, max, max, faceType);
	    break;
	  case Quad:
	    connectQuadToQuad(fp0, fp1, fp2, fp3, cp0, cp1, cp2, cp3, axis1, axis0, axis2, max, max, faceType);
	    break;
	  default:
	    std::cerr << "Internal Error (generateCellsForFace(AxisType, AxisType, AxisType)): Invalid enum value." << std::endl;
	    abort();
	}
      }
  }
  catch (AMRLevel::NoGrid ng) {
    std::cerr << "Error (AMRGridStitcher::generateCellsForFace(AxisType, AxisType, AxisType, MinMaxType): " << ng << std::endl;
  }
}

  template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::generateCellsForEdge(AxisType edgeAxis, AxisType nonEdgeAxis0, MinMaxType nonEdgeAxis0Type, AxisType nonEdgeAxis1, MinMaxType nonEdgeAxis1Type)
{
  AxisType addAxis0132;
  AxisType addAxis0451;
  MinMaxType addAxis0132Type;
  MinMaxType addAxis0451Type;
  if (nonEdgeAxis1Type == min && nonEdgeAxis0Type == min) {
    addAxis0132 = nonEdgeAxis0;
    addAxis0451 = nonEdgeAxis1;
    addAxis0132Type = min;
    addAxis0451Type = min;
  }
  else if (nonEdgeAxis1Type == min && nonEdgeAxis0Type == max) {
    addAxis0132 = nonEdgeAxis1;
    addAxis0451 = nonEdgeAxis0;
    addAxis0132Type = min;
    addAxis0451Type = max;
  }
  if (nonEdgeAxis1Type == max && nonEdgeAxis0Type == min) {
    addAxis0132 = nonEdgeAxis1;
    addAxis0451 = nonEdgeAxis0;
    addAxis0132Type = max;
    addAxis0451Type = min;
  }
  if (nonEdgeAxis1Type == max && nonEdgeAxis0Type == max) {
    addAxis0132 = nonEdgeAxis0;
    addAxis0451 = nonEdgeAxis1;
    addAxis0132Type = max;
    addAxis0451Type = max;
  }

  try {
    for (int edgeIdx=0; edgeIdx<mStitchGrid->dims[edgeAxis]-1; ++edgeIdx) {
      int cIdx0[3];
      cIdx0[edgeAxis] = mOriginInParentLevel[edgeAxis] + edgeIdx/2;
      if (nonEdgeAxis0Type == min) 
	cIdx0[nonEdgeAxis0] = mOriginInParentLevel[nonEdgeAxis0] - 1; 
      else
	cIdx0[nonEdgeAxis0] = mMaxExtInParentLevel[nonEdgeAxis0] + 1;
      if (nonEdgeAxis1Type == min) 
	cIdx0[nonEdgeAxis1] = mOriginInParentLevel[nonEdgeAxis1] - 1; 
      else
	cIdx0[nonEdgeAxis1] = mMaxExtInParentLevel[nonEdgeAxis1] + 1;
      AMRGridCell cp0 = mParentLevel->getGridCell(cIdx0);

      int cIdx1[3];
      cIdx1[edgeAxis] = cIdx0[edgeAxis] + 1;
      cIdx1[nonEdgeAxis0] = cIdx0[nonEdgeAxis0];
      cIdx1[nonEdgeAxis1] = cIdx0[nonEdgeAxis1];
      AMRGridCell cp1 = mParentLevel->getGridCell(cIdx1);

      int cIdx2[3];
      cIdx2[edgeAxis] = cIdx0[edgeAxis];
      if (addAxis0132Type == min)
	cIdx2[addAxis0132] = cIdx0[addAxis0132] + 1;
      else
	cIdx2[addAxis0132] = cIdx0[addAxis0132] - 1;
      cIdx2[addAxis0451] = cIdx0[addAxis0451];
      AMRGridCell cp2 = mParentLevel->getGridCell(cIdx2);

      int cIdx3[3];
      cIdx3[edgeAxis] = cIdx2[edgeAxis] + 1;
      cIdx3[addAxis0132] = cIdx2[addAxis0132];
      cIdx3[addAxis0451] = cIdx2[addAxis0451];
      AMRGridCell cp3 = mParentLevel->getGridCell(cIdx3);

      int cIdx4[3];
      cIdx4[edgeAxis] = cIdx0[edgeAxis];
      cIdx4[addAxis0132] = cIdx0[addAxis0132];
      if (addAxis0451Type == min)
	cIdx4[addAxis0451] = cIdx0[addAxis0451] + 1;
      else
	cIdx4[addAxis0451] = cIdx0[addAxis0451] - 1;
      AMRGridCell cp4 = mParentLevel->getGridCell(cIdx4);

      int cIdx5[3];
      cIdx5[edgeAxis] = cIdx4[edgeAxis] + 1;
      cIdx5[addAxis0132] = cIdx4[addAxis0132];
      cIdx5[addAxis0451] = cIdx4[addAxis0451];
      AMRGridCell cp5 = mParentLevel->getGridCell(cIdx5);

      int eIdx0[3];
      eIdx0[edgeAxis] = edgeIdx;
      if (nonEdgeAxis0Type == min)
	eIdx0[nonEdgeAxis0] = 0;
      else
	eIdx0[nonEdgeAxis0] = mStitchGrid->dims[nonEdgeAxis0]-1;
      if (nonEdgeAxis1Type == min)
	eIdx0[nonEdgeAxis1] = 0;
      else
	eIdx0[nonEdgeAxis1] = mStitchGrid->dims[nonEdgeAxis1]-1;
      AMRGridCell ep0(mStitchGrid, eIdx0[0], eIdx0[1], eIdx0[2]);

      int eIdx1[3];
      eIdx1[edgeAxis] = eIdx0[edgeAxis] + 1;
      eIdx1[nonEdgeAxis0] = eIdx0[nonEdgeAxis0];
      eIdx1[nonEdgeAxis1] = eIdx0[nonEdgeAxis1];
      AMRGridCell ep1(mStitchGrid, eIdx1[0], eIdx1[1], eIdx1[2]);

      if (edgeIdx % 2) {
	connectEdgeSegmentToQuads(ep0, ep1, cp0, cp1, cp2, cp3, cp4, cp5, edgeAxis, addAxis0132, addAxis0132Type, addAxis0451, addAxis0451Type);
      }
      else {
	connectEdgeSegmentToEdges(ep0, ep1, cp0, cp2, cp4, edgeAxis, addAxis0132, addAxis0132Type, addAxis0451, addAxis0451Type);
      }
    }
  }
  catch (AMRLevel::NoGrid ng) {
    std::cerr << "Error (AMRGridStitcher::generateCellsForEdge(AxisType, AxisType, MinMaxType, AxisType, MinMaxType): " << ng << std::endl;
  }
}

  template<class StitchCellHandler>
void AMRGridStitcher<StitchCellHandler>::generateCellsForVertex(AxisType axis0, AxisType axis1, AxisType axis2, MinMaxType axis0VertexType, MinMaxType axis1VertexType, MinMaxType axis2VertexType)
{
  try {
    int fIdx[3];
    fIdx[axis0] = axis0VertexType == min ? 0 : mStitchGrid->dims[axis0] - 1;
    fIdx[axis1] = axis1VertexType == min ? 0 : mStitchGrid->dims[axis1] - 1;
    fIdx[axis2] = axis2VertexType == min ? 0 : mStitchGrid->dims[axis2] - 1;
    AMRGridCell fp(mStitchGrid, fIdx);

    int cp0Idx[3];
    cp0Idx[axis0] = axis0VertexType == min ? mOriginInParentLevel[axis0] - 1 : mMaxExtInParentLevel[axis0] + 1;
    cp0Idx[axis1] = axis1VertexType == min ? mOriginInParentLevel[axis1] - 1 : mMaxExtInParentLevel[axis1] + 1;
    cp0Idx[axis2] = axis2VertexType == min ? mOriginInParentLevel[axis2] - 1 : mMaxExtInParentLevel[axis2] + 1;
    AMRGridCell cp0 = mParentLevel->getGridCell(cp0Idx);

    int cp1Idx[3];
    cp1Idx[axis0] = cp0Idx[axis0];
    cp1Idx[axis1] = axis1VertexType == min ? cp0Idx[axis1] + 1 : cp0Idx[axis1] - 1;
    cp1Idx[axis2] = cp0Idx[axis2];
    AMRGridCell cp1 = mParentLevel->getGridCell(cp1Idx);

    int cp2Idx[3];
    cp2Idx[axis0] = axis0VertexType == min ? cp0Idx[axis0] + 1 : cp0Idx[axis0] - 1;
    cp2Idx[axis1] = cp0Idx[axis1];
    cp2Idx[axis2] = cp0Idx[axis2];
    AMRGridCell cp2 = mParentLevel->getGridCell(cp2Idx);

    int cp3Idx[3];
    cp3Idx[axis0] = axis0VertexType == min ? cp1Idx[axis0] + 1 : cp1Idx[axis0] - 1;
    cp3Idx[axis1] = cp1Idx[axis1];
    cp3Idx[axis2] = cp1Idx[axis2];
    AMRGridCell cp3 = mParentLevel->getGridCell(cp3Idx);

    int cp4Idx[3];
    cp4Idx[axis0] = cp0Idx[axis0];
    cp4Idx[axis1] = cp0Idx[axis1];
    cp4Idx[axis2] = axis2VertexType == min ? cp0Idx[axis2] + 1 : cp0Idx[axis2] - 1;
    AMRGridCell cp4 = mParentLevel->getGridCell(cp4Idx);

    int cp5Idx[3];
    cp5Idx[axis0] = cp1Idx[axis0];
    cp5Idx[axis1] = cp1Idx[axis1];
    cp5Idx[axis2] = axis2VertexType == min ? cp1Idx[axis2] + 1 : cp1Idx[axis2] - 1;
    AMRGridCell cp5 = mParentLevel->getGridCell(cp5Idx);

    int cp6Idx[3];
    cp6Idx[axis0] = cp2Idx[axis0];
    cp6Idx[axis1] = cp0Idx[axis1];
    cp6Idx[axis2] = cp4Idx[axis2];
    AMRGridCell cp6 = mParentLevel->getGridCell(cp6Idx);

    if (cp6.refiningId() == AMRGrid::NoGrid && cp5.refiningId() == AMRGrid::NoGrid && cp4.refiningId() == AMRGrid::NoGrid && cp3.refiningId() == AMRGrid::NoGrid && cp2.refiningId() == AMRGrid::NoGrid && cp1.refiningId() == AMRGrid::NoGrid && cp0.refiningId() == AMRGrid::NoGrid) {
      mStitchCellHandler.handlePyramidCell(cp0, cp2, cp3, cp1, fp);
      mStitchCellHandler.handlePyramidCell(cp0, cp1, cp5, cp4, fp);
      mStitchCellHandler.handlePyramidCell(cp0, cp4, cp6, cp2, fp);
    }
    else if (!(cp3.refiningId() != AMRGrid::NoGrid || cp5.refiningId() != AMRGrid::NoGrid || cp6.refiningId() != AMRGrid::NoGrid)) {
      connectVertexToQuad(fp, cp3, cp1, cp0, cp2);
      connectVertexToQuad(fp, cp5, cp4, cp0, cp1);
      connectVertexToQuad(fp, cp6, cp2, cp0, cp4);
    }
    else if (cp3.refiningId() != AMRGrid::NoGrid) {
      AMRGrid *cp3RefiningGrid = mAMRDataSet.getGridByIndex(cp3.refiningId());
      AxisAlignedBox<int> cp3RefiningCells = cp3RefiningGrid->getRefiningCells(cp3.grid(), cp3.idx(xAxis), cp3.idx(yAxis), cp3.idx(zAxis));
      int cp3RefIdx[3];
      cp3RefIdx[axis0] = axis0VertexType == min ? cp3RefiningCells.min(axis0) : cp3RefiningCells.max(axis0);
      cp3RefIdx[axis1] = axis1VertexType == min ? cp3RefiningCells.min(axis1) : cp3RefiningCells.max(axis1);
      cp3RefIdx[axis2] = axis2VertexType == min ? cp3RefiningCells.max(axis2) : cp3RefiningCells.min(axis2);
      AMRGridCell cp3Ref(cp3RefiningGrid, cp3RefIdx);
      if (axis2VertexType == min)
	connectEdgeSegmentToQuads(cp3Ref, fp, cp0, cp4, cp2, cp6, cp1, cp5, axis2, axis0, axis0VertexType, axis1, axis1VertexType);
      else
	connectEdgeSegmentToQuads(fp, cp3Ref, cp4, cp0, cp5, cp1, cp6, cp2, axis2, axis1, axis1VertexType, axis0, axis0VertexType);
    }
    else if (cp5.refiningId() != AMRGrid::NoGrid) {
      AMRGrid *cp5RefiningGrid = mAMRDataSet.getGridByIndex(cp5.refiningId());
      AxisAlignedBox<int> cp5RefiningCells = cp5RefiningGrid->getRefiningCells(cp5.grid(), cp5.idx(xAxis), cp5.idx(yAxis), cp5.idx(zAxis));
      int cp5RefIdx[3];
      cp5RefIdx[axis0] = axis0VertexType == min ? cp5RefiningCells.max(axis0) : cp5RefiningCells.min(axis0);
      cp5RefIdx[axis1] = axis1VertexType == min ? cp5RefiningCells.min(axis1) : cp5RefiningCells.max(axis1);
      cp5RefIdx[axis2] = axis2VertexType == min ? cp5RefiningCells.min(axis2) : cp5RefiningCells.max(axis2);
      AMRGridCell cp5Ref(cp5RefiningGrid, cp5RefIdx);
      if (axis0VertexType == min)
	connectEdgeSegmentToQuads(cp5Ref, fp, cp0, cp2, cp1, cp3, cp4, cp6, axis0, axis1, axis1VertexType, axis2, axis2VertexType);
      else
	connectEdgeSegmentToQuads(fp, cp5Ref, cp2, cp0, cp6, cp4, cp3, cp1, axis0, axis2, axis2VertexType, axis1, axis1VertexType);
    }
    else if (cp6.refiningId() != AMRGrid::NoGrid) {
      AMRGrid *cp6RefiningGrid = mAMRDataSet.getGridByIndex(cp6.refiningId());
      AxisAlignedBox<int> cp6RefiningCells = cp6RefiningGrid->getRefiningCells(cp6.grid(), cp6.idx(xAxis), cp6.idx(yAxis), cp6.idx(zAxis));
      int cp6RefIdx[3];
      cp6RefIdx[axis0] = axis0VertexType == min ? cp6RefiningCells.min(axis0) : cp6RefiningCells.max(axis0);
      cp6RefIdx[axis1] = axis1VertexType == min ? cp6RefiningCells.max(axis1) : cp6RefiningCells.min(axis1);
      cp6RefIdx[axis2] = axis2VertexType == min ? cp6RefiningCells.min(axis2) : cp6RefiningCells.max(axis2);
      AMRGridCell cp6Ref(cp6RefiningGrid, cp6RefIdx);
      if (axis0VertexType == min)
	connectEdgeSegmentToQuads(cp6Ref, fp, cp0, cp1, cp4, cp5, cp2, cp3, axis1, axis2, axis2VertexType, axis0, axis0VertexType);
      else
	connectEdgeSegmentToQuads(fp, cp6Ref, cp1, cp0, cp3, cp2, cp5, cp4, axis1, axis0, axis0VertexType, axis2, axis2VertexType);
    }
  }
  catch (AMRLevel::NoGrid ng) {
    std::cerr << "Error (AMRGridStitcher::generateCellsForVertex(AxisType, AxisType, AxisType, MinMaxType, MinMaxType, MinMaxType): " << ng << std::endl;
  }
}