using namespace AMDiS;

template<typename ProblemType, bool safeguard>
BaseProblem<ProblemType, safeguard>::BaseProblem(const std::string &name_, bool createProblem) 
  : ProblemInstatBase(name_,NULL),
    prob(NULL),
    dim(1),
    dow(1),
    nTimesteps(-1),
    oldMeshChangeIdx(0),
    oldTimestep(0.0)
{
  // create basic problems
  if (createProblem)
    prob = new ProblemType(name + "->space");
  
  dow = Global::getGeo(WORLD);
}


template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::initialize(Flag initFlag,
					  ProblemStat *adoptProblem,
					  Flag adoptFlag)
{ FUNCNAME("BaseProblem::initialize()");

  TEST_EXIT(prob)("Problem must be created! Add flag createProblem=true to constructor!\n");
  
  prob->initialize(initFlag, adoptProblem, adoptFlag);
  dim = getMesh()->getDim();
}


template<typename ProblemType, bool safeguard>
Flag BaseProblem<ProblemType, safeguard>::initDataFromFile(AdaptInfo *adaptInfo)
{ 
  FUNCNAME("BaseProblem::initDataFromFile()");
  
  Flag initFlag;
  bool readDataFromFile = false;
  Initfile::get(name + "->read data from file", readDataFromFile, 2);
  if (!readDataFromFile)
    return initFlag;
  
  std::string readFormat = "arh";
  Initfile::get(name + "->read format", readFormat, 2);
  if (readFormat != "arh" && readFormat != "dat" && readFormat != "vtu" && readFormat != "multi-vtu") {
    WARNING("You can not read data from formats other than .arh, .dat or .vtu! The .arh-format is selected.\n");
  }
  
  // read data and mesh from arh-files/dat-files
  MSG("read data from file...\n");
  if (readFormat == "arh") 
  {
    std::string filename = "";
    Initfile::get(name + "->value file", filename);
    if (filename.size() == 0) return initFlag;
    if (!file_exists(filename))
      throw(std::runtime_error("The file '" + filename + "' does not exist!"));
    
    std::vector<DOFVector<double>*> solutions;
    for (int i = 0; i < prob->getNumComponents(); i++)
      solutions.push_back(prob->getSolution()->getDOFVector(i));
    ArhReader::read(filename, prob->getMesh(), solutions);
  }
  else if (readFormat == "dat") 
  {
    bool preserveMacroFileInfo = false;
    Parameters::get(prob->getMesh()->getName() + "->preserve macroFileInfo", preserveMacroFileInfo);
    if (prob->getMesh()->getMacroFileInfo() == NULL || !preserveMacroFileInfo)
      throw(std::runtime_error("Yout have to set the flag 'mesh_name->preserve macroFileInfo' to 'true' ("+boost::lexical_cast<std::string>(preserveMacroFileInfo)+"), in order to read .dat-files"));
    
    std::string filename = "";
    for (int i = 0; i < prob->getNumComponents(); i++) {
      Parameters::get(name + "->value file["+boost::lexical_cast<std::string>(i)+"]",filename);
      if (!file_exists(filename))
	throw(std::runtime_error("The file '" + filename + "'does not exist!"));
      ValueReader::readValue(filename,prob->getMesh(),prob->getSolution()->getDOFVector(i),prob->getMesh()->getMacroFileInfo());
    }
  }
  else if (readFormat == "vtu") 
  {
    std::vector<std::string> filenames;
    Initfile::get(name + "->value file", filenames);
    if (filenames.size() == 0)
      return initFlag;
    int arh_idx = -1, vtu_idx = -1;
    for (size_t i = 0; i < filenames.size(); i++) {
      if (!file_exists(filenames[i]))
	throw(std::runtime_error("The file '" + filenames[i] + "' does not exist!"));
      if (filenames[i].find(".vtu") != std::string::npos)
	vtu_idx = i;
      else if (filenames[i].find(".arh") != std::string::npos)
	arh_idx = i;
      else
	throw(std::runtime_error("The file '" + filenames[i] + "' must have the extension .arh or .vtu!"));
    }
    if (arh_idx >= 0) {
      ArhReader::read(filenames[arh_idx], prob->getMesh());
    }
    if (vtu_idx >= 0) {      
      std::vector<int> readComponents;
      Initfile::get(name + "->read components", readComponents);
      if (readComponents.size() == 0)
	for (int i = 0; i < prob->getNumComponents(); i++)
	  readComponents.push_back(i);
	
      std::vector<DOFVector<double>*> solutions;
      std::vector<std::string> names;
      for (size_t i = 0; i < readComponents.size(); i++) {
	solutions.push_back(prob->getSolution()->getDOFVector(readComponents[i]));
	names.push_back(prob->getComponentName(readComponents[i]));
      }
      VtuReader::readValue(filenames[vtu_idx], prob->getMesh(), solutions, names);
    } else
      throw(std::runtime_error("You have to specify a .vtu file!"));
  }
  else if (readFormat == "multi-vtu") 
  {
    size_t numFiles = 0;
    Initfile::get(name + "->number of files", numFiles);
    for (size_t n = 0; n < numFiles; n++) {
      std::vector<std::string> filenames;
      Initfile::get(name + "->value file[" + boost::lexical_cast<std::string>(n) + "]", filenames);
      if (filenames.size() == 0)
	return initFlag;
      int arh_idx = -1, vtu_idx = -1;
      for (size_t i = 0; i < filenames.size(); i++) {
	if (!file_exists(filenames[i]))
	  throw(std::runtime_error("The file '" + filenames[i] + "' does not exist!"));
	if (filenames[i].find(".vtu") != std::string::npos)
	  vtu_idx = i;
	else if (filenames[i].find(".arh") != std::string::npos)
	  arh_idx = i;
	else
	  throw(std::runtime_error("The file '" + filenames[i] + "' must have the extension .arh or .vtu!"));
      }
      if (arh_idx >= 0) {
	ArhReader::read(filenames[arh_idx], prob->getMesh());
      }
      if (vtu_idx >= 0) {      
	std::vector<int> readComponents;
	Initfile::get(name + "->read components[" + boost::lexical_cast<std::string>(n) + "]", readComponents);
	if (readComponents.size() == 0)
	  for (int i = 0; i < prob->getNumComponents(); i++)
	    readComponents.push_back(i);
	  
	std::vector<DOFVector<double>*> solutions;
	std::vector<std::string> names;
	for (size_t i = 0; i < readComponents.size(); i++) {
	  solutions.push_back(prob->getSolution()->getDOFVector(readComponents[i]));
	  names.push_back(prob->getComponentName(readComponents[i]));
	}
	VtuReader::readValue(filenames[vtu_idx], prob->getMesh(), solutions, names);
      } else
	throw(std::runtime_error("You have to specify a .vtu file!"));
    }
  }
  else {
    throw(std::runtime_error("Parameter 'read data from file' set to 'true', but no parameter 'read format' specified!"));
  }
  
  bool readPvdFromFile = false;
  Initfile::get(name + "->read pvd from file", readPvdFromFile, 2);
  if (readPvdFromFile) {
    vector<FileWriterInterface*>& fileWriters = prob->getFileWriterList();
    for (size_t i = 0; i < fileWriters.size(); i++) {
      FileWriterTemplated<double>* fileWriter = dynamic_cast<FileWriterTemplated<double>*>(fileWriters[i]);
      if (fileWriter)
	initFileWriterFromFile(adaptInfo, *fileWriter);
    }
  }
  
  initFlag.setFlag(DATA_ADOPTED);
  initFlag.setFlag(MESH_ADOPTED);

  return initFlag;
}


template<typename ProblemType, bool safeguard>
template<typename T>
void BaseProblem<ProblemType, safeguard>::initFileWriterFromFile(AdaptInfo* adaptInfo, FileWriterTemplated<T>& fileWriter, bool keep_all)
{
  using namespace pugi;
      
  if(!boost::filesystem::exists(fileWriter.getFilename() + ".pvd"))
    return;
  
  xml_document vtu;
  if(!vtu.load_file((fileWriter.getFilename() + ".pvd").c_str()))
    throw(std::runtime_error("Could not load pvd file! Error in xml structure."));

  xml_node VTKFile = vtu.child("VTKFile");
  xml_node Collection = VTKFile.child("Collection");
  
  std::vector<std::pair<double, std::string> >& paraviewAnimationFrames = fileWriter.getParaviewAnimationFrames();

  for (xml_node DataSet = Collection.child("DataSet"); DataSet; DataSet = DataSet.next_sibling("DataSet")) {
    std::string timestepStr = DataSet.attribute("timestep").value();
    std::string fileStr = DataSet.attribute("file").value();
    double time = boost::lexical_cast<double>(timestepStr);
    if (keep_all || time <= adaptInfo->getStartTime())
      paraviewAnimationFrames.push_back(std::make_pair(time, fileStr));
  }
}
  

template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::beginIteration(AdaptInfo *adaptInfo) 
{ FUNCNAME("BaseProblem::beginIteration()");
	
  MSG("\n");
  MSG(("[[ <"+name+"> iteration ]]\n").c_str());
    
  // assemble some blocks only once, if some conditions are fullfilled
  #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  bool fixedMatrixCell = false;
  #else
  bool fixedMatrixCell = adaptInfo->getTimestepNumber() > 0
    && std::abs(adaptInfo->getTimestep() - oldTimestep) < DBL_TOL
    && oldMeshChangeIdx == getMesh()->getChangeIndex();
  #endif
  for(size_t i = 0; i < fixedMatrixTimestep.size(); ++i) {
    prob->setAssembleMatrixOnlyOnce(
      fixedMatrixTimestep[i].first, 
      fixedMatrixTimestep[i].second,
      fixedMatrixCell);
  }
}


template<typename ProblemType, bool safeguard>
Flag BaseProblem<ProblemType, safeguard>::oneIteration(AdaptInfo *adaptInfo, 
                                Flag toDo) 
{ 
  FUNCNAME("BaseProblem::oneIteration()");

  Flag flag, markFlag;
  
  if (toDo.isSet(MARK))
    markFlag = prob->markElements(adaptInfo);
  else
    markFlag = 3;

  // refine
  if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
    flag = prob->refineMesh(adaptInfo);  
  // coarsen
  if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
    flag |= prob->coarsenMesh(adaptInfo);

  #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  MPI::COMM_WORLD.Barrier();
  #endif

#ifndef NONLIN_PROBLEM
  if (toDo.isSet(BUILD))
    prob->buildAfterCoarsen(adaptInfo, markFlag, true, true);
#endif

  if (toDo.isSet(SOLVE)) {
    if (safeguard) {
      try {
	prob->solve(adaptInfo);
      } catch(std::exception &e) {
	WARNING(("Could not solve system. Simulation will be stoped here! ERROR: "
	  + std::string(e.what()) + "\n").c_str());
	adaptInfo->setTime(adaptInfo->getEndTime());
	adaptInfo->setTimestepNumber(adaptInfo->getNumberOfTimesteps());
	return flag;
      } catch(...) {
	WARNING("Could not solve system. Simulation will be stoped here! Unknown ERROR\n");
	adaptInfo->setTime(adaptInfo->getEndTime());
	adaptInfo->setTimestepNumber(adaptInfo->getNumberOfTimesteps());
	return flag;
      }
    } else {
      prob->solve(adaptInfo);
    }
  }

  if (toDo.isSet(ESTIMATE))
    prob->estimate(adaptInfo);

  oldTimestep = adaptInfo->getTimestep();
  oldMeshChangeIdx = getMesh()->getChangeIndex();

  return flag;
}


template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::endIteration(AdaptInfo *adaptInfo) 
{ FUNCNAME("BaseProblem::endIteration()");
	
  MSG("\n");
  MSG(("[[ end of <"+name+"> iteration ]]\n").c_str());
}


template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::closeTimestep(AdaptInfo *adaptInfo)
{
  writeFiles(adaptInfo, false);
}


template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::addTimeOperator(ProblemStat *prob, int i, int j)
{
  Operator *op = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
    op->addZeroOrderTerm(new Simple_ZOT);
  Operator *opRhs = new Operator(prob->getFeSpace(i));
    opRhs->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(j), NULL));

  prob->addMatrixOperator(op, i, j, getInvTau(), getInvTau());
  prob->addVectorOperator(opRhs, i, getInvTau(), getInvTau());
}


template<typename ProblemType, bool safeguard>
void BaseProblem<ProblemType, safeguard>::addTimeOperator(RosenbrockStationary *prob, int i, int j)
{
  prob->addTimeOperator(i,j);
}