text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    TestUnstructuredGridGeometryFilter.cxx
  
    Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
    All rights reserved.
    See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
  
       This software is distributed WITHOUT ANY WARRANTY; without even
       the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
       PURPOSE.  See the above copyright notice for more information.
  
  =========================================================================*/
  // This test covers the vtkUnstructuredGridGeometryFilter class on all
  // types of cells. 
  // The command line arguments are:
  // -I        => run in interactive mode; unless this is used, the program will
  //              not allow interaction and exit
  // -D <path> => path to the data; the data should be in <path>/Data/
  
  
  // If READ_FILE is defined, the unstructured grid is read from a file otherwise
  // it is created cell by cell.
  //#define READ_FILE
  
  // If WRITE_RESULT is defined, the result of the surface filter is saved.
  //#define WRITE_RESULT
  
  // If USE_SHRINK is defined, each face is shrink to easily detect bad faces.
  #define USE_SHRINK
  
  // If FAST_GEOMETRY is defined, a vtkDataSetSurfaceFilter is used instead
  // of a vtkGeometryFilter at the end of the pipeline
  //#define FAST_GEOMETRY
  
  // If USE_CULLING is defined, backface culling is used to detect any bad
  // ordering of points defining a face.
  #define USE_CULLING
  
  // The configuration for the regression test is:
  // USE_SHRINK and USE_CULLING
  
  #include "vtkActor.h"
  #include "vtkCellData.h"
  #include "vtkTestUtilities.h"
  #include "vtkRegressionTestImage.h"
  #include "vtkRenderer.h"
  #include "vtkRenderWindow.h"
  #include "vtkRenderWindowInteractor.h"
  #include "vtkXMLUnstructuredGridReader.h"
  #include <assert.h>
  #include "vtkLookupTable.h"
  #include "vtkPolyData.h"
  #include "vtkUnstructuredGridGeometryFilter.h"
  #ifdef FAST_GEOMETRY
  # include "vtkDataSetSurfaceFilter.h"
  #else
  # include "vtkGeometryFilter.h"
  #endif
  #include "vtkPolyDataMapper.h"
  #include "vtkTimerLog.h"
  #include "vtkCamera.h"
  #include "vtkPointData.h"
  #include "vtkProperty.h"
  
  #ifndef READ_FILE
  # include "vtkUnstructuredGrid.h"
  # include "vtkFloatArray.h"
  # include "vtkIdTypeArray.h"
  # include "vtkVertex.h"
  # include "vtkPolyVertex.h"
  # include "vtkLine.h"
  # include "vtkPolyLine.h"
  # include "vtkQuadraticEdge.h"
  # include "vtkTriangle.h"
  # include "vtkTriangleStrip.h"
  # include "vtkPolygon.h"
  # include "vtkPixel.h"
  # include "vtkQuad.h"
  # include "vtkQuadraticTriangle.h"
  # include "vtkQuadraticQuad.h"
  # include "vtkBiQuadraticQuad.h"
  # include "vtkQuadraticLinearQuad.h"
  # include "vtkTetra.h"
  # include "vtkVoxel.h"
  # include "vtkHexahedron.h"
  # include "vtkWedge.h"
  # include "vtkPyramid.h"
  # include "vtkPentagonalPrism.h"
  # include "vtkHexagonalPrism.h"
  # include "vtkQuadraticTetra.h"
  # include "vtkQuadraticHexahedron.h"
  # include "vtkQuadraticWedge.h"
  # include "vtkQuadraticPyramid.h"
  # include "vtkTriQuadraticHexahedron.h"
  # include "vtkQuadraticLinearWedge.h"
  # include "vtkBiQuadraticQuadraticWedge.h"
  # include "vtkBiQuadraticQuadraticHexahedron.h"
  #endif
  
  #ifdef WRITE_RESULT
  # include "vtkXMLUnstructuredGridWriter.h"
  #endif // #ifdef WRITE_RESULT
  
  #ifdef USE_SHRINK
  # include "vtkShrinkFilter.h"
  #endif
  
  int TestUnstructuredGridGeometryFilter(int argc, char* argv[])
  {
    // Standard rendering classes
    vtkRenderer *renderer = vtkRenderer::New();
    vtkRenderWindow *renWin = vtkRenderWindow::New();
    renWin->AddRenderer(renderer);
    vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
    iren->SetRenderWindow(renWin);
    
  #ifdef READ_FILE
    // Load the mesh geometry and data from a file
    vtkXMLUnstructuredGridReader *reader = vtkXMLUnstructuredGridReader::New();
    char *cfname=vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/quadraticTetra01.vtu");
  
    reader->SetFileName( cfname );
    delete[] cfname;
  
    // Force reading
    reader->Update();
  #else
    // Create an unstructured grid.
    vtkUnstructuredGrid *grid=vtkUnstructuredGrid::New();
    vtkPoints *points=vtkPoints::New();
    vtkFloatArray *scalars=vtkFloatArray::New();
    scalars->SetName("ramp");
    vtkIdTypeArray *cellIds=vtkIdTypeArray::New();
    cellIds->SetName("cellIds");
    
    float scalar=0.0;
    const float scalarStep=0.1;
    vtkIdType cellId=0;
    double xOffset=0.0;
    double yOffset=0.0;
    vtkIdType pointId=0;
    
    // About 60 cells.
    grid->Allocate(60,60);
    
    // 0D: vertex
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkVertex *vertex=vtkVertex::New();
    vertex->GetPointIds()->SetId(0,pointId);
    ++pointId;
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(vertex->GetCellType(),vertex->GetPointIds());
    vertex->Delete();
    
    // 0D: polyvertex
    xOffset+=1.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkPolyVertex *polyVertex=vtkPolyVertex::New();
    polyVertex->GetPointIds()->SetNumberOfIds(2);
    polyVertex->GetPointIds()->SetId(0,pointId);
    ++pointId;
    polyVertex->GetPointIds()->SetId(1,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polyVertex->GetCellType(),polyVertex->GetPointIds());
    polyVertex->Delete();
    
    // 1D: line, polyline and quadratic edge
    yOffset+=2.0;
  //  xOffset+=1.0;
    xOffset=0.0;
    // 1D: line
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkLine *line=vtkLine::New();
    line->GetPointIds()->SetId(0,pointId);
    ++pointId;
    line->GetPointIds()->SetId(1,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(line->GetCellType(),line->GetPointIds());
    line->Delete();
    
    // 1D: polyline
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkPolyLine *polyLine=vtkPolyLine::New();
    polyLine->GetPointIds()->SetNumberOfIds(3);
    polyLine->GetPointIds()->SetId(0,pointId);
    ++pointId;
    polyLine->GetPointIds()->SetId(1,pointId);
    ++pointId;
    polyLine->GetPointIds()->SetId(2,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polyLine->GetCellType(),polyLine->GetPointIds());
    polyLine->Delete();
    
    // 1D: quadratic edge
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkQuadraticEdge *quadEdge=vtkQuadraticEdge::New();
    quadEdge->GetPointIds()->SetId(0,pointId);
    ++pointId;
    quadEdge->GetPointIds()->SetId(1,pointId);
    ++pointId;
    quadEdge->GetPointIds()->SetId(2,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadEdge->GetCellType(),quadEdge->GetPointIds());
    quadEdge->Delete();
    
    
    // 2D: triangle, triangle strip, polygon (triangle, quad, pentagon, hexagon),
    // pixel, quad.
    // quadratic quad, biquadratic quad, quadratic linear quad
    
    // 2D: triangle
    yOffset+=3.0;
    xOffset=0.0;
  //  xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkTriangle *triangle=vtkTriangle::New();
    triangle->GetPointIds()->SetId(0,pointId);
    ++pointId;
    triangle->GetPointIds()->SetId(1,pointId);
    ++pointId;
    triangle->GetPointIds()->SetId(2,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(triangle->GetCellType(),triangle->GetPointIds());
    triangle->Delete();
    
    // 2D: triangle strip
    // vertices are placed like that (first triangle is 0-1-2):
    // 0 2 4
    // 1 3 5
    // Be careful: the figure in the text book shows:
    // 1 3 5
    // 0 2 4
    // which create normals in the wrong way.
    
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,-0.2);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
     
    points->InsertNextPoint(xOffset+2.0,yOffset+1.0,-2.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
   
    vtkTriangleStrip *triangleStrip=vtkTriangleStrip::New();
    triangleStrip->GetPointIds()->SetNumberOfIds(6);
    int i=0;
    while(i<6)
      {
      triangleStrip->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(triangleStrip->GetCellType(),
                         triangleStrip->GetPointIds());
    triangleStrip->Delete();
    
    // 2D: polygon-triangle
    xOffset+=3.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkPolygon *polygon=vtkPolygon::New();
    polygon->GetPointIds()->SetNumberOfIds(3);
    polygon->GetPointIds()->SetId(0,pointId);
    ++pointId;
    polygon->GetPointIds()->SetId(1,pointId);
    ++pointId;
    polygon->GetPointIds()->SetId(2,pointId);
    ++pointId;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polygon->GetCellType(),polygon->GetPointIds());
    polygon->Delete();
    
    // 2D: polygon-quad
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.8,yOffset+0.8,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    polygon=vtkPolygon::New();
    polygon->GetPointIds()->SetNumberOfIds(4);
    i=0;
    while(i<4)
      {
      polygon->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polygon->GetCellType(),polygon->GetPointIds());
    polygon->Delete();
    
    // 2D: polygon-pentagon
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.2,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    polygon=vtkPolygon::New();
    polygon->GetPointIds()->SetNumberOfIds(5);
    i=0;
    while(i<5)
      {
      polygon->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polygon->GetCellType(),polygon->GetPointIds());
    polygon->Delete();
    
    // 2D: polygon-hexagon
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.2,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.1,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    polygon=vtkPolygon::New();
    polygon->GetPointIds()->SetNumberOfIds(6);
    i=0;
    while(i<6)
      {
      polygon->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(polygon->GetCellType(),polygon->GetPointIds());
    polygon->Delete();
    
    // 2D: pixel
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkPixel *pixel=vtkPixel::New();
    i=0;
    while(i<4)
      {
      pixel->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
  
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pixel->GetCellType(),pixel->GetPointIds());
    pixel->Delete();
    
     // 2D: quad
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.8,yOffset+0.8,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    vtkQuad *quad=vtkQuad::New();
    i=0;
    while(i<4)
      {
      quad->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quad->GetCellType(),quad->GetPointIds());
    quad->Delete();
    
    // 2D: quadratic triangle
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.3,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.3,yOffset+1.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.2,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkQuadraticTriangle *quadraticTriangle=vtkQuadraticTriangle::New();
    i=0;
    while(i<6)
      {
      quadraticTriangle->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticTriangle->GetCellType(),
                         quadraticTriangle->GetPointIds());
    quadraticTriangle->Delete();
    
    // 2D: quadratic quad
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.8,yOffset+0.8,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset-0.2,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.2,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+0.7,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.3,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkQuadraticQuad *quadraticQuad=vtkQuadraticQuad::New();
    i=0;
    while(i<8)
      {
      quadraticQuad->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticQuad->GetCellType(),
                         quadraticQuad->GetPointIds());
    quadraticQuad->Delete();
    
    // 2D: biquadratic quad: add a center point  vtkBiQuadraticQuad.h
    
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.8,yOffset+0.8,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset-0.2,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.2,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+0.7,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.3,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.4,yOffset+0.4,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkBiQuadraticQuad *biquadraticQuad=vtkBiQuadraticQuad::New();
    i=0;
    while(i<9)
      {
      biquadraticQuad->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biquadraticQuad->GetCellType(),
                         biquadraticQuad->GetPointIds());
    biquadraticQuad->Delete();
    
    // 2D: quadratic linear quad, no center, no mid-edge on sides
    
    xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.8,yOffset+0.8,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset-0.2,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+0.7,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    vtkQuadraticLinearQuad *quadraticLinearQuad=vtkQuadraticLinearQuad::New();
    i=0;
    while(i<6)
      {
      quadraticLinearQuad->GetPointIds()->SetId(i,pointId);
      ++pointId;
      ++i;
      }
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticLinearQuad->GetCellType(),
                         quadraticLinearQuad->GetPointIds());
    quadraticLinearQuad->Delete();
    
    // 3D: tetra, voxel, hexahedron, wedge, pyramid, pentagonal prism,
    // hexagonal prism,
    // quadratic tetra, quadratic hexa, quadratic wedge, quadratic pyramid,
    // triquadratic hexa, quadratic linear wedge, biquadratic quadratic wedge,
    // biquadratic quadratic hexa
    
    // 3D: tetra: 2 tetra with one common face
    yOffset+=3.0;
    xOffset=0.0;
  //  xOffset+=2.0;
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.25,yOffset+0.3,-2.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkTetra *tetra=vtkTetra::New();
    tetra->GetPointIds()->SetId(0,pointId);
    tetra->GetPointIds()->SetId(1,pointId+1);
    tetra->GetPointIds()->SetId(2,pointId+2);
    tetra->GetPointIds()->SetId(3,pointId+3);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(tetra->GetCellType(),tetra->GetPointIds());
    tetra->Delete();
    
    tetra=vtkTetra::New();
    tetra->GetPointIds()->SetId(0,pointId);
    tetra->GetPointIds()->SetId(1,pointId+2);
    tetra->GetPointIds()->SetId(2,pointId+1);
    tetra->GetPointIds()->SetId(3,pointId+4);
    
    pointId+=5;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(tetra->GetCellType(),tetra->GetPointIds());
    tetra->Delete();
  
    // 3D: voxel: 2 voxels with one common face
    xOffset+=2.0;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkVoxel *voxel=vtkVoxel::New();
    voxel->GetPointIds()->SetId(0,pointId);
    voxel->GetPointIds()->SetId(1,pointId+1);
    voxel->GetPointIds()->SetId(2,pointId+2);
    voxel->GetPointIds()->SetId(3,pointId+3);
    voxel->GetPointIds()->SetId(4,pointId+4);
    voxel->GetPointIds()->SetId(5,pointId+5);
    voxel->GetPointIds()->SetId(6,pointId+6);
    voxel->GetPointIds()->SetId(7,pointId+7);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(voxel->GetCellType(),voxel->GetPointIds());
    voxel->Delete();
    
    voxel=vtkVoxel::New();
    voxel->GetPointIds()->SetId(0,pointId+4);
    voxel->GetPointIds()->SetId(1,pointId+5);
    voxel->GetPointIds()->SetId(2,pointId+6);
    voxel->GetPointIds()->SetId(3,pointId+7);
    voxel->GetPointIds()->SetId(4,pointId+8);
    voxel->GetPointIds()->SetId(5,pointId+9);
    voxel->GetPointIds()->SetId(6,pointId+10);
    voxel->GetPointIds()->SetId(7,pointId+11);
    
    pointId+=12;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(voxel->GetCellType(),voxel->GetPointIds());
    voxel->Delete();
  
    // 3D: hexahedron: 2 hexahedra with one common face
    xOffset+=2.0;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkHexahedron *hexahedron=vtkHexahedron::New();
    hexahedron->GetPointIds()->SetId(0,pointId);
    hexahedron->GetPointIds()->SetId(1,pointId+1);
    hexahedron->GetPointIds()->SetId(2,pointId+2);
    hexahedron->GetPointIds()->SetId(3,pointId+3);
    hexahedron->GetPointIds()->SetId(4,pointId+4);
    hexahedron->GetPointIds()->SetId(5,pointId+5);
    hexahedron->GetPointIds()->SetId(6,pointId+6);
    hexahedron->GetPointIds()->SetId(7,pointId+7);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(hexahedron->GetCellType(),hexahedron->GetPointIds());
    hexahedron->Delete();
    
    hexahedron=vtkHexahedron::New();
    hexahedron->GetPointIds()->SetId(0,pointId+4);
    hexahedron->GetPointIds()->SetId(1,pointId+5);
    hexahedron->GetPointIds()->SetId(2,pointId+6);
    hexahedron->GetPointIds()->SetId(3,pointId+7);
    hexahedron->GetPointIds()->SetId(4,pointId+8);
    hexahedron->GetPointIds()->SetId(5,pointId+9);
    hexahedron->GetPointIds()->SetId(6,pointId+10);
    hexahedron->GetPointIds()->SetId(7,pointId+11);
    
    pointId+=12;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(hexahedron->GetCellType(),hexahedron->GetPointIds());
    hexahedron->Delete();
    
    // 3D: wedge: 3 wedges, some share a quad face, some share a triangle face
    xOffset+=2.0;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+0.9,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    vtkWedge *wedge=vtkWedge::New();
    wedge->GetPointIds()->SetId(0,pointId);
    wedge->GetPointIds()->SetId(1,pointId+1);
    wedge->GetPointIds()->SetId(2,pointId+2);
    wedge->GetPointIds()->SetId(3,pointId+3);
    wedge->GetPointIds()->SetId(4,pointId+4);
    wedge->GetPointIds()->SetId(5,pointId+5);
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(wedge->GetCellType(),wedge->GetPointIds());
    wedge->Delete();
    
    wedge=vtkWedge::New();
    wedge->GetPointIds()->SetId(0,pointId+3);
    wedge->GetPointIds()->SetId(1,pointId+4);
    wedge->GetPointIds()->SetId(2,pointId+5);
    wedge->GetPointIds()->SetId(3,pointId+6);
    wedge->GetPointIds()->SetId(4,pointId+7);
    wedge->GetPointIds()->SetId(5,pointId+8);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(wedge->GetCellType(),wedge->GetPointIds());
    wedge->Delete();
    
    wedge=vtkWedge::New();
    wedge->GetPointIds()->SetId(0,pointId+2);
    wedge->GetPointIds()->SetId(1,pointId+1);
    wedge->GetPointIds()->SetId(2,pointId+9);
    wedge->GetPointIds()->SetId(3,pointId+5);
    wedge->GetPointIds()->SetId(4,pointId+4);
    wedge->GetPointIds()->SetId(5,pointId+10);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(wedge->GetCellType(),wedge->GetPointIds());
    wedge->Delete();
    
    pointId+=11;
    
    // 3D: pyramid: 3 pyramids, some share the base quad face,
    // some share a triangle face
    xOffset+=2.0;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.2);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,-1.2);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset-1.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-0.1);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-0.9);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkPyramid *pyramid=vtkPyramid::New();
    pyramid->GetPointIds()->SetId(0,pointId);
    pyramid->GetPointIds()->SetId(1,pointId+1);
    pyramid->GetPointIds()->SetId(2,pointId+2);
    pyramid->GetPointIds()->SetId(3,pointId+3);
    pyramid->GetPointIds()->SetId(4,pointId+4);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pyramid->GetCellType(),pyramid->GetPointIds());
    pyramid->Delete();
    
    pyramid=vtkPyramid::New();
    pyramid->GetPointIds()->SetId(0,pointId+3);
    pyramid->GetPointIds()->SetId(1,pointId+2);
    pyramid->GetPointIds()->SetId(2,pointId+1);
    pyramid->GetPointIds()->SetId(3,pointId);
    pyramid->GetPointIds()->SetId(4,pointId+5);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pyramid->GetCellType(),pyramid->GetPointIds());
    pyramid->Delete();
    
    pyramid=vtkPyramid::New();
    pyramid->GetPointIds()->SetId(0,pointId+1);
    pyramid->GetPointIds()->SetId(1,pointId+6);
    pyramid->GetPointIds()->SetId(2,pointId+7);
    pyramid->GetPointIds()->SetId(3,pointId+2);
    pyramid->GetPointIds()->SetId(4,pointId+4);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pyramid->GetCellType(),pyramid->GetPointIds());
    pyramid->Delete();
    
    pointId+=8;
    
    // 3D: pentagonal prism: a wedge with a pentagonal base
    // Becareful, base face ordering is different from wedge...
    
    xOffset+=4.0;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 15
    points->InsertNextPoint(xOffset+2.5,yOffset+0.75,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.5,yOffset+2.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+2.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18
    points->InsertNextPoint(xOffset+2.5,yOffset+0.75,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.5,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+2.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    vtkPentagonalPrism *pentagonalPrism=vtkPentagonalPrism::New();
    pentagonalPrism->GetPointIds()->SetId(0,pointId);
    pentagonalPrism->GetPointIds()->SetId(1,pointId+4);
    pentagonalPrism->GetPointIds()->SetId(2,pointId+3);
    pentagonalPrism->GetPointIds()->SetId(3,pointId+2);
    pentagonalPrism->GetPointIds()->SetId(4,pointId+1);
    pentagonalPrism->GetPointIds()->SetId(5,pointId+5);
    pentagonalPrism->GetPointIds()->SetId(6,pointId+9);
    pentagonalPrism->GetPointIds()->SetId(7,pointId+8);
    pentagonalPrism->GetPointIds()->SetId(8,pointId+7);
    pentagonalPrism->GetPointIds()->SetId(9,pointId+6);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pentagonalPrism->GetCellType(),
                         pentagonalPrism->GetPointIds());
    pentagonalPrism->Delete();
    
    pentagonalPrism=vtkPentagonalPrism::New();
    pentagonalPrism->GetPointIds()->SetId(0,pointId+5);
    pentagonalPrism->GetPointIds()->SetId(1,pointId+9);
    pentagonalPrism->GetPointIds()->SetId(2,pointId+8);
    pentagonalPrism->GetPointIds()->SetId(3,pointId+7);
    pentagonalPrism->GetPointIds()->SetId(4,pointId+6);
    pentagonalPrism->GetPointIds()->SetId(5,pointId+10);
    pentagonalPrism->GetPointIds()->SetId(6,pointId+14);
    pentagonalPrism->GetPointIds()->SetId(7,pointId+13);
    pentagonalPrism->GetPointIds()->SetId(8,pointId+12);
    pentagonalPrism->GetPointIds()->SetId(9,pointId+11);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pentagonalPrism->GetCellType(),
                         pentagonalPrism->GetPointIds());
    pentagonalPrism->Delete();
    
    pentagonalPrism=vtkPentagonalPrism::New();
    pentagonalPrism->GetPointIds()->SetId(0,pointId+2);
    pentagonalPrism->GetPointIds()->SetId(1,pointId+3);
    pentagonalPrism->GetPointIds()->SetId(2,pointId+17);
    pentagonalPrism->GetPointIds()->SetId(3,pointId+16);
    pentagonalPrism->GetPointIds()->SetId(4,pointId+15);
    pentagonalPrism->GetPointIds()->SetId(5,pointId+7);
    pentagonalPrism->GetPointIds()->SetId(6,pointId+8);
    pentagonalPrism->GetPointIds()->SetId(7,pointId+20);
    pentagonalPrism->GetPointIds()->SetId(8,pointId+19);
    pentagonalPrism->GetPointIds()->SetId(9,pointId+18);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(pentagonalPrism->GetCellType(),
                         pentagonalPrism->GetPointIds());
    pentagonalPrism->Delete();
    
    pointId+=21;
    
    // 3D: hexagonal prism: a wedge with an hexagonal base.
    // Becareful, base face ordering is different from wedge...
    
    xOffset+=4.0;
    
    // hexagon
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // hexagon
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // hexagon
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.5,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18
    points->InsertNextPoint(xOffset+2.5,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    points->InsertNextPoint(xOffset+3.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.5,yOffset+1.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+1.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 22
    points->InsertNextPoint(xOffset+2.5,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    points->InsertNextPoint(xOffset+3.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+2.5,yOffset+1.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.5,yOffset+1.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
  
    vtkHexagonalPrism *hexagonalPrism=vtkHexagonalPrism::New();
    hexagonalPrism->GetPointIds()->SetId(0,pointId);
    hexagonalPrism->GetPointIds()->SetId(1,pointId+5);
    hexagonalPrism->GetPointIds()->SetId(2,pointId+4);
    hexagonalPrism->GetPointIds()->SetId(3,pointId+3);
    hexagonalPrism->GetPointIds()->SetId(4,pointId+2);
    hexagonalPrism->GetPointIds()->SetId(5,pointId+1);
    
    hexagonalPrism->GetPointIds()->SetId(6,pointId+6);
    hexagonalPrism->GetPointIds()->SetId(7,pointId+11);
    hexagonalPrism->GetPointIds()->SetId(8,pointId+10);
    hexagonalPrism->GetPointIds()->SetId(9,pointId+9);
    hexagonalPrism->GetPointIds()->SetId(10,pointId+8);
    hexagonalPrism->GetPointIds()->SetId(11,pointId+7);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(hexagonalPrism->GetCellType(),
                         hexagonalPrism->GetPointIds());
    hexagonalPrism->Delete();
  
    hexagonalPrism=vtkHexagonalPrism::New();
    hexagonalPrism->GetPointIds()->SetId(0,pointId+6);
    hexagonalPrism->GetPointIds()->SetId(1,pointId+11);
    hexagonalPrism->GetPointIds()->SetId(2,pointId+10);
    hexagonalPrism->GetPointIds()->SetId(3,pointId+9);
    hexagonalPrism->GetPointIds()->SetId(4,pointId+8);
    hexagonalPrism->GetPointIds()->SetId(5,pointId+7);
    
    hexagonalPrism->GetPointIds()->SetId(6,pointId+12);
    hexagonalPrism->GetPointIds()->SetId(7,pointId+17);
    hexagonalPrism->GetPointIds()->SetId(8,pointId+16);
    hexagonalPrism->GetPointIds()->SetId(9,pointId+15);
    hexagonalPrism->GetPointIds()->SetId(10,pointId+14);
    hexagonalPrism->GetPointIds()->SetId(11,pointId+13);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(hexagonalPrism->GetCellType(),
                         hexagonalPrism->GetPointIds());
    hexagonalPrism->Delete();
    
    hexagonalPrism=vtkHexagonalPrism::New();
    hexagonalPrism->GetPointIds()->SetId(0,pointId+2);
    hexagonalPrism->GetPointIds()->SetId(1,pointId+3);
    hexagonalPrism->GetPointIds()->SetId(2,pointId+21);
    hexagonalPrism->GetPointIds()->SetId(3,pointId+20);
    hexagonalPrism->GetPointIds()->SetId(4,pointId+19);
    hexagonalPrism->GetPointIds()->SetId(5,pointId+18);
    
    hexagonalPrism->GetPointIds()->SetId(6,pointId+8);
    hexagonalPrism->GetPointIds()->SetId(7,pointId+9);
    hexagonalPrism->GetPointIds()->SetId(8,pointId+25);
    hexagonalPrism->GetPointIds()->SetId(9,pointId+24);
    hexagonalPrism->GetPointIds()->SetId(10,pointId+23);
    hexagonalPrism->GetPointIds()->SetId(11,pointId+22);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(hexagonalPrism->GetCellType(),
                         hexagonalPrism->GetPointIds());
    hexagonalPrism->Delete();
    
    pointId+=26;
    
    // 3D: quadratic tetra: 2 tetra with one common face
    yOffset+=3.0;
    xOffset=0.0;
  //  xOffset+=2.0;
    
    // corner points
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.25,yOffset+0.3,-2.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on common face
    points->InsertNextPoint(xOffset+0.25,yOffset+0.6,0.0); // y=0.5->0.6 (concave)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.3,yOffset+1.5,0.0); // x=0.25->0.3 (convex)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset-0.2,yOffset+1.0,0.0); //x=0.0->-0.2 (convex)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on other edges for the first tetra
    points->InsertNextPoint(xOffset+0.0,yOffset+0.3,0.5); //y=0.25->0.3 (concave)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.4,yOffset+0.75,0.5); //x=0.25->0.4 (convex)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.25,1.0); //z=0.5->1.0 (convex)
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on other edges for the second tetra
    
    points->InsertNextPoint(xOffset+0.125,yOffset+0.15,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.125,yOffset+1.15,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.375,yOffset+0.65,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    vtkQuadraticTetra *quadraticTetra=vtkQuadraticTetra::New();
    quadraticTetra->GetPointIds()->SetId(0,pointId);
    quadraticTetra->GetPointIds()->SetId(1,pointId+1);
    quadraticTetra->GetPointIds()->SetId(2,pointId+2);
    quadraticTetra->GetPointIds()->SetId(3,pointId+3);
    quadraticTetra->GetPointIds()->SetId(4,pointId+5);
    quadraticTetra->GetPointIds()->SetId(5,pointId+6);
    quadraticTetra->GetPointIds()->SetId(6,pointId+7);
    quadraticTetra->GetPointIds()->SetId(7,pointId+8);
    quadraticTetra->GetPointIds()->SetId(8,pointId+9);
    quadraticTetra->GetPointIds()->SetId(9,pointId+10);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticTetra->GetCellType(),
                         quadraticTetra->GetPointIds());
    quadraticTetra->Delete();
    
    quadraticTetra=vtkQuadraticTetra::New();
    quadraticTetra->GetPointIds()->SetId(0,pointId);
    quadraticTetra->GetPointIds()->SetId(1,pointId+2);
    quadraticTetra->GetPointIds()->SetId(2,pointId+1);
    quadraticTetra->GetPointIds()->SetId(3,pointId+4);
    quadraticTetra->GetPointIds()->SetId(4,pointId+7);
    quadraticTetra->GetPointIds()->SetId(5,pointId+6);
    quadraticTetra->GetPointIds()->SetId(6,pointId+5);
    quadraticTetra->GetPointIds()->SetId(7,pointId+11);
    quadraticTetra->GetPointIds()->SetId(8,pointId+12);
    quadraticTetra->GetPointIds()->SetId(9,pointId+13);
    
    pointId+=14;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticTetra->GetCellType(),
                         quadraticTetra->GetPointIds());
    quadraticTetra->Delete();
  
    // 3D: quadratic hexahedron: 2 with a common face
    xOffset+=2.0;
    
    // a face (back): 0-3
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (common): 4-7
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (front): 8-11
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the back face: 12-15
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the common face: 16-19
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // mid-points between the back and common face: 20-23
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    
    // mid-points on the front face: 24-27
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points between the common face and the front face: 28-31
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    vtkQuadraticHexahedron *quadraticHexahedron=vtkQuadraticHexahedron::New();
    quadraticHexahedron->GetPointIds()->SetId(0,pointId);
    quadraticHexahedron->GetPointIds()->SetId(1,pointId+1);
    quadraticHexahedron->GetPointIds()->SetId(2,pointId+2);
    quadraticHexahedron->GetPointIds()->SetId(3,pointId+3);
    
    quadraticHexahedron->GetPointIds()->SetId(4,pointId+4);
    quadraticHexahedron->GetPointIds()->SetId(5,pointId+5);
    quadraticHexahedron->GetPointIds()->SetId(6,pointId+6);
    quadraticHexahedron->GetPointIds()->SetId(7,pointId+7);
    
    quadraticHexahedron->GetPointIds()->SetId(8,pointId+12);
    quadraticHexahedron->GetPointIds()->SetId(9,pointId+13);
    quadraticHexahedron->GetPointIds()->SetId(10,pointId+14);
    quadraticHexahedron->GetPointIds()->SetId(11,pointId+15);
    
    quadraticHexahedron->GetPointIds()->SetId(12,pointId+16);
    quadraticHexahedron->GetPointIds()->SetId(13,pointId+17);
    quadraticHexahedron->GetPointIds()->SetId(14,pointId+18);
    quadraticHexahedron->GetPointIds()->SetId(15,pointId+19);
    
    quadraticHexahedron->GetPointIds()->SetId(16,pointId+20);
    quadraticHexahedron->GetPointIds()->SetId(17,pointId+21);
    quadraticHexahedron->GetPointIds()->SetId(18,pointId+22);
    quadraticHexahedron->GetPointIds()->SetId(19,pointId+23);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticHexahedron->GetCellType(),
                         quadraticHexahedron->GetPointIds());
    quadraticHexahedron->Delete();
    
    quadraticHexahedron=vtkQuadraticHexahedron::New();
    quadraticHexahedron->GetPointIds()->SetId(0,pointId+4);
    quadraticHexahedron->GetPointIds()->SetId(1,pointId+5);
    quadraticHexahedron->GetPointIds()->SetId(2,pointId+6);
    quadraticHexahedron->GetPointIds()->SetId(3,pointId+7);
    
    quadraticHexahedron->GetPointIds()->SetId(4,pointId+8);
    quadraticHexahedron->GetPointIds()->SetId(5,pointId+9);
    quadraticHexahedron->GetPointIds()->SetId(6,pointId+10);
    quadraticHexahedron->GetPointIds()->SetId(7,pointId+11);
    
    quadraticHexahedron->GetPointIds()->SetId(8,pointId+16);
    quadraticHexahedron->GetPointIds()->SetId(9,pointId+17);
    quadraticHexahedron->GetPointIds()->SetId(10,pointId+18);
    quadraticHexahedron->GetPointIds()->SetId(11,pointId+19);
    
    quadraticHexahedron->GetPointIds()->SetId(12,pointId+24);
    quadraticHexahedron->GetPointIds()->SetId(13,pointId+25);
    quadraticHexahedron->GetPointIds()->SetId(14,pointId+26);
    quadraticHexahedron->GetPointIds()->SetId(15,pointId+27);
    
    quadraticHexahedron->GetPointIds()->SetId(16,pointId+28);
    quadraticHexahedron->GetPointIds()->SetId(17,pointId+29);
    quadraticHexahedron->GetPointIds()->SetId(18,pointId+30);
    quadraticHexahedron->GetPointIds()->SetId(19,pointId+31);
    
    
    pointId+=32;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticHexahedron->GetCellType(),
                         quadraticHexahedron->GetPointIds());
    quadraticHexahedron->Delete();
  
    // 3D: quadratic wedge: 3 wedges, some share a quadratic quad face, some
    // share a quadratic triangle face
    xOffset+=2.0;
    
    // corner points
    // triangle face of the first wedge
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 2
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 4
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 5
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge
    // 6
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 8
    points->InsertNextPoint(xOffset+0.5,yOffset+0.9,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // remaining vertices of the third wedge
    // 9
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points
    
    // triangle face of the first wedge: id=11
    // 11
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 12
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 13
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face: id=14
    // 14
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 15
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 16
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the first wedge: id=17
    // 17
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 19
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge: id=20
    // 20
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 21
    points->InsertNextPoint(xOffset+0.7,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 22
    points->InsertNextPoint(xOffset+0.3,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the second wedge: id=23
    // 23
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 24
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 25
    points->InsertNextPoint(xOffset+0.5,yOffset+0.95,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // first triangle face of the third wedge: id=26
    // 1+9
    // 26
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 27
    // 2+9
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    // second triangle face of the third wedge: id=28
    // 28
    // 4+10
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 29
    // 5+10
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the third wedge: id=30
    // 30
    // 9+10
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    
    vtkQuadraticWedge *quadraticWedge=vtkQuadraticWedge::New();
    quadraticWedge->GetPointIds()->SetId(0,pointId);
    quadraticWedge->GetPointIds()->SetId(1,pointId+1);
    quadraticWedge->GetPointIds()->SetId(2,pointId+2);
    quadraticWedge->GetPointIds()->SetId(3,pointId+3);
    quadraticWedge->GetPointIds()->SetId(4,pointId+4);
    quadraticWedge->GetPointIds()->SetId(5,pointId+5);
    
    quadraticWedge->GetPointIds()->SetId(6,pointId+11);
    quadraticWedge->GetPointIds()->SetId(7,pointId+12);
    quadraticWedge->GetPointIds()->SetId(8,pointId+13);
    quadraticWedge->GetPointIds()->SetId(9,pointId+14);
    quadraticWedge->GetPointIds()->SetId(10,pointId+15);
    quadraticWedge->GetPointIds()->SetId(11,pointId+16);
    
    quadraticWedge->GetPointIds()->SetId(12,pointId+17);
    quadraticWedge->GetPointIds()->SetId(13,pointId+18);
    quadraticWedge->GetPointIds()->SetId(14,pointId+19);
    
    
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticWedge->GetCellType(),
                         quadraticWedge->GetPointIds());
    quadraticWedge->Delete();
  
    // this wedge shares a triangle face
    quadraticWedge=vtkQuadraticWedge::New();
    quadraticWedge->GetPointIds()->SetId(0,pointId+3);
    quadraticWedge->GetPointIds()->SetId(1,pointId+4);
    quadraticWedge->GetPointIds()->SetId(2,pointId+5);
    quadraticWedge->GetPointIds()->SetId(3,pointId+6);
    quadraticWedge->GetPointIds()->SetId(4,pointId+7);
    quadraticWedge->GetPointIds()->SetId(5,pointId+8);
    
    quadraticWedge->GetPointIds()->SetId(6,pointId+14);
    quadraticWedge->GetPointIds()->SetId(7,pointId+15);
    quadraticWedge->GetPointIds()->SetId(8,pointId+16);
    quadraticWedge->GetPointIds()->SetId(9,pointId+20);
    quadraticWedge->GetPointIds()->SetId(10,pointId+21);
    quadraticWedge->GetPointIds()->SetId(11,pointId+22);
    
    quadraticWedge->GetPointIds()->SetId(12,pointId+23);
    quadraticWedge->GetPointIds()->SetId(13,pointId+24);
    quadraticWedge->GetPointIds()->SetId(14,pointId+25);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticWedge->GetCellType(),
                         quadraticWedge->GetPointIds());
    quadraticWedge->Delete();
  
    // this wedge shares a quad face
    quadraticWedge=vtkQuadraticWedge::New();
    quadraticWedge->GetPointIds()->SetId(0,pointId+2);
    quadraticWedge->GetPointIds()->SetId(1,pointId+1);
    quadraticWedge->GetPointIds()->SetId(2,pointId+9);
    quadraticWedge->GetPointIds()->SetId(3,pointId+5);
    quadraticWedge->GetPointIds()->SetId(4,pointId+4);
    quadraticWedge->GetPointIds()->SetId(5,pointId+10);
    
    quadraticWedge->GetPointIds()->SetId(6,pointId+12);
    quadraticWedge->GetPointIds()->SetId(7,pointId+26);
    quadraticWedge->GetPointIds()->SetId(8,pointId+27);
    quadraticWedge->GetPointIds()->SetId(9,pointId+15);
    quadraticWedge->GetPointIds()->SetId(10,pointId+28);
    quadraticWedge->GetPointIds()->SetId(11,pointId+29);
    
    quadraticWedge->GetPointIds()->SetId(12,pointId+19);
    quadraticWedge->GetPointIds()->SetId(13,pointId+18);
    quadraticWedge->GetPointIds()->SetId(14,pointId+30);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticWedge->GetCellType(),
                         quadraticWedge->GetPointIds());
    quadraticWedge->Delete();
  
    pointId+=31;
    
     // 3D: quadratic pyramid: 3 pyramids, some share the base quad face,
    // some share a triangle face
    
    xOffset+=2.0;
    
    // quad face
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.2);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 2
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,-1.2);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // top vertex
    // 4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // bottom vertex
    // 5
    points->InsertNextPoint(xOffset+0.5,yOffset-1.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other quad base
    // 6
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-0.1);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-0.9);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points
    
    // 8=(0+1)/2
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.1);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 9=(1+2)/2
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10=(2+3)/2
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,-1.1);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 11=(0+3)/2
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    // 12=(0+4)/2
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,-0.15);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 13=(1+4)/2
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,-0.25);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 14=(2+4)/2
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,-0.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // 15=(3+4)/2
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,-0.85);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 16=(3+5)/2
    points->InsertNextPoint(xOffset+0.25,yOffset-0.5,-0.85);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 17=(2+5)/2
    points->InsertNextPoint(xOffset+0.75,yOffset-0.5,-0.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18=(1+5)/2
    points->InsertNextPoint(xOffset+0.75,yOffset-0.5,-0.25);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 19=(0+5)/2
    points->InsertNextPoint(xOffset+0.25,yOffset-0.5,-0.15);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 20=(1+6)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+0.0,-0.05);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 21=(6+7)/2
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 22=(2+7)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+0.0,-0.95);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 23=(4+6)/2
    points->InsertNextPoint(xOffset+1.25,yOffset+0.5,-0.3);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 24=(4+7)/2
    points->InsertNextPoint(xOffset+1.25,yOffset+0.5,-0.7);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    vtkQuadraticPyramid *quadraticPyramid=vtkQuadraticPyramid::New();
    quadraticPyramid->GetPointIds()->SetId(0,pointId);
    quadraticPyramid->GetPointIds()->SetId(1,pointId+1);
    quadraticPyramid->GetPointIds()->SetId(2,pointId+2);
    quadraticPyramid->GetPointIds()->SetId(3,pointId+3);
    quadraticPyramid->GetPointIds()->SetId(4,pointId+4);
    
    quadraticPyramid->GetPointIds()->SetId(5,pointId+8);
    quadraticPyramid->GetPointIds()->SetId(6,pointId+9);
    quadraticPyramid->GetPointIds()->SetId(7,pointId+10);
    quadraticPyramid->GetPointIds()->SetId(8,pointId+11);
    
    quadraticPyramid->GetPointIds()->SetId(9,pointId+12);
    quadraticPyramid->GetPointIds()->SetId(10,pointId+13);
    quadraticPyramid->GetPointIds()->SetId(11,pointId+14);
    quadraticPyramid->GetPointIds()->SetId(12,pointId+15);
    
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticPyramid->GetCellType(),
                         quadraticPyramid->GetPointIds());
    quadraticPyramid->Delete();
    
    quadraticPyramid=vtkQuadraticPyramid::New();
    quadraticPyramid->GetPointIds()->SetId(0,pointId+3);
    quadraticPyramid->GetPointIds()->SetId(1,pointId+2);
    quadraticPyramid->GetPointIds()->SetId(2,pointId+1);
    quadraticPyramid->GetPointIds()->SetId(3,pointId);
    quadraticPyramid->GetPointIds()->SetId(4,pointId+5);
    
    quadraticPyramid->GetPointIds()->SetId(5,pointId+10);
    quadraticPyramid->GetPointIds()->SetId(6,pointId+9);
    quadraticPyramid->GetPointIds()->SetId(7,pointId+8);
    quadraticPyramid->GetPointIds()->SetId(8,pointId+11);
    
    quadraticPyramid->GetPointIds()->SetId(9,pointId+16);
    quadraticPyramid->GetPointIds()->SetId(10,pointId+17);
    quadraticPyramid->GetPointIds()->SetId(11,pointId+18);
    quadraticPyramid->GetPointIds()->SetId(12,pointId+19);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticPyramid->GetCellType(),
                         quadraticPyramid->GetPointIds());
    quadraticPyramid->Delete();
    
    quadraticPyramid=vtkQuadraticPyramid::New();
    quadraticPyramid->GetPointIds()->SetId(0,pointId+1);
    quadraticPyramid->GetPointIds()->SetId(1,pointId+6);
    quadraticPyramid->GetPointIds()->SetId(2,pointId+7);
    quadraticPyramid->GetPointIds()->SetId(3,pointId+2);
    quadraticPyramid->GetPointIds()->SetId(4,pointId+4);
    
    quadraticPyramid->GetPointIds()->SetId(5,pointId+20);
    quadraticPyramid->GetPointIds()->SetId(6,pointId+21);
    quadraticPyramid->GetPointIds()->SetId(7,pointId+22);
    quadraticPyramid->GetPointIds()->SetId(8,pointId+9);
    
    quadraticPyramid->GetPointIds()->SetId(9,pointId+13);
    quadraticPyramid->GetPointIds()->SetId(10,pointId+23);
    quadraticPyramid->GetPointIds()->SetId(11,pointId+24);
    quadraticPyramid->GetPointIds()->SetId(12,pointId+14);
    
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticPyramid->GetCellType(),
                         quadraticPyramid->GetPointIds());
    quadraticPyramid->Delete();
    
    pointId+=25;
    
    // 3D: triquadratic hexahedron: 2 with a common face
    xOffset+=2.0;
    
    // a face (back): 0-3
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 2
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (common): 4-7
    // 4
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 5
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 6
    points->InsertNextPoint(xOffset+0.9,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+0.1,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (front): 8-11
    // 8
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 9
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 11
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the back face: 12-15
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the common face: 16-19
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // mid-points between the back and common face: 20-23
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    
    // mid-points on the front face: 24-27
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points between the common face and the front face: 28-31
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // face-centered points
    // 32=(0+1+4+5)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    // 33=(1+2+5+6)/4
    points->InsertNextPoint(xOffset+0.95,yOffset+1.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 34=(2+3+6+7)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 35=(0+3+4+7)/4
    points->InsertNextPoint(xOffset+0.05,yOffset+1.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 36=(0+1+2+3)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // 37=(4+5+6+7)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    
    // other hexa
    // 38=(4+5+8+9)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // 39=(5+6+9+10)/4
    points->InsertNextPoint(xOffset+0.95,yOffset+1.0,3.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 40=(6+7+10+11)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 41=(4+7+8+11)/4
    points->InsertNextPoint(xOffset+0.05,yOffset+1.0,3.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // 42=(8+9+10+11)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    vtkTriQuadraticHexahedron *triQuadraticHexahedron;
    triQuadraticHexahedron=vtkTriQuadraticHexahedron::New();
    triQuadraticHexahedron->GetPointIds()->SetId(0,pointId);
    triQuadraticHexahedron->GetPointIds()->SetId(1,pointId+1);
    triQuadraticHexahedron->GetPointIds()->SetId(2,pointId+2);
    triQuadraticHexahedron->GetPointIds()->SetId(3,pointId+3);
    
    triQuadraticHexahedron->GetPointIds()->SetId(4,pointId+4);
    triQuadraticHexahedron->GetPointIds()->SetId(5,pointId+5);
    triQuadraticHexahedron->GetPointIds()->SetId(6,pointId+6);
    triQuadraticHexahedron->GetPointIds()->SetId(7,pointId+7);
    
    triQuadraticHexahedron->GetPointIds()->SetId(8,pointId+12);
    triQuadraticHexahedron->GetPointIds()->SetId(9,pointId+13);
    triQuadraticHexahedron->GetPointIds()->SetId(10,pointId+14);
    triQuadraticHexahedron->GetPointIds()->SetId(11,pointId+15);
    
    triQuadraticHexahedron->GetPointIds()->SetId(12,pointId+16);
    triQuadraticHexahedron->GetPointIds()->SetId(13,pointId+17);
    triQuadraticHexahedron->GetPointIds()->SetId(14,pointId+18);
    triQuadraticHexahedron->GetPointIds()->SetId(15,pointId+19);
    
    triQuadraticHexahedron->GetPointIds()->SetId(16,pointId+20);
    triQuadraticHexahedron->GetPointIds()->SetId(17,pointId+21);
    triQuadraticHexahedron->GetPointIds()->SetId(18,pointId+22);
    triQuadraticHexahedron->GetPointIds()->SetId(19,pointId+23);
    
    // before: 32,33,34,35,36,37
    triQuadraticHexahedron->GetPointIds()->SetId(20,pointId+35);
    triQuadraticHexahedron->GetPointIds()->SetId(21,pointId+33);
    triQuadraticHexahedron->GetPointIds()->SetId(22,pointId+32);
    triQuadraticHexahedron->GetPointIds()->SetId(23,pointId+34);
    triQuadraticHexahedron->GetPointIds()->SetId(24,pointId+36);
    triQuadraticHexahedron->GetPointIds()->SetId(25,pointId+37);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(triQuadraticHexahedron->GetCellType(),
                         triQuadraticHexahedron->GetPointIds());
    triQuadraticHexahedron->Delete();
    
    triQuadraticHexahedron=vtkTriQuadraticHexahedron::New();
    triQuadraticHexahedron->GetPointIds()->SetId(0,pointId+4);
    triQuadraticHexahedron->GetPointIds()->SetId(1,pointId+5);
    triQuadraticHexahedron->GetPointIds()->SetId(2,pointId+6);
    triQuadraticHexahedron->GetPointIds()->SetId(3,pointId+7);
    
    triQuadraticHexahedron->GetPointIds()->SetId(4,pointId+8);
    triQuadraticHexahedron->GetPointIds()->SetId(5,pointId+9);
    triQuadraticHexahedron->GetPointIds()->SetId(6,pointId+10);
    triQuadraticHexahedron->GetPointIds()->SetId(7,pointId+11);
    
    triQuadraticHexahedron->GetPointIds()->SetId(8,pointId+16);
    triQuadraticHexahedron->GetPointIds()->SetId(9,pointId+17);
    triQuadraticHexahedron->GetPointIds()->SetId(10,pointId+18);
    triQuadraticHexahedron->GetPointIds()->SetId(11,pointId+19);
    
    triQuadraticHexahedron->GetPointIds()->SetId(12,pointId+24);
    triQuadraticHexahedron->GetPointIds()->SetId(13,pointId+25);
    triQuadraticHexahedron->GetPointIds()->SetId(14,pointId+26);
    triQuadraticHexahedron->GetPointIds()->SetId(15,pointId+27);
    
    triQuadraticHexahedron->GetPointIds()->SetId(16,pointId+28);
    triQuadraticHexahedron->GetPointIds()->SetId(17,pointId+29);
    triQuadraticHexahedron->GetPointIds()->SetId(18,pointId+30);
    triQuadraticHexahedron->GetPointIds()->SetId(19,pointId+31);
    
    // before: 38,39,40,41,37,42
    triQuadraticHexahedron->GetPointIds()->SetId(20,pointId+41);
    triQuadraticHexahedron->GetPointIds()->SetId(21,pointId+39);
    triQuadraticHexahedron->GetPointIds()->SetId(22,pointId+38);
    triQuadraticHexahedron->GetPointIds()->SetId(23,pointId+40);
    triQuadraticHexahedron->GetPointIds()->SetId(24,pointId+37);
    triQuadraticHexahedron->GetPointIds()->SetId(25,pointId+42);
   
    pointId+=43;
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(triQuadraticHexahedron->GetCellType(),
                         triQuadraticHexahedron->GetPointIds());
    triQuadraticHexahedron->Delete();
  
    // 3D: quadratic linear wedge: 3 wedges, some share a quadratic linear quad
    // face, some share a quadratic triangle face
    // NOTE: ordering is different from linear wedge or quadratic wedge
    // linear or quad: triangle 0-1-2 points outside, 3-4-5 points inside
    // here: 0-1-2 points inside, 3-4-5 points outside
    // 
    xOffset+=2.0;
    
    // corner points
    // triangle face of the first wedge
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 2
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 4
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 5
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge
    // 6
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 8
    points->InsertNextPoint(xOffset+0.5,yOffset+0.9,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // remaining vertices of the third wedge
    // 9
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points
    
    // triangle face of the first wedge: id=11
    // 11
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 12
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 13
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face: id=14
    // 14
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 15
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 16
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge: id=17
    // 17
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18
    points->InsertNextPoint(xOffset+0.7,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 19
    points->InsertNextPoint(xOffset+0.3,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // first triangle face of the third wedge: id=20
    // 1+9
    // 20
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 21
    // 2+9
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    // second triangle face of the third wedge: id=22
    // 22
    // 4+10
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 23
    // 5+10
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
      
    vtkQuadraticLinearWedge *quadraticLinearWedge=vtkQuadraticLinearWedge::New();
    quadraticLinearWedge->GetPointIds()->SetId(0,pointId);
    quadraticLinearWedge->GetPointIds()->SetId(1,pointId+1);
    quadraticLinearWedge->GetPointIds()->SetId(2,pointId+2);
    quadraticLinearWedge->GetPointIds()->SetId(3,pointId+3);
    quadraticLinearWedge->GetPointIds()->SetId(4,pointId+4);
    quadraticLinearWedge->GetPointIds()->SetId(5,pointId+5);
    
    quadraticLinearWedge->GetPointIds()->SetId(6,pointId+11);
    quadraticLinearWedge->GetPointIds()->SetId(7,pointId+12);
    quadraticLinearWedge->GetPointIds()->SetId(8,pointId+13);
    quadraticLinearWedge->GetPointIds()->SetId(9,pointId+14);
    quadraticLinearWedge->GetPointIds()->SetId(10,pointId+15);
    quadraticLinearWedge->GetPointIds()->SetId(11,pointId+16);  
    
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticLinearWedge->GetCellType(),
                         quadraticLinearWedge->GetPointIds());
    quadraticLinearWedge->Delete();
  
    // this wedge shares a triangle face
    quadraticLinearWedge=vtkQuadraticLinearWedge::New();
    quadraticLinearWedge->GetPointIds()->SetId(0,pointId+3);
    quadraticLinearWedge->GetPointIds()->SetId(1,pointId+4);
    quadraticLinearWedge->GetPointIds()->SetId(2,pointId+5);
    quadraticLinearWedge->GetPointIds()->SetId(3,pointId+6);
    quadraticLinearWedge->GetPointIds()->SetId(4,pointId+7);
    quadraticLinearWedge->GetPointIds()->SetId(5,pointId+8);
    
    quadraticLinearWedge->GetPointIds()->SetId(6,pointId+14);
    quadraticLinearWedge->GetPointIds()->SetId(7,pointId+15);
    quadraticLinearWedge->GetPointIds()->SetId(8,pointId+16);
    quadraticLinearWedge->GetPointIds()->SetId(9,pointId+17);
    quadraticLinearWedge->GetPointIds()->SetId(10,pointId+18);
    quadraticLinearWedge->GetPointIds()->SetId(11,pointId+19); 
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticLinearWedge->GetCellType(),
                         quadraticLinearWedge->GetPointIds());
    quadraticLinearWedge->Delete();
  
    // this wedge shares a quad face
    quadraticLinearWedge=vtkQuadraticLinearWedge::New();
    quadraticLinearWedge->GetPointIds()->SetId(0,pointId+2);
    quadraticLinearWedge->GetPointIds()->SetId(1,pointId+1);
    quadraticLinearWedge->GetPointIds()->SetId(2,pointId+9);
    quadraticLinearWedge->GetPointIds()->SetId(3,pointId+5);
    quadraticLinearWedge->GetPointIds()->SetId(4,pointId+4);
    quadraticLinearWedge->GetPointIds()->SetId(5,pointId+10);
    
    quadraticLinearWedge->GetPointIds()->SetId(6,pointId+12);
    quadraticLinearWedge->GetPointIds()->SetId(7,pointId+20);
    quadraticLinearWedge->GetPointIds()->SetId(8,pointId+21);
    quadraticLinearWedge->GetPointIds()->SetId(9,pointId+15);
    quadraticLinearWedge->GetPointIds()->SetId(10,pointId+22);
    quadraticLinearWedge->GetPointIds()->SetId(11,pointId+23); 
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(quadraticLinearWedge->GetCellType(),
                         quadraticLinearWedge->GetPointIds());
    quadraticLinearWedge->Delete();
  
    pointId+=24;
    
    // 3D: bi quadratic quadratic wedge: 3 wedges, some share a biquadratic
    // quad face, some share a quadratic triangle face
    xOffset+=2.0;
    
    // corner points
    // triangle face of the first wedge
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 2
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 4
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 5
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge
    // 6
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 8
    points->InsertNextPoint(xOffset+0.5,yOffset+0.9,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // remaining vertices of the third wedge
    // 9
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points
    
    // triangle face of the first wedge: id=11
    // 11
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 12
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 13
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // common triangle face: id=14
    // 14
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 15
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 16
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the first wedge: id=17
    // 17
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 18
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 19
    points->InsertNextPoint(xOffset+0.5,yOffset+1.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // triangle face of the second wedge: id=20
    // 20
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 21
    points->InsertNextPoint(xOffset+0.7,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 22
    points->InsertNextPoint(xOffset+0.3,yOffset+0.5,-1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the second wedge: id=23
    // 23
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 24
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 25
    points->InsertNextPoint(xOffset+0.5,yOffset+0.95,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // first triangle face of the third wedge: id=26
    // 1+9
    // 26
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 27
    // 2+9
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,1.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    // second triangle face of the third wedge: id=28
    // 28
    // 4+10
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 29
    // 5+10
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other edges of the third wedge: id=30
    // 30
    // 9+10
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // face-centered points
    // 31=(0+1+3+4)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 32=(1+2+4+5)/4
    points->InsertNextPoint(xOffset+0.75,yOffset+0.5,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 33=(0+2+3+5)/4
    points->InsertNextPoint(xOffset+0.25,yOffset+0.5,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 34=(3+4+6+7)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.05,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 35=(4+5+7+8)/4
    points->InsertNextPoint(xOffset+0.725,yOffset+0.5,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    // 36=(3+5+6+8)/4
    points->InsertNextPoint(xOffset+0.275,yOffset+0.5,-0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 37=(1+4+9+10)/4
    points->InsertNextPoint(xOffset+1.0,yOffset+0.5,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 38=(2+5+9+10)/4
    points->InsertNextPoint(xOffset+0.75,yOffset+1.0,0.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    vtkBiQuadraticQuadraticWedge *biQuadraticQuadraticWedge;
    biQuadraticQuadraticWedge=vtkBiQuadraticQuadraticWedge::New();
    biQuadraticQuadraticWedge->GetPointIds()->SetId(0,pointId);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(1,pointId+1);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(2,pointId+2);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(3,pointId+3);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(4,pointId+4);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(5,pointId+5);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(6,pointId+11);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(7,pointId+12);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(8,pointId+13);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(9,pointId+14);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(10,pointId+15);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(11,pointId+16);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(12,pointId+17);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(13,pointId+18);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(14,pointId+19);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(15,pointId+31);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(16,pointId+32);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(17,pointId+33);
   
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticWedge->GetCellType(),
                         biQuadraticQuadraticWedge->GetPointIds());
    biQuadraticQuadraticWedge->Delete();
  
    // this wedge shares a triangle face
    biQuadraticQuadraticWedge=vtkBiQuadraticQuadraticWedge::New();
    biQuadraticQuadraticWedge->GetPointIds()->SetId(0,pointId+3);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(1,pointId+4);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(2,pointId+5);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(3,pointId+6);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(4,pointId+7);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(5,pointId+8);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(6,pointId+14);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(7,pointId+15);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(8,pointId+16);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(9,pointId+20);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(10,pointId+21);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(11,pointId+22);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(12,pointId+23);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(13,pointId+24);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(14,pointId+25);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(15,pointId+34);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(16,pointId+35);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(17,pointId+36);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticWedge->GetCellType(),
                         biQuadraticQuadraticWedge->GetPointIds());
    biQuadraticQuadraticWedge->Delete();
  
    // this wedge shares a quad face
    biQuadraticQuadraticWedge=vtkBiQuadraticQuadraticWedge::New();
    biQuadraticQuadraticWedge->GetPointIds()->SetId(0,pointId+2);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(1,pointId+1);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(2,pointId+9);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(3,pointId+5);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(4,pointId+4);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(5,pointId+10);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(6,pointId+12);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(7,pointId+26);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(8,pointId+27);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(9,pointId+15);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(10,pointId+28);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(11,pointId+29);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(12,pointId+19);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(13,pointId+18);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(14,pointId+30);
    
    biQuadraticQuadraticWedge->GetPointIds()->SetId(15,pointId+32);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(16,pointId+37);
    biQuadraticQuadraticWedge->GetPointIds()->SetId(17,pointId+38);
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticWedge->GetCellType(),
                         biQuadraticQuadraticWedge->GetPointIds());
    biQuadraticQuadraticWedge->Delete();
  
    pointId+=39;
    
    // 3D: biquadraticquadratic hexahedron: 3 with a some with a common face
    // with no center point, and other with a face witha center point.
    xOffset+=2.0;
    
    // a face (back): 0-3
    // 0
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 1
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 2
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 3
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (common): 4-7
    // 4
    points->InsertNextPoint(xOffset+0.1,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 5
    points->InsertNextPoint(xOffset+0.9,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 6
    points->InsertNextPoint(xOffset+0.9,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 7
    points->InsertNextPoint(xOffset+0.1,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // a face (front): 8-11
    // 8
    points->InsertNextPoint(xOffset+0.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 9
    points->InsertNextPoint(xOffset+1.0,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 10
    points->InsertNextPoint(xOffset+1.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 11
    points->InsertNextPoint(xOffset+0.0,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the back face: 12-15
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points on the common face: 16-19
    points->InsertNextPoint(xOffset+0.5,yOffset+0.1,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.9,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+1.9,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.1,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // mid-points between the back and common face: 20-23
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    
    // mid-points on the front face: 24-27
    points->InsertNextPoint(xOffset+0.5,yOffset+0.0,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+1.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.5,yOffset+2.0,5.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.0,yOffset+1.0,4.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // mid-points between the common face and the front face: 28-31
    points->InsertNextPoint(xOffset+0.05,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.95,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    points->InsertNextPoint(xOffset+0.05,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // face-centered points
    // 32=(0+1+4+5)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.05,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    // 33=(1+2+5+6)/4
    points->InsertNextPoint(xOffset+0.95,yOffset+1.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 34=(2+3+6+7)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.95,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 35=(0+3+4+7)/4
    points->InsertNextPoint(xOffset+0.05,yOffset+1.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // other hexa
    // 36=(4+5+8+9)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+0.05,3.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
  
    // 37=(5+6+9+10)/4
    points->InsertNextPoint(xOffset+0.95,yOffset+1.0,3.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 38=(6+7+10+11)/4
    points->InsertNextPoint(xOffset+0.5,yOffset+1.95,4.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 39=(4+7+8+11)/4
    points->InsertNextPoint(xOffset+0.05,yOffset+1.0,3.75);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // Third hexahedron
    // 40
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // 41
    points->InsertNextPoint(xOffset+2.0,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 42
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 43
    points->InsertNextPoint(xOffset+2.0,yOffset+2.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // face-edge: 44,45,46
    // 44=(1+40)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+0.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // 45=(40+41)/2
    points->InsertNextPoint(xOffset+2.0,yOffset+1.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 46=(1+40)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+2.0,0.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // face-edge: 47,48,49
    // 47=(5+42)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+0.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    // 48=(42+43)/2
    points->InsertNextPoint(xOffset+2.0,yOffset+1.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 49=(6+43)/2
    points->InsertNextPoint(xOffset+1.5,yOffset+2.0,3.0);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // cross-face-edge: 50,51
    // 50=(42+40)/2
    points->InsertNextPoint(xOffset+2.0,yOffset+0.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    // 51=(41+43)/2
    points->InsertNextPoint(xOffset+2.0,yOffset+2.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    
    
    // center: 52,53,54
    // 52=(40+41+42+43)/4
    points->InsertNextPoint(xOffset+2.0,yOffset+1.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 53=(1+5+40+42)/4
    points->InsertNextPoint(xOffset+1.5,yOffset+0.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
    // 54=(2+6+41+43)/4
    points->InsertNextPoint(xOffset+1.5,yOffset+2.0,1.5);
    scalars->InsertNextValue(scalar);
    scalar+=scalarStep;
   
    
    
    
    vtkBiQuadraticQuadraticHexahedron *biQuadraticQuadraticHexahedron;
    biQuadraticQuadraticHexahedron=vtkBiQuadraticQuadraticHexahedron::New();
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(0,pointId);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(1,pointId+1);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(2,pointId+2);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(3,pointId+3);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(4,pointId+4);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(5,pointId+5);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(6,pointId+6);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(7,pointId+7);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(8,pointId+12);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(9,pointId+13);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(10,pointId+14);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(11,pointId+15);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(12,pointId+16);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(13,pointId+17);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(14,pointId+18);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(15,pointId+19);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(16,pointId+20);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(17,pointId+21);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(18,pointId+22);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(19,pointId+23);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(20,pointId+35);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(21,pointId+33);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(22,pointId+32);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(23,pointId+34);
   
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticHexahedron->GetCellType(),
                         biQuadraticQuadraticHexahedron->GetPointIds());
    biQuadraticQuadraticHexahedron->Delete();
    
    biQuadraticQuadraticHexahedron=vtkBiQuadraticQuadraticHexahedron::New();
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(0,pointId+4);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(1,pointId+5);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(2,pointId+6);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(3,pointId+7);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(4,pointId+8);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(5,pointId+9);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(6,pointId+10);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(7,pointId+11);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(8,pointId+16);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(9,pointId+17);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(10,pointId+18);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(11,pointId+19);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(12,pointId+24);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(13,pointId+25);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(14,pointId+26);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(15,pointId+27);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(16,pointId+28);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(17,pointId+29);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(18,pointId+30);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(19,pointId+31);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(20,pointId+39);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(21,pointId+37);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(22,pointId+36);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(23,pointId+38);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticHexahedron->GetCellType(),
                         biQuadraticQuadraticHexahedron->GetPointIds());
    biQuadraticQuadraticHexahedron->Delete();
  
    biQuadraticQuadraticHexahedron=vtkBiQuadraticQuadraticHexahedron::New();
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(0,pointId+1);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(1,pointId+40);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(2,pointId+41);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(3,pointId+2);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(4,pointId+5);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(5,pointId+42);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(6,pointId+43);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(7,pointId+6);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(8,pointId+44);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(9,pointId+45);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(10,pointId+46);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(11,pointId+13);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(12,pointId+47);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(13,pointId+48);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(14,pointId+49);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(15,pointId+17);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(16,pointId+21);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(17,pointId+50);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(18,pointId+51);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(19,pointId+22);
    
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(20,pointId+33);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(21,pointId+52);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(22,pointId+53);
    biQuadraticQuadraticHexahedron->GetPointIds()->SetId(23,pointId+54);
    
    cellIds->InsertNextValue(cellId);
    ++cellId;
    grid->InsertNextCell(biQuadraticQuadraticHexahedron->GetCellType(),
                         biQuadraticQuadraticHexahedron->GetPointIds());
    biQuadraticQuadraticHexahedron->Delete();
  
    pointId+=55;
    
    
    grid->SetPoints(points);
    grid->GetPointData()->SetScalars(scalars);
    grid->GetCellData()->SetScalars(cellIds);
    points->Delete();
    scalars->Delete();
    cellIds->Delete();
  #endif
    
    // Create the filter
    vtkUnstructuredGridGeometryFilter *geom;
    geom=vtkUnstructuredGridGeometryFilter::New();
  #ifdef READ_FILE
    geom->SetInputConnection(0,reader->GetOutputPort(0));
  #else
    geom->SetInput(grid);
  #endif
    geom->Update(); //So that we can call GetRange() on the scalars
    
    // Check that all the cells are not 3D.
    // TODO
  #ifdef USE_SHRINK
    vtkShrinkFilter *shrink=vtkShrinkFilter::New();
    shrink->SetShrinkFactor(0.5);
    shrink->SetInputConnection(0,geom->GetOutputPort(0));
  #endif
  #ifdef FAST_GEOMETRY
    vtkDataSetSurfaceFilter *linearGeom=vtkDataSetSurfaceFilter::New();
  #else
    vtkGeometryFilter *linearGeom=vtkGeometryFilter::New();
  #endif
  #ifdef USE_SHRINK
    linearGeom->SetInputConnection(0,shrink->GetOutputPort(0));
  #else
    linearGeom->SetInputConnection(0,geom->GetOutputPort(0));
  #endif
    linearGeom->Update(); //So that we can call GetRange() on the scalars
   
    // This creates a blue to red lut.
    vtkLookupTable *lut = vtkLookupTable::New(); 
    lut->SetHueRange (0.667, 0.0);
  
    vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
    mapper->SetLookupTable(lut);
    mapper->SetInputConnection(0, linearGeom->GetOutputPort(0) );
    
   if(linearGeom->GetOutput()->GetPointData()!=0)
      {
      if(linearGeom->GetOutput()->GetPointData()->GetScalars()!=0)
        {
        mapper->SetScalarRange( linearGeom->GetOutput()->GetPointData()->
                                GetScalars()->GetRange());
        }
      }
  
    vtkActor *actor = vtkActor::New();
  //  cout<<"prop="<<actor->GetProperty()->GetBackfaceCulling()<<endl;
  #ifdef USE_CULLING
    actor->GetProperty()->SetBackfaceCulling(1);
  #endif
    actor->SetMapper(mapper);
    renderer->AddActor(actor);
    
  #ifdef WRITE_RESULT
    // Save the result of the filter in a file
    vtkXMLUnstructuredGridWriter *writer3d=vtkXMLUnstructuredGridWriter::New();
    writer3d->SetInputConnection(0,geom->GetOutputPort(0));
    writer3d->SetFileName("surface3d.vtu");
    writer3d->SetDataModeToAscii();
    writer3d->Write();
    writer3d->Delete();
  #endif // #ifdef WRITE_RESULT
    
    // Standard testing code.
    renderer->SetBackground(0.5,0.5,0.5);
    renWin->SetSize(300,300);
  //  vtkCamera *cam=renderer->GetActiveCamera();
  //  renderer->ResetCamera();
  //  cam->Azimuth(180);
    
    renWin->Render();
    int retVal = vtkRegressionTestImage( renWin );
    if ( retVal == vtkRegressionTester::DO_INTERACTOR)
      {
      iren->Start();
      }
  
    // Cleanup
    renderer->Delete();
    renWin->Delete();
    iren->Delete();
    
    mapper->Delete();
    actor->Delete();
    linearGeom->Delete();
    geom->Delete();
    lut->Delete();
  #ifdef USE_SHRINK
    shrink->Delete();
  #endif
  #ifdef READ_FILE
    reader->Delete();
  #else
    grid->Delete();
  #endif
    
    return !retVal;
  }
  