#include "GR_config.h"
#include "GR_assert.h"
#include "GR_CellCV.h"
#ifdef BROKEN_INLINE

CellCV::CellCV(const int iNVerts, Vert** const ppVIn)
{
  assert(iNVerts >= 2 && iNVerts <= 8);
  iNV = iNVerts;
  ppVVerts = new Vert*[iNVerts];
  int i;
  if (ppVIn)
    for (i = 0; i < iNVerts; i++) {
      ppVVerts[i] = ppVIn[i];
      assert2(ppVVerts[i] != pVInvalidVert, "Uninitialized vert");
    }
  else
    for (i = 0; i < iNVerts; i++)
      ppVVerts[i] = pVInvalidVert;
}

int CellCV::iFullCheck() const
{
  int iRetVal = qValid();
  int iNVerts = iNumVerts();
  for (int i = 0; iRetVal && i < iNVerts; i++) {
    Vert* pVThisVert = pVVert(i);
    iRetVal = iRetVal && pVThisVert->qValid();
  }
  return (iRetVal);
}

bool CellCV::qHasVert(const Vert* pV) const
{
  assert(pV->qValid());
  for (int i = 0; i < iNV; i++)
    if (ppVVerts[i] == pV) return(true);
  return(false);
}

void CellCV::vAssign(const int iNVerts, Vert** const ppVIn)
{
  assert(iNVerts >= 2 && iNVerts <= 8);
  if (iNV != iNVerts) {
    delete [] ppVVerts;
    iNV = iNVerts;
    ppVVerts = new Vert*[iNVerts];
  }
  int i;
  for (i = 0; i < iNVerts; i++) {
    ppVVerts[i] = ppVIn[i];
    assert2(ppVVerts[i] != pVInvalidVert, "Uninitialized vert");
  }
}

void CellCV::vReplaceVert(Vert* const pVOld, Vert* const pVNew)
{
  int i, iRep = -1;
  for (i = 0; i < iNV; i++) {
    if (ppVVerts[i] == pVOld) iRep = i;
    assert(ppVVerts[i] != pVNew);
  }
  assert(iRep != -1);
  ppVVerts[iRep] = pVNew;
}

#endif

double CellCV::dSize() const
{
  // Different stuff for different types of CellCV!
  int iNF = iNumFaces();
  switch (iNV) {
  case 3:
    // Triangle.
    {
      assert(iNF == 3);
      const Vert *pV0 = pVVert(0);
      const Vert *pV1 = pVVert(1);
      const Vert *pV2 = pVVert(2);
      int iDim = pV0->iSpaceDimen();
      double adEdge01[] = {0,0,0};
      double adEdge12[] = {0,0,0};
      for (int i = 0; i < iDim; i++) {
	adEdge01[i] = pV1->adCoords()[i] - pV0->adCoords()[i];
	adEdge12[i] = pV2->adCoords()[i] - pV1->adCoords()[i];
      }
      if (iDim == 2) {
	return (0.5 * fabs(dCROSS2D(adEdge01, adEdge12)));
      }
      else {
	double adRes[3];
	vCROSS3D(adEdge01, adEdge12, adRes);
	return 0.5 * dMAG3D(adRes);
      } // Triangle in 2D or 3D
    }
  case 4:
    assert(iNF == 4);
    // Tet or quad.  No way to tell, for quads in a surface mesh!
    // Assume it's a tet if it's 3D.
    if (pVVert(0)->iSpaceDimen() == 2) {
      assert(0);
      return -1;
    }
    else {
      const Vert *pV0 = pVVert(0);
      const Vert *pV1 = pVVert(1);
      const Vert *pV2 = pVVert(2);
      const Vert *pV3 = pVVert(3);

      double adVec1[] = { pV1->dX() - pV0->dX(),
			  pV1->dY() - pV0->dY(),
			  pV1->dZ() - pV0->dZ() };
      double adVec2[] = { pV2->dX() - pV0->dX(),
			  pV2->dY() - pV0->dY(),
			  pV2->dZ() - pV0->dZ() };
      double adVec3[] = { pV3->dX() - pV0->dX(),
			  pV3->dY() - pV0->dY(),
			  pV3->dZ() - pV0->dZ() };
      double adCross[3];
      vCROSS3D(adVec2, adVec3, adCross);
      return dDOT3D(adCross, adVec1) / 6.;
    }
  case 5:
    // Pyramid
    assert(iNF == 5);
    assert(0);
    return -1;
  case 6:
    // Prism
    assert(iNF == 5);
    assert(0);
    return -1;
  case 8:
    // Hex
    assert(iNF == 6);
    assert(0);
    return -1;
  default:
    assert(0);
    return -1;
  }
  // Should never reach here, but some compilers whine anyway.
//   assert(0);
//   return -1;
}

void CellCV::vAllDihed(double adDihed[], int* const piNDihed) const
{
  // Different stuff for different types of CellCV!
  int iNF = iNumFaces();
  if (iNV == 3) {
    assert(iNF == 3);
    // Triangle.
    const Vert *pV0 = pVVert(0);
    const Vert *pV1 = pVVert(1);
    const Vert *pV2 = pVVert(2);
    int iDim = pV0->iSpaceDimen();
    double adEdge01[] = {0,0,0};
    double adEdge12[] = {0,0,0};
    double adEdge20[] = {0,0,0};
    for (int i = 0; i < iDim; i++) {
      adEdge01[i] = pV1->adCoords()[i] - pV0->adCoords()[i];
      adEdge12[i] = pV2->adCoords()[i] - pV1->adCoords()[i];
      adEdge20[i] = pV0->adCoords()[i] - pV2->adCoords()[i];
    }
    vNORMALIZE3D(adEdge01);
    vNORMALIZE3D(adEdge12);
    vNORMALIZE3D(adEdge20);
    double dArg0 = -dDOT3D(adEdge20, adEdge01);
    double dArg1 = -dDOT3D(adEdge01, adEdge12);
    double dArg2 = -dDOT3D(adEdge12, adEdge20);

    double dScale = 180. / M_PI;
    adDihed[0] = GR_acos (dArg0) * dScale; // angle @ 0
    adDihed[1] = GR_acos (dArg1) * dScale; // angle @ 1
    adDihed[2] = GR_acos (dArg2) * dScale; // angle @ 2
    *piNDihed = 3;
  } // Triangle in 2D or 3D
  else if (iNV == 4) {
    // Tet or quad.  No way to tell, for quads in a surface mesh!  So
    // assume tets in 3D. 
    if (pVVert(0)->iSpaceDimen() == 2) {
      // Quad in 2D
    }
    else {
      // Tet in 3D
      const TetCellCV *pTCCV = reinterpret_cast<const TetCellCV*>(this);
      //      (pVVert(0), pVVert(1), pVVert(2), pVVert(3));
      pTCCV->TetCellCV::vAllDihed(adDihed, piNDihed);
      assert(*piNDihed == 6);
    }
  }
  else if (iNV == 5) {
    // Pyramid
    assert(0);
  }
  else if (iNV == 6) {
    // Prism
    assert(0);
  }
  else if (iNV == 8) {
    // Hex
    assert(0);
  }
}

void CellCV::vAllSolid(double adDihed[], int* const piNDihed) const
{
  // Different stuff for different types of CellCV, and we have to
  // deduce which type we have!
  int iNF = iNumFaces();
  if (iNV == 3) {
    assert(iNF == 3);
    // Triangle.
    vAllDihed(adDihed, piNDihed);
  } // Triangle in 2D or 3D
  else if (iNV == 4) {
    // Tet or quad.  No way to tell, for quads in a surface mesh!
    assert(0);
  }
  else if (iNV == 5) {
    // Pyramid
    assert(0);
  }
  else if (iNV == 6) {
    // Prism
    assert(0);
  }
  else if (iNV == 8) {
    // Hex
    assert(0);
  }
}
