text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    vtkBridgeCell.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.
  
  =========================================================================*/
  // .NAME vtkBridgeCell - Implementation of vtkGenericAdaptorCell
  // .SECTION Description
  // It is just an example that show how to implement the Generic. It is also
  // used for testing and evaluating the Generic.
  // .SECTION See Also
  // vtkGenericAdaptorCell, vtkBridgeDataSet
  
  
  #include "vtkBridgeCell.h"
  
  #include <cassert>
  
  #include "vtkBridgeCellIterator.h"
  #include "vtkObjectFactory.h"
  #include "vtkBridgeDataSet.h"
  #include "vtkDataSet.h"
  #include "vtkDoubleArray.h"
  #include "vtkGenericAttributeCollection.h"
  #include "vtkBridgeAttribute.h"
  #include "vtkContourValues.h"
  #include "vtkPoints.h"
  #include "vtkCellArray.h"
  #include "vtkDataSetAttributes.h"
  #include "vtkBridgePointIterator.h"
  
  // All that stuff is for InterpolationFunction()
  
  #include "vtkEmptyCell.h"
  #include "vtkVertex.h"
  #include "vtkPolyVertex.h"
  #include "vtkLine.h"
  #include "vtkPolyLine.h"
  #include "vtkTriangle.h"
  #include "vtkTriangleStrip.h"
  #include "vtkQuad.h"
  #include "vtkPixel.h"
  #include "vtkPolygon.h"
  #include "vtkTetra.h"
  #include "vtkHexahedron.h"
  #include "vtkVoxel.h"
  #include "vtkWedge.h"
  #include "vtkPyramid.h"
  
  #include "vtkQuadraticEdge.h"
  #include "vtkQuadraticTriangle.h"
  #include "vtkQuadraticQuad.h"
  #include "vtkQuadraticTetra.h"
  #include "vtkQuadraticHexahedron.h"
  #include "vtkBiQuadraticTriangle.h"
  #include "vtkConvexPointSet.h"
  # include "vtkPentagonalPrism.h"
  # include "vtkHexagonalPrism.h"
  # include "vtkQuadraticWedge.h"
  # include "vtkQuadraticPyramid.h"
  
  
  
  vtkStandardNewMacro(vtkBridgeCell);
  
  //-----------------------------------------------------------------------------
  void vtkBridgeCell::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Unique identification number of the cell over the whole
  // data set. This unique key may not be contiguous.
  vtkIdType vtkBridgeCell::GetId()
  {
    return this->Id;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Does `this' a cell of a dataset? (otherwise, it is a boundary cell)
  int vtkBridgeCell::IsInDataSet()
  {
    return this->BoolIsInDataSet;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Type of the current cell.
  // \post (result==VTK_HIGHER_ORDER_EDGE)||
  //       (result==VTK_HIGHER_ORDER_TRIANGLE)||
  //       (result==VTK_HIGHER_ORDER_TETRAHEDRON)
  int vtkBridgeCell::GetType()
  {
    int result=0;
    switch(this->Cell->GetCellType())
      {
      case VTK_TRIANGLE:
      case VTK_QUADRATIC_TRIANGLE:
      case VTK_BIQUADRATIC_TRIANGLE:
        result=VTK_HIGHER_ORDER_TRIANGLE;
        break;
      case VTK_QUAD:
      case VTK_QUADRATIC_QUAD:
        result=VTK_HIGHER_ORDER_QUAD;
        break;
      case VTK_TETRA:
      case VTK_QUADRATIC_TETRA:
        result=VTK_HIGHER_ORDER_TETRAHEDRON;
        break;
      case VTK_VOXEL:
      case VTK_HEXAHEDRON:
      case VTK_QUADRATIC_HEXAHEDRON:
        result=VTK_HIGHER_ORDER_HEXAHEDRON;
        break;
      case VTK_WEDGE:
      case VTK_QUADRATIC_WEDGE:
        result=VTK_HIGHER_ORDER_WEDGE;
        break;
      case VTK_PYRAMID:
      case VTK_QUADRATIC_PYRAMID:
        result=VTK_HIGHER_ORDER_PYRAMID;
        break;
      case VTK_PENTAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      case VTK_HEXAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      default:
        assert("check: impossible case" && 0);
        break;
      }
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Topological dimension of the current cell.
  // \post valid_result: result>=0 && result<=3
  int vtkBridgeCell::GetDimension()
  {
    int result=this->Cell->GetCellDimension();
    assert("post: valid_result" && (result>=0)&&(result<=3));
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Interpolation order of the geometry.
  // \post positive_result: result>=0
  int vtkBridgeCell::GetGeometryOrder()
  {
    int result;
    if(this->Cell->IsLinear())
      {
      result=1;
      }
    else
      {
      result=2; // GetOrder() is missing in vtkCell...
      }
    assert("post: positive_result" && result>=0);
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Does the cell have no higher-order interpolation for geometry?
  // \post definition: result==(GetGeometryOrder()==1)
  int vtkBridgeCell::IsGeometryLinear()
  {
    int result=this->Cell->IsLinear();
    assert("post: definition" && result==(GetGeometryOrder()==1));
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Interpolation order of attribute `a' on the cell (may differ by cell).
  // \pre a_exists: a!=0
  // \post positive_result: result>=0
  int vtkBridgeCell::GetAttributeOrder(vtkGenericAttribute *vtkNotUsed(a))
  {
    int result=this->GetGeometryOrder();
    assert("post: positive_result" && result>=0);
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Does the attribute `a' have no higher-order interpolation for the cell?
  // \pre a_exists: a!=0
  // \post definition: result==(GetAttributeOrder()==1)
  int vtkBridgeCell::IsAttributeLinear(vtkGenericAttribute *a)
  {
    (void)a; // The attribute order is the order of the geometry.
    int result=this->IsGeometryLinear();
    assert("post: definition" && result==(GetAttributeOrder(a)==1));
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Is the cell primary (i.e. not composite) ?
  int vtkBridgeCell::IsPrimary()
  {
    return this->Cell->IsPrimaryCell();
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Number of points that compose the cell.
  // \post positive_result: result>=0
  int vtkBridgeCell::GetNumberOfPoints()
  {
    int result=this->Cell->GetNumberOfPoints();
    assert("post: positive_result" && result>=0);
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Number of boundaries of dimension `dim' (or all dimensions less than
  // GetDimension() if -1) of the cell.
  // \pre valid_dim_range: (dim==-1) || ((dim>=0)&&(dim<GetDimension()))
  // \post positive_result: result>=0
  int vtkBridgeCell::GetNumberOfBoundaries(int dim)
  {
    assert("pre: valid_dim_range" && ((dim==-1) ||((dim>=0)&&(dim<GetDimension()))));
  
    int result=0;
    if( (dim==0) && (this->GetDimension()>1) )
      {
      result += this->Cell->GetNumberOfPoints();
      if(!this->Cell->IsLinear())
        { // Old cell API treats mid-edge nodes as vertices; subtract those out:
        result -= this->Cell->GetNumberOfEdges();
        }
      }
    if( ((dim==-1) && (this->GetDimension()>1)) || (dim==1) )
      {
      result=result+this->Cell->GetNumberOfEdges();
      }
    if( ((dim==-1) && (this->GetDimension()>2)) || (dim==2) )
      {
      result=result+this->Cell->GetNumberOfFaces();
      }
  
    assert("post: positive_result" && result>=0);
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Accumulated number of DOF nodes of the current cell. A DOF node is
  // a component of cell with a given topological dimension. e.g.: a triangle
  // has 7 DOF: 1 face, 3 edges, 3 vertices. An hexahedron has 27 DOF:
  // 1 region, 6 faces, 12 edges, 8 vertices.
  // \post valid_result: result==GetNumberOfBoundaries(-1)+1
  int vtkBridgeCell::GetNumberOfDOFNodes()
  {
    return this->GetNumberOfBoundaries(-1)+1;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return the points of cell into `it'.
  // \pre it_exists: it!=0
  void vtkBridgeCell::GetPointIterator(vtkGenericPointIterator *it)
  {
    assert("pre: it_exists" && it!=0);
    static_cast<vtkBridgePointIterator *>(it)->InitWithCell(this);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Create an empty cell iterator.
  // \post result_exists: result!=0
  vtkGenericCellIterator *vtkBridgeCell::NewCellIterator()
  {
    vtkGenericCellIterator *result=vtkBridgeCellIterator::New();
    assert("post: result_exists" && result!=0);
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return in `boundaries' the cells of dimension `dim' (or all dimensions
  // less than GetDimension() if -1) that are part of the boundary of the cell.
  // \pre valid_dim_range: (dim==-1) || ((dim>=0)&&(dim<GetDimension()))
  // \pre boundaries_exist: boundaries!=0
  void vtkBridgeCell::GetBoundaryIterator(vtkGenericCellIterator *boundaries,
                                          int dim)
  {
    assert("pre: valid_dim_range" && ((dim==-1) ||((dim>=0)&&(dim<GetDimension()))));
    assert("pre: boundaries_exist" && boundaries!=0);
    static_cast<vtkBridgeCellIterator *>(boundaries)->InitWithCellBoundaries(this,dim);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Number of cells (dimension>boundary->GetDimension()) of the dataset
  // that share the boundary `boundary' of `this'.
  // `this' IS NOT INCLUDED.
  // \pre boundary_exists: boundary!=0
  // \pre real_boundary: !boundary->IsInDataSet()
  // \pre cell_of_the_dataset: IsInDataSet()
  // \pre boundary: HasBoundary(boundary)
  // \post positive_result: result>=0
  int vtkBridgeCell::CountNeighbors(vtkGenericAdaptorCell *boundary)
  {
    assert("pre: boundary_exists" && boundary!=0);
    assert("pre: real_boundary" && !boundary->IsInDataSet());
    assert("pre: cell_of_the_dataset" && IsInDataSet());
  
    vtkIdList *cells=vtkIdList::New();
    vtkBridgeCell *b=static_cast<vtkBridgeCell *>(boundary);
    vtkIdList *pts=b->Cell->GetPointIds();
    this->DataSet->Implementation->GetCellNeighbors(this->Id,pts,cells);
    int result=cells->GetNumberOfIds();
    cells->Delete();
  
    assert("post: positive_result" && result>=0);
  
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // \pre large_enough: GetDimension()>=2
  // \pre right_size: sizeof(sharing)==GetNumberOfBoundaries(1);
  void vtkBridgeCell::CountEdgeNeighbors(int *sharing)
  {
    assert("pre: large_enough" && this->GetDimension()>=2);
  
    vtkIdType c=this->Cell->GetNumberOfEdges();
    vtkIdList *cells=vtkIdList::New();
    vtkIdType i=0;
    vtkCell *edge;
    vtkIdList *pts;
  
    while(i<c)
      {
      edge=this->Cell->GetEdge(i); // edge is deleted automatically by this->Cell
      pts=edge->GetPointIds();
      this->DataSet->Implementation->GetCellNeighbors(this->Id,pts,cells);
      sharing[i]=cells->GetNumberOfIds();
      ++i;
      }
    cells->Delete();
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Put into `neighbors' the cells (dimension>boundary->GetDimension())
  // of the dataset that share the boundary `boundary' of `this'.
  // `this' IS NOT INCLUDED.
  // \pre boundary_exists: boundary!=0
  // \pre real_boundary: !boundary->IsInDataSet()
  // \pre cell_of_the_dataset: IsInDataSet()
  // \pre boundary: HasBoundary(boundary)
  // \pre neighbors_exist: neighbors!=0
  void vtkBridgeCell::GetNeighbors(vtkGenericAdaptorCell *boundary,
                                   vtkGenericCellIterator *neighbors)
  {
    assert("pre: boundary_exists" && boundary!=0);
    assert("pre: real_boundary" && !boundary->IsInDataSet());
    assert("pre: cell_of_the_dataset" && IsInDataSet());
    assert("pre: neighbors_exist" && neighbors!=0);
  
    vtkIdList *cells=vtkIdList::New();
    vtkIdList *pts=static_cast<vtkBridgeCell *>(boundary)->Cell->GetPointIds();
    this->DataSet->Implementation->GetCellNeighbors(this->Id,pts,cells);
  
    static_cast<vtkBridgeCellIterator *>(neighbors)->InitWithCells(cells,this->DataSet);
  
    cells->Delete();
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Compute the closest boundary of the current sub-cell `subId' for point
  // `pcoord' (in parametric coordinates) in `boundary', and return whether
  // the point is inside the cell or not. `boundary' is of dimension
  // GetDimension()-1.
  // \pre positive_subId: subId>=0
  int vtkBridgeCell::FindClosestBoundary(int subId,
                                         double pcoords[3],
                                         vtkGenericCellIterator* &boundary)
  {
    assert("pre: positive_subId" && subId>=0);
  
    vtkIdList *pts=vtkIdList::New();
    int result=this->Cell->CellBoundary(subId,pcoords,pts);
    static_cast<vtkBridgeCellIterator *>(boundary)->InitWithPoints(this->Cell->Points,pts,this->GetDimension()-1,0); // id of the boundary always 0?
    pts->Delete();
  
    return result;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Is `x' inside the current cell? It also evaluate parametric coordinates
  // `pcoords', sub-cell id `subId' (0 means primary cell), distance squared
  // to the sub-cell in `dist2' and closest corner point `closestPoint'.
  // `dist2' and `closestPoint' are not evaluated if `closestPoint'==0.
  // If a numerical error occurred, -1 is returned and all other results
  // should be ignored.
  // \post valid_result: result==-1 || result==0 || result==1
  // \post positive_distance: result!=-1 implies (closestPoint!=0 implies
  //                                               dist2>=0)
  int vtkBridgeCell::EvaluatePosition(double x[3],
                                      double *closestPoint,
                                      int &subId,
                                      double pcoords[3],
                                      double &dist2)
  {
    this->AllocateWeights();
    int result=this->Cell->EvaluatePosition(x,closestPoint,subId,pcoords,dist2,
                                            this->Weights);
  
    if(result)
      {
      // clamp pcoords
      int i=0;
      while(i<3)
        {
        if(pcoords[i]<0)
          {
          pcoords[i]=0;
          }
        else if(pcoords[i]>1)
          {
          pcoords[i]=1;
          }
        ++i;
        }
      }
  
    assert("post: valid_result" && (result==-1 || result==0 || result==1));
    assert("post: positive_distance" && (!(result!=-1) || (!(closestPoint!=0)||dist2>=0))); // A=>B: !A || B
    return result;
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Determine global coordinates `x' from sub-cell `subId' and parametric
  // coordinates `pcoords' in the cell.
  // \pre positive_subId: subId>=0
  // \pre clamped_pcoords: (0<=pcoords[0])&&(pcoords[0]<=1)&&(0<=pcoords[1])
  // &&(pcoords[1]<=1)&&(0<=pcoords[2])&&(pcoords[2]<=1)
  void vtkBridgeCell::EvaluateLocation(int subId,
                                       double pcoords[3],
                                       double x[3])
  {
    assert("pre: positive_subId" && subId>=0);
    assert("pre: clamped_pcoords" && (0<=pcoords[0])&&(pcoords[0]<=1)
               &&(0<=pcoords[1])&&(pcoords[1]<=1)&&(0<=pcoords[2])
               &&(pcoords[2]<=1));
  
    this->AllocateWeights();
    this->Cell->EvaluateLocation(subId,pcoords,x,this->Weights);
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Interpolate the attribute `a' at local position `pcoords' of the cell into
  // `val'.
  // \pre a_exists: a!=0
  // \pre a_is_point_centered: a->GetCentering()==vtkPointCentered
  // \pre clamped_point: pcoords[0]>=0 && pcoords[0]<=1 && pcoords[1]>=0 &&
  //                     pcoords[1]<=1 && pcoords[2]>=0 && pcoords[2]<=1
  // \pre val_exists: val!=0
  // \pre valid_size: sizeof(val)==a->GetNumberOfComponents()
  void vtkBridgeCell::InterpolateTuple(vtkGenericAttribute *a, double pcoords[3],
                                       double *val)
  {
    assert("pre: a_exists" && a!=0);
    assert("pre: a_is_point_centered" && a->GetCentering()==vtkPointCentered);
    assert("pre: clamped_point" && (pcoords[0]>=0 && pcoords[0]<=1
               && pcoords[1]>=0 && pcoords[1]<=1 && pcoords[2]>=0
                                    && pcoords[2]<=1));
    assert("pre: val_exists" && val!=0);
  
    vtkBridgeAttribute *ba = static_cast<vtkBridgeAttribute *>(a);
  
    int componentCount = a->GetNumberOfComponents();
    int ptCount = this->GetNumberOfPoints();
  
    if(a->GetCentering() == vtkPointCentered)
      {
      this->AllocateWeights();
      this->InterpolationFunctions(pcoords, this->Weights);
  
      memset(val,0, sizeof(double)*componentCount);
      for(int pt = 0; pt<ptCount; ++pt)
        {
        ba->Data->GetArray(ba->AttributeNumber)->
          GetTuple(this->Cell->GetPointId(pt),ba->InternalTuple);
        for(int component = 0; component<componentCount; ++component)
          {
          val[component] += ba->InternalTuple[component]*this->Weights[pt];
          }
        }
      }
    else // cell centered
      {
      // not need to interpolate
      ba->Data->GetArray(ba->AttributeNumber)->GetTuple(this->GetId(),val);
      }
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Interpolate the whole collection of attributes `c' at local position
  // `pcoords' of the cell into `val'. Only point centered attributes are
  // taken into account.
  // \pre c_exists: c!=0
  // \pre clamped_point: pcoords[0]>=0 && pcoords[0]<=1 && pcoords[1]>=0 &&
  //                     pcoords[1]<=1 && pcoords[2]>=0 && pcoords[2]<=1
  // \pre val_exists: val!=0
  // \pre valid_size: sizeof(val)==c->GetNumberOfPointCenteredComponents()
  void vtkBridgeCell::InterpolateTuple(vtkGenericAttributeCollection *c,
                                       double pcoords[3],
                                       double *val)
  {
    assert("pre: c_exists" && c!=0);
    assert("pre: clamped_point" && (pcoords[0]>=0 && pcoords[0]<=1
               && pcoords[1]>=0 && pcoords[1]<=1 && pcoords[2]>=0
                                    && pcoords[2]<=1));
    assert("pre: val_exists" && val!=0);
  
  ///  assert("check: used!" && 0);
  
    double *p=val;
    int i=0;
    int count=c->GetNumberOfAttributes();
    while(i<count)
      {
      if(c->GetAttribute(i)->GetCentering()==vtkPointCentered)
        {
        this->InterpolateTuple(c->GetAttribute(i),pcoords,p);
        p=p+c->GetAttribute(i)->GetNumberOfComponents();
        }
      ++i;
      }
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Is there an intersection between the current cell and the ray (`p1',`p2')
  // according to a tolerance `tol'? If true, `x' is the global intersection,
  // `t' is the parametric coordinate for the line, `pcoords' are the
  // parametric coordinates for cell. `subId' is the sub-cell where
  // the intersection occurs.
  // \pre positive_tolerance: tol>0
  int vtkBridgeCell::IntersectWithLine(double p1[3],
                                       double p2[3],
                                       double tol,
                                       double &t,
                                       double x[3],
                                       double pcoords[3],
                                       int &subId)
  {
    return this->Cell->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Compute derivatives `derivs' of the attribute `attribute' (from its
  // values at the corner points of the cell) given sub-cell `subId' (0 means
  // primary cell) and parametric coordinates `pcoords'.
  // Derivatives are in the x-y-z coordinate directions for each data value.
  // \pre positive_subId: subId>=0
  // \pre clamped_pcoords: (0<=pcoords[0])&&(pcoords[0]<=1)&&(0<=pcoords[1])
  // &&(pcoords[1]<=1)&&(0<=pcoords[2])%%(pcoords[2]<=1)
  // \pre attribute_exists: attribute!=0
  // \pre derivs_exists: derivs!=0
  // \pre valid_size: sizeof(derivs)>=attribute->GetNumberOfComponents()*3
  void vtkBridgeCell::Derivatives(int subId,
                                  double pcoords[3],
                                  vtkGenericAttribute *attribute,
                                  double *derivs)
  {
    double *tuples =
      new double[attribute->GetNumberOfComponents()*this->GetNumberOfPoints()];
    attribute->GetTuple(this->InternalIterator,tuples);
    this->Cell->Derivatives(subId,pcoords,tuples,
                            attribute->GetNumberOfComponents(),derivs);
    delete [] tuples;
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Compute the bounding box of the current cell in `bounds' in global
  // coordinates.
  // THREAD SAFE
  void vtkBridgeCell::GetBounds(double bounds[6])
  {
    this->Cell->GetBounds(bounds);
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Return the bounding box of the current cell in global coordinates.
  // NOT THREAD SAFE
  // \post result_exists: result!=0
  // \post valid_size: sizeof(result)>=6
  double *vtkBridgeCell::GetBounds()
  {
    return this->Cell->GetBounds();
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Bounding box diagonal squared of the current cell.
  // \post positive_result: result>=0
  double vtkBridgeCell::GetLength2()
  {
    return this->Cell->GetLength2();
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Center of the current cell in parametric coordinates `pcoords'.
  // If the current cell is a composite, the return value is the sub-cell id
  // that the center is in.
  // \post valid_result: (result>=0) && (IsPrimary() implies result==0)
  int vtkBridgeCell::GetParametricCenter(double pcoords[3])
  {
    return this->Cell->GetParametricCenter(pcoords);
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Distance of the parametric coordinate `pcoords' to the current cell.
  // If inside the cell, a distance of zero is returned. This is used during
  // picking to get the correct cell picked. (The tolerance will occasionally
  // allow cells to be picked who are not really intersected "inside" the
  // cell.)
  // \post positive_result: result>=0
  double vtkBridgeCell::GetParametricDistance(double pcoords[3])
  {
    return this->Cell->GetParametricDistance(pcoords);
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Return a contiguous array of parametric coordinates of the points defining
  // the current cell. In other words, (px,py,pz, px,py,pz, etc..) The
  // coordinates are ordered consistent with the definition of the point
  // ordering for the cell. Note that 3D parametric coordinates are returned
  // no matter what the topological dimension of the cell. It includes the DOF
  // nodes.
  // \post valid_result_exists: ((IsPrimary()) && (result!=0)) ||
  //                             ((!IsPrimary()) && (result==0))
  //                     result!=0 implies sizeof(result)==GetNumberOfPoints()
  double *vtkBridgeCell::GetParametricCoords()
  {
    return this->Cell->GetParametricCoords();
  }
  
  // For the internals of the tesselation algorithm (the hash table in particular)
  // Is the face `faceId' of the current cell on a exterior boundary of the
  // dataset or not?
  // \pre 3d: GetDimension()==3
  //----------------------------------------------------------------------------
  int vtkBridgeCell::IsFaceOnBoundary(vtkIdType faceId)
  {
    assert("pre: 3d" && this->GetDimension()==3);
  
    // result=CountNeighbors(boundary(faceId))==0;
  
    vtkCell *face = this->Cell->GetFace(faceId);
    vtkIdList *cells = vtkIdList::New(); // expensive
    this->DataSet->Implementation->
      GetCellNeighbors(this->Id,face->GetPointIds(),cells);
  
    int result=cells->GetNumberOfIds()==0;
    cells->Delete(); // expensive
  #if 0
    if(this->GetType()==VTK_QUADRATIC_TETRA)
      {
      if(result)
        {
        cout<<"************************************************ boundary"<<endl;
        }
      else
        {
        cout<<"************************************************ NOT boundary"<<endl;
        }
  //    assert(result);
      }
  #endif
    return result;
  }
  
  // Is the cell on the exterior boundary of the dataset?
  // \pre 2d: GetDimension()==2
  //----------------------------------------------------------------------------
  int vtkBridgeCell::IsOnBoundary()
  {
    assert("pre: 2d" && this->GetDimension()==2);
  //  assert("check: TODO" && 0);
    return 1;
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Put into `id' the list of ids the point of the cell.
  // \pre id_exists: id!=0
  // \pre valid_size: sizeof(id)==GetNumberOfPoints();
  void vtkBridgeCell::GetPointIds(vtkIdType *id)
  {
    vtkIdType i=0;
    vtkIdList *l=this->Cell->GetPointIds();
    vtkIdType c=this->GetNumberOfBoundaries(0);
    while(i<c)
      {
       id[i]=l->GetId(i);
      ++i;
      }
  }
  //----------------------------------------------------------------------------
  // Description:
  // Return the ids of the vertices defining face `faceId'.
  // \pre is_3d: this->GetDimension()==3
  // \pre valid_faceId_range: faceId>=0 && faceId<this->GetNumberOfBoundaries(2)
  // \post result_exists: result!=0
  // \post valid_size: sizeof(result)>=GetNumberOfVerticesOnFace(faceId)
  int *vtkBridgeCell::GetFaceArray(int faceId)
  {
    assert("pre: is_3d" && this->GetDimension()==3);
    assert("pre: valid_faceId_range" && faceId>=0
      && faceId<this->GetNumberOfBoundaries(2));
  
    int *result = 0;
  
    switch(this->GetType())
      {
      case VTK_HIGHER_ORDER_TETRAHEDRON:
        result = vtkTetra::GetFaceArray(faceId);
        break;
      case VTK_HIGHER_ORDER_HEXAHEDRON:
        if(this->Cell->GetCellType()==VTK_VOXEL)
          {
          result = vtkVoxel::GetFaceArray(faceId);
          }
        else
          {
          result = vtkHexahedron::GetFaceArray(faceId);
          }
        break;
      case VTK_HIGHER_ORDER_WEDGE:
        result = vtkWedge::GetFaceArray(faceId);
        break;
      case VTK_HIGHER_ORDER_PYRAMID:
        result = vtkPyramid::GetFaceArray(faceId);
        break;
      case VTK_PENTAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      case VTK_HEXAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      default:
        assert("check: impossible case" && 0);
        break;
      }
    return result;
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Return the number of vertices defining face `faceId'
  // \pre is_3d: this->GetDimension()==3
  // \pre valid_faceId_range: faceId>=0 && faceId<this->GetNumberOfBoundaries(2)
  // \post positive_result: && result>0
  int vtkBridgeCell::GetNumberOfVerticesOnFace(int faceId)
  {
    assert("pre: is_3d" && this->GetDimension()==3);
    assert("pre: valid_faceId_range" && faceId>=0
      && faceId<this->GetNumberOfBoundaries(2));
  
    int result = 0;
  
    switch(this->GetType())
      {
      case VTK_HIGHER_ORDER_TETRAHEDRON:
        result = 3;
        break;
      case VTK_HIGHER_ORDER_HEXAHEDRON:
        result = 4;
        break;
      case  VTK_HIGHER_ORDER_WEDGE:
        if(faceId <= 1) // triangle face
          {
          result = 3;
          }
        else // quad face
          {
          result = 4;
          }
        break;
      case VTK_HIGHER_ORDER_PYRAMID:
        if( faceId == 0)  // base
          {
          result = 4;
          }
        else // side
          {
          result = 3;
          }
        break;
  #if 0 // TODO
      case VTK_PENTAGONAL_PRISM:
        if(faceId<=1)
          {
          result=4;
          }
        else
          {
          result=3;
          }
        break;
      case VTK_HEXAGONAL_PRISM:
         if(faceId<=1)
          {
          result=6;
          }
        else
          {
          result=4;
          }
        break;
  #endif
      default:
        assert("check: impossible case" && 0);
        break;
      }
  
    assert("post: positive_result" && result>0);
    return result;
  }
  
  // copy/paste of vtkTriangle.cxx
  static int triangleEdges[3][2] = { {0,1}, {1,2}, {2,0} };
  static int quadEdges[4][2] = { {0,1}, {1,2}, {3,2}, {0,3} };
  
  //----------------------------------------------------------------------------
  // Description:
  // Return the ids of the vertices defining edge `edgeId'.
  // \pre valid_dimension: this->GetDimension()>=2
  // \pre valid_edgeId_range: edgeId>=0 && edgeId<this->GetNumberOfBoundaries(1)
  // \post result_exists: result!=0
  // \post valid_size: sizeof(result)==2
  int *vtkBridgeCell::GetEdgeArray(int edgeId)
  {
    assert("pre: valid_dimension" && this->GetDimension()>=2);
    assert("pre: valid_faceId_range" && edgeId>=0
      && edgeId<this->GetNumberOfBoundaries(1));
  
    int *result = 0;
  
    switch(this->GetType())
      {
      case VTK_HIGHER_ORDER_TRIANGLE:
        result = triangleEdges[edgeId];
        break;
      case VTK_HIGHER_ORDER_QUAD:
        result = quadEdges[edgeId];
        break;
      case VTK_HIGHER_ORDER_TETRAHEDRON:
        result = vtkTetra::GetEdgeArray(edgeId);
        break;
      case VTK_HIGHER_ORDER_HEXAHEDRON:
        if(this->Cell->GetCellType()==VTK_VOXEL)
          {
          result = vtkVoxel::GetEdgeArray(edgeId);
          }
        else
          {
          result = vtkHexahedron::GetEdgeArray(edgeId);
          }
        break;
      case VTK_HIGHER_ORDER_WEDGE:
        result = vtkWedge::GetEdgeArray(edgeId);
        break;
      case VTK_HIGHER_ORDER_PYRAMID:
        result = vtkPyramid::GetEdgeArray(edgeId);
        break;
      case VTK_PENTAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      case VTK_HEXAGONAL_PRISM:
        assert("check: TODO" && 0);
        break;
      default:
        assert("check: impossible case" && 0);
        break;
      }
  
    return result;
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Used internally for the Bridge.
  // Initialize the cell from a dataset `ds' and `cellid'.
  // \pre ds_exists: ds!=0
  // \pre valid_cellid: (cellid>=0) && (cellid<ds->GetNumberOfCells())
  void vtkBridgeCell::Init(vtkBridgeDataSet *ds,
                           vtkIdType cellid)
  {
    assert("pre: ds_exists" && ds!=0);
    assert("pre: valid_cellid" && (cellid>=0)
           && (cellid<ds->GetNumberOfCells()));
  
    vtkSetObjectBodyMacro(DataSet,vtkBridgeDataSet,ds);
    vtkCell *tmp = ds->Implementation->GetCell(cellid);
    vtkSetObjectBodyMacro(Cell,vtkCell,tmp);
    this->Id = cellid;
    this->BoolIsInDataSet = 1;
    if(this->InternalIterator == 0)
      {
      this->InternalIterator = vtkBridgeCellIterator::New();
      }
    this->InternalIterator->InitWithOneCell(this);
  
    this->InternalIterator->Begin();
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Used internally for the Bridge.
  // Initialize the cell from a cell `c' and an `id'.
  // \pre c_exists: c!=0
  void vtkBridgeCell::InitWithCell(vtkCell *c, vtkIdType id)
  {
    assert("pre: c_exists" && c!=0);
  
    vtkSetObjectBodyMacro(DataSet,vtkBridgeDataSet,0);
    this->Id = id;
  
    // warning: do directly vtkSetObjectBodyMacro(Cell,vtkCell,c->NewInstance())
    // add memory leak because the text "c->NewInstance()" is copied several
    // time in the macro...
  
    if(this->Cell)
      {
      this->Cell->Delete();
      }
    this->Cell = c->NewInstance();
  
    this->Cell->DeepCopy(c);
    this->BoolIsInDataSet=0;
  
    if(this->InternalIterator==0)
      {
      this->InternalIterator=vtkBridgeCellIterator::New();
      }
    this->InternalIterator->InitWithOneCell(this);
    this->InternalIterator->Begin();
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Recursive copy of `other' into `this'.
  // \pre other_exists: other!=0
  // \pre other_differ: this!=other
  void vtkBridgeCell::DeepCopy(vtkBridgeCell *other)
  {
    assert("pre: other_exists" && other!=0);
    assert("pre: other_differ" && this!=other);
  
    vtkCell *tmp;
  
    if(this->InternalIterator==0)
      {
      this->InternalIterator=vtkBridgeCellIterator::New();
      }
    this->Id = other->Id;
    this->BoolIsInDataSet = other->BoolIsInDataSet;
    if(other->BoolIsInDataSet)
      {
      vtkSetObjectBodyMacro(DataSet,vtkBridgeDataSet,other->DataSet);
      tmp = this->DataSet->Implementation->GetCell(this->Id);
      vtkSetObjectBodyMacro(Cell,vtkCell,tmp);
      this->InternalIterator->InitWithOneCell(this);
      this->InternalIterator->Begin();
      }
    else
      {
      vtkSetObjectBodyMacro(DataSet,vtkBridgeDataSet,0);
      tmp = other->Cell->NewInstance();
      vtkSetObjectBodyMacro(Cell,vtkCell,tmp);
      this->Cell->Delete(); // because newinstance+macro=2 ref
      this->Cell->DeepCopy(other->Cell);
      this->InternalIterator->InitWithOneCell(this);
      this->InternalIterator->Begin();
      }
    this->Modified();
  }
  
  //----------------------------------------------------------------------------
  vtkBridgeCell::vtkBridgeCell()
  {
    this->DataSet = 0;
    this->InternalIterator = 0; // we cannot create the cell iterator here
    // because we will have an infinite recursion: a cell creates a
    // celliterator which creates a cell, which creates a celliterator ...
    this->Cell = 0;
    this->BoolIsInDataSet = 0;
    this->Id = -1000; // magic ?
  
    this->Weights = 0;
    this->WeightsCapacity = 0;
  }
  
  //----------------------------------------------------------------------------
  vtkBridgeCell::~vtkBridgeCell()
  {
    vtkSetObjectBodyMacro(DataSet,vtkBridgeDataSet,0);
    vtkSetObjectBodyMacro(InternalIterator,vtkBridgeCellIterator,0);
    vtkSetObjectBodyMacro(Cell,vtkCell,0);
  
    if(this->Weights)
      {
      delete[] this->Weights;
      }
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Allocate an array for the weights, only if it does not exist yet or if
  // the capacity is too small.
  void vtkBridgeCell::AllocateWeights()
  {
    if( this->Weights != 0
      && this->WeightsCapacity < this->GetNumberOfPoints() )
      {
      delete[] this->Weights;
      this->Weights = 0;
      }
    if(this->Weights == 0)
      {
      this->Weights = new double[this->GetNumberOfPoints()];
      this->WeightsCapacity = this->GetNumberOfPoints();
      }
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Compute the weights for parametric coordinates `pcoords'.
  void vtkBridgeCell::InterpolationFunctions(double pcoords[3], double *weights)
  {
    this->Cell->InterpolateFunctions(pcoords, weights);
  }
  