#include "ProblemImplicit.h"

namespace AMDiS {
  namespace calcspace {
    void readDofvec(std::istream& in, DOFVector< double >* vec, Mesh* mesh) 
    {
      long size = mesh->getNumberOfVertices();
      std::string buf;
      getline(in, buf);
      while(buf != "vertex values: solution") getline(in, buf);
      for(long i = 0; i < size; i++)
      {
        in >> buf;
        (*vec)[i] = atof(buf.c_str());
      }
    }

    void readR(std::istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, r, mesh);
      std::vector< double >& vecR = r->getVector();
      std::vector< double >& vecPhi1 = phi1->getVector();
      std::vector< double >& vecPhi2 = phi2->getVector();
      std::vector< double >& vecLevelSet = levelSet->getVector();

      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("something went wrong\n");
      for (int i = vecR.size()-1; i >= 0; --i) 
      {
        vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi2[i] = Phi2(vecR[i], eps);
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

    void readPhi1(istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi1, mesh);
      std::vector< double >& vecR = r->getVector();
      std::vector< double >& vecPhi1 = phi1->getVector();
      std::vector< double >& vecPhi2 = phi2->getVector();
      std::vector< double >& vecLevelSet = levelSet->getVector();

      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("something went wrong\n");
      for (int i = vecR.size()-1; i>= 0; --i) 
      {
        vecR[i] = Phi1ToR(vecPhi1[i], eps);
        //vecPhi2[i] = Phi2(vecR[i], eps);
        vecPhi2[i] = 1 - vecPhi1[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }

    void readPhi2(istream& inStream, double eps, 
		    DOFVector< double >* r, DOFVector< double >* phi1, 
		    DOFVector< double >* phi2, DOFVector< double >* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi2, mesh);
      std::vector< double >& vecR = r->getVector();
      std::vector< double >& vecPhi1 = phi1->getVector();
      std::vector< double >& vecPhi2 = phi2->getVector();
      std::vector< double >& vecLevelSet = levelSet->getVector();

      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("something went wrong\n");
      for (int i = vecR.size()-1; i >= 0; --i) 
      {
        vecR[i] = Phi2ToR(vecPhi2[i], eps);
        //vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi1[i] = 1 - vecPhi2[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }


  }

  bool ProblemImplicitScal::createImplicitMesh() 
  {
    std::string meshFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;

    std::string dofFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

    double eps(-1);
    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
    if (eps<0)
      return false;

    int serType(-1);
    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
    if (serType<0)
      return false;

    TEST_EXIT(mesh != NULL)("the mesh was not created\n");

    std::ifstream inStream(meshFilename.c_str());
    mesh->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin

    createFESpace();

    r = new DOFVector< double >(getFESpace(), "r");
    phi1 = new DOFVector<double>(getFESpace(), "phi1");
    phi2 = new DOFVector<double>(getFESpace(), "phi2");
    levelSet = new DOFVector< double >(getFESpace(), "levelSet");
    inStream.open(dofFilename.c_str());
    switch (serType) {
    case 0:
      calcspace::readR(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 1:
      calcspace::readPhi1(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 2:
      calcspace::readPhi2(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

  void ProblemImplicitScal::initialize(Flag initFlag, ProblemScal* adaptProblem, Flag adoptFlag)
  {
     FUNCNAME("ImplicitScal::initialize()");
     ProblemScal::initialize(CREATE_MESH);
     createImplicitMesh();
     initFlag = initFlag & ~CREATE_MESH;
     ProblemScal::initialize(initFlag);
  }

  bool ProblemImplicitVec::createImplicitMesh()
  {
    //check each mesh if it's an implicit one
    r.resize(meshes.size());
    phi1.resize(meshes.size());
    phi2.resize(meshes.size());
    levelSet.resize(meshes.size());
    for(int i=0;i<meshes.size();++i) {
      createImplicitMesh(i);
    }
  }

  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
    std::string path=name + "->implicit mesh[" + boost::lexical_cast< std::string >(p) + "]->";
    std::string meshFilename("");
    GET_PARAMETER(0, path + "mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;
    std::string dofFilename("");
    GET_PARAMETER(0,path + "dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

    double eps(-1);
    GET_PARAMETER(0, path + "eps", "%d", &eps);
    if(eps < 0)
        return false;
    int serType(-1);
    GET_PARAMETER(0, path + "type", "%d", &serType);
    if (serType < 0)
      return false;
    TEST_EXIT(meshes[p] != NULL)("the mesh was not created\n");
    std::ifstream inStream(meshFilename.c_str());
    meshes[p]->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin
    createFESpace(NULL);
    r[p] = new DOFVector< double >(getFESpace(p), "r");
    phi1[p] = new DOFVector<double>(getFESpace(p), "phi1");
    phi2[p] = new DOFVector<double>(getFESpace(p), "phi2");
    levelSet[p] = new DOFVector< double >(getFESpace(p), "levelSet");
    inStream.open(dofFilename.c_str());
    switch (serType) 
    {
      case 0:
        calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
        break;
      case 1:
	calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
        break;
      case 2:
	calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
        break;
      default:
        WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, Flag adoptFlag)
  {
     FUNCNAME("ImplicitScal::initialize()");
     ProblemVec::initialize(CREATE_MESH);
     createImplicitMesh();
     initFlag = initFlag & ~CREATE_MESH;
     ProblemVec::initialize(initFlag);    
  }
}