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

  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 "vtkBridgeAttribute.h"
#include "vtkBridgeCellIterator.h"
#include "vtkBridgeDataSet.h"
#include "vtkBridgePointIterator.h"
#include "vtkCellArray.h"
#include "vtkContourValues.h"
#include "vtkDataSet.h"
#include "vtkDataSetAttributes.h"
#include "vtkDoubleArray.h"
#include "vtkGenericAttributeCollection.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"

// All that stuff is for InterpolationFunction()

#include "vtkEmptyCell.h"
#include "vtkHexahedron.h"
#include "vtkLine.h"
#include "vtkPixel.h"
#include "vtkPolyLine.h"
#include "vtkPolyVertex.h"
#include "vtkPolygon.h"
#include "vtkPyramid.h"
#include "vtkQuad.h"
#include "vtkTetra.h"
#include "vtkTriangle.h"
#include "vtkTriangleStrip.h"
#include "vtkVertex.h"
#include "vtkVoxel.h"
#include "vtkWedge.h"

#include "vtkBiQuadraticTriangle.h"
#include "vtkConvexPointSet.h"
#include "vtkHexagonalPrism.h"
#include "vtkPentagonalPrism.h"
#include "vtkQuadraticEdge.h"
#include "vtkQuadraticHexahedron.h"
#include "vtkQuadraticPyramid.h"
#include "vtkQuadraticQuad.h"
#include "vtkQuadraticTetra.h"
#include "vtkQuadraticTriangle.h"
#include "vtkQuadraticWedge.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)
vtkTypeBool 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 != nullptr);
  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 != nullptr);
  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 != nullptr);
  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 != nullptr);
  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 != nullptr);
  assert("pre: real_boundary" && !boundary->IsInDataSet());
  assert("pre: cell_of_the_dataset" && IsInDataSet());
  assert("pre: neighbors_exist" && neighbors != nullptr);

  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(
  const 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 != nullptr) || 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 != nullptr);
  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 != nullptr);

  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 != nullptr);
  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 != nullptr);

  ///  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(const 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 tessellation 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)
const vtkIdType* vtkBridgeCell::GetFaceArray(vtkIdType faceId)
{
  assert("pre: is_3d" && this->GetDimension() == 3);
  assert("pre: valid_faceId_range" && faceId >= 0 && faceId < this->GetNumberOfBoundaries(2));

  const vtkIdType* result = nullptr;

  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 constexpr vtkIdType triangleEdges[3][2] = { { 0, 1 }, { 1, 2 }, { 2, 0 } };
static constexpr vtkIdType 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
const vtkIdType* vtkBridgeCell::GetEdgeArray(vtkIdType edgeId)
{
  assert("pre: valid_dimension" && this->GetDimension() >= 2);
  assert("pre: valid_faceId_range" && edgeId >= 0 && edgeId < this->GetNumberOfBoundaries(1));

  const vtkIdType* result = nullptr;

  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 != nullptr);
  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 == nullptr)
  {
    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 != nullptr);

  vtkSetObjectBodyMacro(DataSet, vtkBridgeDataSet, static_cast<vtkBridgeDataSet*>(nullptr));
  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 == nullptr)
  {
    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 != nullptr);
  assert("pre: other_differ" && this != other);

  vtkCell* tmp;

  if (this->InternalIterator == nullptr)
  {
    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, static_cast<vtkBridgeDataSet*>(nullptr));
    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 = nullptr;
  this->InternalIterator = nullptr; // 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 = nullptr;
  this->BoolIsInDataSet = 0;
  this->Id = -1000; // magic ?

  this->Weights = nullptr;
  this->WeightsCapacity = 0;
}

//------------------------------------------------------------------------------
vtkBridgeCell::~vtkBridgeCell()
{
  vtkSetObjectBodyMacro(DataSet, vtkBridgeDataSet, static_cast<vtkBridgeDataSet*>(nullptr));
  vtkSetObjectBodyMacro(
    InternalIterator, vtkBridgeCellIterator, static_cast<vtkBridgeCellIterator*>(nullptr));
  vtkSetObjectBodyMacro(Cell, vtkCell, static_cast<vtkCell*>(nullptr));

  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 != nullptr && this->WeightsCapacity < this->GetNumberOfPoints())
  {
    delete[] this->Weights;
    this->Weights = nullptr;
  }
  if (this->Weights == nullptr)
  {
    this->Weights = new double[this->GetNumberOfPoints()];
    this->WeightsCapacity = this->GetNumberOfPoints();
  }
}

//------------------------------------------------------------------------------
// Description:
// Compute the weights for parametric coordinates `pcoords'.
void vtkBridgeCell::InterpolationFunctions(const double pcoords[3], double* weights)
{
  this->Cell->InterpolateFunctions(pcoords, weights);
}