using namespace AMDiS;

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


template<typename ProblemType>
void BaseProblem<ProblemType>::initialize(Flag initFlag,
					  ProblemStat *adoptProblem,
					  Flag adoptFlag)
{ FUNCNAME("BaseProblem::initialize()");
  
  prob->initialize(initFlag, adoptProblem, adoptFlag);
  dim = getMesh()->getDim();
}


template<typename ProblemType>
Flag BaseProblem<ProblemType>::initDataFromFile(AdaptInfo *adaptInfo)
{ 
  FUNCNAME("BaseProblem::initDataFromFile()");
  
  Flag initFlag;
  bool readDataFromFile = false, readArhFiles = false, readDatFiles = 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") {
    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 (size_t 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 (size_t 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<DOFVector<double>*> solutions;
      std::vector<std::string> names;
      for (size_t i = 0; i < prob->getNumComponents(); i++) {
	solutions.push_back(prob->getSolution()->getDOFVector(i));
	names.push_back(prob->getComponentName(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!"));
  }
  
  initFlag.setFlag(DATA_ADOPTED);
  initFlag.setFlag(MESH_ADOPTED);

  return initFlag;
};


template<typename ProblemType>
void BaseProblem<ProblemType>::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>
Flag BaseProblem<ProblemType>::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->getSolver()->setMultipleRhs(false);
    prob->buildAfterCoarsen(adaptInfo, markFlag, true, true);
  }
  #endif

  if (toDo.isSet(SOLVE)) {
    if (secureIteration) {
      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>
void BaseProblem<ProblemType>::endIteration(AdaptInfo *adaptInfo) 
{ FUNCNAME("BaseProblem::endIteration()");
	
  MSG("\n");
  MSG(("[[ end of <"+name+"> iteration ]]\n").c_str());
};


template<typename ProblemType>
void BaseProblem<ProblemType>::closeTimestep(AdaptInfo *adaptInfo)
{ FUNCNAME("BaseProblem::closeTimestep()");

  int outputPeriod = 1;
  Initfile::get("user parameter->write every ith timestep", outputPeriod);
  if (adaptInfo->getTimestepNumber() % outputPeriod == 0
    && adaptInfo->getStartTime() < adaptInfo->getEndTime()) {
    writeFiles(adaptInfo, false);
  }
};



template<typename ProblemType>
void BaseProblem<ProblemType>::addTimeOperator(ProblemStat *prob, int i, int j)
{ FUNCNAME("BaseProblem::addTimeOperator()");

  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>
void BaseProblem<ProblemType>::addTimeOperator(RosenbrockStationary *prob, int i, int j)
{ FUNCNAME("BaseProblem::addTimeOperator()");

  prob->addTimeOperator(i,j);
};