#include <stdio.h>
#include "GR_misc.h"
#include "GR_events.h"
#include "GR_Face.h"
#include "GR_Cell.h"
#include "GR_BFace.h"
#include "GR_Vertex.h"
#include "GR_Mesh2D.h"
void vWriteFEA(Mesh2D& OutMesh,
	       const char strBaseFileName[],
	       const int iOrder,
	       const char strExtraFileSuffix[])
{
  SUMAA_LOG_EVENT_BEGIN(OUTPUT);
  assert(OutMesh.qSimplicial());
  FILE *pFOutFile = fopen("/dev/zero", "r");
  char strFileName[1024];
  OutMesh.vPurge();
  sprintf(strFileName, "%s%s.fea", strBaseFileName, strExtraFileSuffix);

  /*newfile fea*/
  printf("Opening output file %s\n", strFileName);
  if (pFOutFile) fclose(pFOutFile);
  pFOutFile = fopen(strFileName, "w");
  if (NULL == pFOutFile)
    vFatalError("Couldn't open output file for writing",
                "2d mesh output");

  /*nverts ncells nfaces nbfaces*/
  fprintf(pFOutFile, "%d %d %d %d\n", OutMesh.iNumVerts(),
	  OutMesh.iNumCells(), OutMesh.iNumFaces(), OutMesh.iNumBdryFaces());

  for (int iV = 0; iV < OutMesh.iNumVerts(); iV++) 
    /*verts: coords*/
    fprintf(pFOutFile, "%.16g %.16g\n", OutMesh.pVVert(iV)->dX(), OutMesh.pVVert(iV)->dY());


  for (int iC = 0; iC < OutMesh.iNumCells(); iC++) {
    Cell* pC = OutMesh.pCCell(iC);
    for (int iV = 0; iV < pC->iNumVerts(); iV++) {
      fprintf(pFOutFile, "%d ", OutMesh.iVertIndex(pC->pVVert(iV)));
    }
    /*cells: verts*/
    fprintf(pFOutFile, "\n");
  }

  fclose(pFOutFile);

  // Need to write the bdry data
  if (iOrder > 2) {
    sprintf(strFileName, "%s%s.bpts", strBaseFileName, strExtraFileSuffix);

    vMessage(1,"Opening boundary point output file %s\n", strFileName);
    pFOutFile = fopen(strFileName, "w");
    if (NULL == pFOutFile)
      vFatalError("Couldn't open boundary point output file for writing",
		  "2d mesh output");

    int iNumPoints = iOrder+1;
    double adS[iNumPoints];
    for (int ii = 0; ii < iNumPoints; ii++) {
      adS[ii] = double(ii)/(iNumPoints-1);
    }

    for (int iBF = 0; iBF < OutMesh.iNumBdryFaces(); iBF++) {
      BFace* pBF = OutMesh.pBFBFace(iBF);
      assert((pBF->eType() == Cell::eBdryEdge) ||
	     (pBF->eType() == Cell::eIntBdryEdge));	
	
      BdryPatch2D *pBP = dynamic_cast<BdryPatch2D*>(pBF->pPatchPointer());
      
      double adPointLoc[2];
      
      const double* const adVertLocA = pBF->pVVert(0)->adCoords();
      const double* const adVertLocB = pBF->pVVert(1)->adCoords();

      fprintf(pFOutFile, "%d\t%d\t%d\t%s\n", iBF,
	      OutMesh.iVertIndex(pBF->pVVert(0)),
	      OutMesh.iVertIndex(pBF->pVVert(1)),
	      pBF->iBdryCond() == 1 ? "wall" : "far field");
      
      for (int ii = 0; ii < iNumPoints; ii++) {
 	pBP->vGeodesicDistanceRatioBetweenTwoPoints(adS[ii], adVertLocA,
						    adVertLocB, adPointLoc);
	fprintf(pFOutFile, " %.16g %.16g\n", adPointLoc[0], adPointLoc[1]);
      }
      fprintf(pFOutFile, "\n");
    }
  }

  fclose(pFOutFile);

  SUMAA_LOG_EVENT_END(OUTPUT);
}
