FireDomain.cpp

/**
 * @file FireDomain.cpp
 * @brief Implements the methods of the FireDomain class
 * @copyright Copyright (C) 2025 ForeFire, Fire Team, SPE, CNRS/Universita di Corsica.
 * @license This program is free software; See LICENSE file for details. (See LICENSE file).
 * @author Jean‑Baptiste Filippi — 2025
 */

 #include "FireDomain.h"
 #include "BurningMapLayer.h"
 
 #include <sys/stat.h>
 #include "RosLayer.h"
 
 namespace libforefire{
 //  On garde 1 domaine pour faire le parallele - domain 1
 //  
 //  
 //  
 //  Test - Domain 1 est celui qui a les fronts - il propage et il a des ghosts de sous domaines.
 //  Il faut qu'in conaisse la taille totale..Il lit ça dans le NC qui a la taille du domaine physique 
 // attention /etc/beegfs/beegfs-client.conf tuneRemoteFSync=true
 //  (pas de halos externes, faire attention  aux points de flux)
 //  Domain 1 va lire des fichiers de vents U/V a tous les updates
 //  Domain N va écrire des BMaps à tous les updates
 //  Recuperer les extents de tous les sous domaines... 
 //  On fait plus de chaines....
 
	 // Static variables
	 long ForeFireAtom::instanceNRCount = 0;
 
	 const double FireDomain::endChain = -1.;
	 const double FireDomain::endCom = -10.;
	 const double FireDomain::noCom = -100.;
 
 
	 const string FireDomain::altitude = "altitude";
	 const FFPoint FireDomain::outPoint
	 = FFPoint(numeric_limits<double>::infinity(), numeric_limits<double>::infinity(),0);
 
	 bool FireDomain::outputs = false;
	 bool FireDomain::commandOutputs = false;
	 bool FireDomain::recycleNodes = false;
	 bool FireDomain::recycleFronts = false;
	 
	 std::list<FireDomain::distributedDomainInfo*> FireDomain::parallelDispatchDomains;
	 std::list<FireNode*> FireDomain::createdNodes;
	 std::list<FireNode*> FireDomain::trashNodes;
	 std::list<FireFront*> FireDomain::trashFronts;
 
	 FireFrontData* FireDomain::mainFrontBackup;
	 size_t FireDomain::atmoIterNumber = 0;
	 FireDomain::FrontDepthScheme FireDomain::fdScheme = FireDomain::normalDir;
 
	 const size_t FireDomain::NUM_MAX_PROPMODELS;
	 const size_t FireDomain::NUM_MAX_FLUXMODELS;
 
	 PropagationModel* FireDomain::propModelsTable[FireDomain::NUM_MAX_PROPMODELS];
	 FluxModel* FireDomain::fluxModelsTable[FireDomain::NUM_MAX_FLUXMODELS];
 
	 FireDomain::PropModelMap& FireDomain::prop_instantiatorMap(){
		 static FireDomain::PropModelMap* inst = new FireDomain::PropModelMap;
		 return *inst;
	 }
 
	 FireDomain::FluxModelMap& FireDomain::flux_instantiatorMap(){
		 static FireDomain::FluxModelMap* inst = new FireDomain::FluxModelMap;
		 return *inst;
	 }
 
	 // Default constructor
	 FireDomain::FireDomain() : ForeFireAtom(0.) {
		 cout<<"Trying to instantiate an empty FireDomain, not relevant"<<endl;
	 }	
 
	 FireDomain::FireDomain(const double& t
							, FFPoint& swc, FFPoint& nec)
	 : ForeFireAtom(t), SWCorner(swc), NECorner(nec) {
		 isFireActive = false;
		 refLatitude = 0.;
		 refLongitude = 0.;
		 atmosphericCoupling = false;
		 propagationSpeedAdjustmentFactor =1;
		 //parallelDispatchDomains.size() = 0;
		 params = SimulationParameters::GetInstance();
		 numIterationAtmoModel = 0;
		 // Maximum time-step for Firenodes is not constrained
		 dtMax = numeric_limits<double>::infinity();
		 double readRes = 0;
 
 
		 if (params->getParameter("runmode") == "masterMNH"){
 
			  #ifdef MPI_COUPLING
			  #else
				 readMultiDomainMetadata();
			 #endif

			 atmosphericCoupling = true;
	 
			 if(parallelDispatchDomains.size() >0){
 
			 
			 distributedDomainInfo* fit = parallelDispatchDomains.front();
			 double Mswx = fit->SWCorner->getX();
			 double Mswy = fit->SWCorner->getY();
			 double Mnex = fit->NECorner->getX();
			 double Mney = fit->NECorner->getY();
  
			 readRes = (Mnex-Mswx)/fit->atmoNX;
			 list<distributedDomainInfo*>::iterator it;
			 for (it = parallelDispatchDomains.begin(); it != parallelDispatchDomains.end(); it++)
				 {
				  
						 Mswx = std::min(Mswx,(*it)->SWCorner->getX());
						 Mswy = std::min(Mswy,(*it)->SWCorner->getY());
						 Mnex = std::max(Mnex,(*it)->NECorner->getX());
						 Mney = std::max(Mney,(*it)->NECorner->getY());
			  
			 }
			SWCorner.setLoc(Mswx,Mswy);
			NECorner.setLoc(Mnex,Mney); 

			atmoNX = (Mnex-Mswx)/readRes;
			atmoNY = (Mney-Mswy)/readRes;

			params->setSize("atmoNX",atmoNX);
			params->setSize("atmoNY",atmoNY);
			
 
			 for (it = parallelDispatchDomains.begin(); it != parallelDispatchDomains.end(); it++)
				 {
				  
						 (*it)->refNX = size_t(((*it)->SWCorner->getX()-Mswx)/readRes);
						 (*it)->refNY = size_t(((*it)->SWCorner->getY()-Mswy)/readRes);
					 //	cout<<(*it)->ID<<" HAVE A RESOLUTION OF "<< readRes <<" location "<<(*it)->refNX<<":"<<(*it)->refNY<<" coords:"<<(*it)->SWCorner->print() << ":"<<(*it)->NECorner->print()<<endl;
						 
			 }
					 
			 }else{
				 cout<<"ERROR big problem master mode but 0 domains read, rerun for domain file caching"<<endl;
				 atmoNX = params->getSize("atmoNX");
				 atmoNY = params->getSize("atmoNY");
	  
			 }
 
 
	 
		 }else{
			 params->setParameter("runmode", "standalone");
			 atmoNX = params->getSize("atmoNX");
			 atmoNY = params->getSize("atmoNY");
		 }
		 atmoNX = params->getSize("atmoNX");
		 atmoNY = params->getSize("atmoNY");
		 atmoNZ = params->getSize("atmoNZ");
		 // computing the cells' mesh
		 double cellsMeshX[atmoNX+1];
		 double cellsMeshY[atmoNY+1];
		 cellsMeshX[0] = SWCorner.getX();
		 cellsMeshY[0] = SWCorner.getY();
		 double dx, dy;
		 dx = (NECorner.getX()-SWCorner.getX())/atmoNX;
		 dy = (NECorner.getY()-SWCorner.getY())/atmoNY;
		 for ( size_t i = 1; i < atmoNX; i++ ) {
			 cellsMeshX[i] = cellsMeshX[i-1] + dx;
		 }
		 for ( size_t j = 1; j < atmoNY; j++ ) {
			 cellsMeshY[j] = cellsMeshY[j-1] + dy;
		 }
		 cellsMeshX[atmoNX] = NECorner.getX();
		 cellsMeshY[atmoNY] = NECorner.getY();
 
		 // simulation won't be run in parallel
		 params->setInt("parallel", 0);
 
		 // simulation won't be be run with mnh
		
 
		 // reference time
		 int year = params->getInt("refYear");
		 int yday = params->getInt("refDay");
		 
		 if (yday == 0)
		 {
			 year = params->getInt("year");
			 int month = params->getInt("month");
			 int day = params->getInt("day");
			 yday = getDayNumber(year, month, day);
		 }
		 
		 // Common initialization for all constructors
		 commonInitialization(cellsMeshX,cellsMeshY , year, yday);
 
	 }
 
	 FireDomain::FireDomain(const int& mpirank
							, const int& year, const int& month
							, const int& day, const double& t
							, const double& lat, const double& lon
							, const int& mdimx, const double* meshx
							, const int& mdimy, const double* meshy
							, const int& mdimz, const double& dt)
	 : ForeFireAtom(t), refLatitude(lat), refLongitude(lon) {
 
		 getNewID(mpirank);
		 isFireActive = false;
		 numIterationAtmoModel = 0;
		 params = SimulationParameters::GetInstance();
 
		 // Maximum time-step for Firenodes is constrained by the atmospheric model
		 dtMax = dt;
	 
		 /* simulation is run with mnh */
		 params->setParameter("runmode", "coupled");
		 params->setInt("parallel", 1);
		 params->setInt("mpirank", mpirank);
 
		 atmosphericCoupling = true;
 
		 /* Definition of the mesh */
		 atmoNX = (size_t) mdimx;
		 params->setSize("atmoNX", atmoNX);
		 atmoNY = (size_t) mdimy;
		 params->setSize("atmoNY", atmoNY);
		 atmoNZ = (size_t) mdimz;
		 params->setSize("atmoNZ", atmoNZ);
 
	 string ffOutputsPattern(params->getParameter("caseDirectory")+'/'+params->getParameter("fireOutputDirectory")+'/'+params->getParameter("outputFiles")+"."+params->getParameter("mpirank"));
	 SimulationParameters::GetInstance()->setParameter("ffOutputsPattern", ffOutputsPattern);
 
	 
 
		 // computing the position of the physical and numerical corners
		 SWCorner = FFPoint(meshx[0], meshy[0],0);
		 NECorner = FFPoint(2.*meshx[atmoNX-1]-meshx[atmoNX-2]
							, 2.*meshy[atmoNY-1]-meshy[atmoNY-2],0);
 
		 bool dumpDomainInfo = true;
		 #ifdef MPI_COUPLING
			 dumpDomainInfo = false;
		 #endif
		 string ffDomPattern(params->getParameter("caseDirectory")+'/'+params->getParameter("PPath")+'/'+params->getParameter("mpirank")+".domainData");
		 if (dumpDomainInfo){
			 ofstream FileOut(ffDomPattern.c_str(), ios_base::binary);
 
			 FileOut.write(reinterpret_cast<const char*>(&atmoNX), sizeof(size_t));
			 FileOut.write(reinterpret_cast<const char*>(&atmoNY), sizeof(size_t));
			 FileOut.write(reinterpret_cast<const char*>(&SWCorner.x), sizeof(double));
			 FileOut.write(reinterpret_cast<const char*>(&SWCorner.y), sizeof(double));
			 FileOut.write(reinterpret_cast<const char*>(&NECorner.x), sizeof(double));
			 FileOut.write(reinterpret_cast<const char*>(&NECorner.y), sizeof(double));
			 
			 FileOut.flush();   
			 FileOut.rdbuf()->pubsync(); 		
			 FileOut.close();
		 }
 
		 // computing the cells' mesh
		 double cellsMeshX[atmoNX+1];
		 double cellsMeshY[atmoNY+1];
		 for ( size_t i = 0; i < atmoNX; i++ ) {
			 cellsMeshX[i] = meshx[i];
		 }
		 cellsMeshX[atmoNX] = NECornerX();
		 for ( size_t j = 0; j < atmoNY; j++ ) {
			 cellsMeshY[j] = meshy[j];
		 }
		 cellsMeshY[atmoNY] = NECornerY();
 
		 int yday = getDayNumber(year, month, day);
		 
		 // Common initialization for all constructors
		 commonInitialization(cellsMeshX,cellsMeshY, year, yday);
 
	 }
 
	 FireDomain::~FireDomain() {
		// Deleting nodes
		while (!createdNodes.empty()) {
			if (createdNodes.back()) {
				delete createdNodes.back();
			}
			createdNodes.pop_back();
		}
	
		// Deleting fronts
		if (domainFront) {
			delete domainFront;
			domainFront = nullptr;
		}
		while (!trashFronts.empty()) {
			if (trashFronts.back()) {
				delete trashFronts.back();
			}
			trashFronts.pop_back();
		}
	
		// Deleting cells
		if (cells) {
			for (size_t i = 0; i < atmoNX; ++i) {
				if (cells[i]) {
					delete[] cells[i];
				}
			}
			delete[] cells;
			cells = nullptr;
		}
	
		if (trashCell) {
			delete trashCell;
			trashCell = nullptr;
		}
	
		// Deleting frontiers
		while (!frontiers.empty()) {
			if (frontiers.back()) {
				delete frontiers.back();
			}
			frontiers.pop_back();
		}
		while (!infrontiers.empty()) {
			if (infrontiers.back()) {
				delete infrontiers.back();
			}
			infrontiers.pop_back();
		}
	
		// Deleting data broker and layers
		if (dataBroker) {
			delete dataBroker;
			dataBroker = nullptr;
		}
		if (propagativeLayer) {
			//delete propagativeLayer;
			propagativeLayer = nullptr;
		}
	
		if ( mainFrontBackup != 0 ) delete mainFrontBackup;
	}
 
	 void FireDomain::backupState(){
		 /* Backup of the simulation */
		 if ( mainFrontBackup != 0 ) delete mainFrontBackup;
		 mainFrontBackup = new FireFrontData(domainFront);
	 }
 inline size_t fileSize(const std::string& name) {
   struct stat buffer;   
	   if(stat(name.c_str(), &buffer) != 0) {
		 return 0;
	 }
	 return buffer.st_size;   
 }
 void FireDomain::setPerimeterResolution(double res) {
	this->perimeterResolution = res;
 }
 void FireDomain::setSpatialIncrement(double inc) {
	 this->spatialIncrement = inc;
 }
 void FireDomain::pushMultiDomainMetadataInList(size_t id, double lastTime, size_t atmoNX, size_t atmoNY, double nswx, double nswy, double nnex, double nney) {
	 distributedDomainInfo *currentInfo = new distributedDomainInfo;
	 currentInfo->ID = id;
	 currentInfo->lastTime = lastTime;
	 currentInfo->atmoNX = atmoNX;
	 currentInfo->atmoNY = atmoNY;
	 currentInfo->SWCorner = new FFPoint(nswx, nswy, 0.);
	 currentInfo->NECorner = new FFPoint(nnex, nney, 0.);
	 
	 parallelDispatchDomains.push_back(currentInfo);
 }
 void FireDomain::readMultiDomainMetadata(){
		 if (getDomainID()!=0) return;
 
		 size_t numReadDom = 1;
		 string ffDomPattern(params->getParameter("caseDirectory")+'/'+params->getParameter("PPath")+'/');
		  size_t domFileSize = fileSize(ffDomPattern+std::to_string(numReadDom)+".domainData");
		 while(domFileSize > 0){
			  numReadDom++;
			 domFileSize = fileSize(ffDomPattern+std::to_string(numReadDom)+".domainData");
		 }
		 numReadDom -=1;
  
		 //parallelDispatchDomains.size() = numReadDom;
 
		 while(numReadDom > 0){
 
			 ifstream FileIn((ffDomPattern+std::to_string(numReadDom)+".domainData").c_str(), ios_base::binary);
 
			 size_t anx = 0;
			 size_t any = 0;
			 double nswx = 0;
			 double nswy = 0;
			 double nnex = 0;
			 double nney = 0;
	 
			 FileIn.read((char *)&anx, sizeof(size_t));
			 FileIn.read((char *)&any, sizeof(size_t));
			 FileIn.read((char *)&nswx, sizeof(double));
			 FileIn.read((char *)&nswy, sizeof(double));
			 FileIn.read((char *)&nnex, sizeof(double));
			 FileIn.read((char *)&nney, sizeof(double));
			 FileIn.close();
			 pushMultiDomainMetadataInList(numReadDom,0,anx,any,nswx,nswy,nnex,nney);
 
			 numReadDom--;
		 }	
	 
	 
	 } 
 
	 void FireDomain::restoreValidState(){
		 /* Restoring a previously saved valid state */
		 if ( mainFrontBackup == 0 ){
			 debugOutput<<getDomainID()<<": PROBLEM, tried to return to a valid state"
			 <<" while no backup state is available"<<endl;
		 } else {
			 debugOutput<<getDomainID()<<": RECONSTRUCTING STATE"<<endl;
		 }
		 /* Trashing the current state */
		 trashFrontsAndNodes();
		 /* Restoring the valid state */
		 mainFrontBackup->reconstructState(domainFront);
	 }
 
	 void FireDomain::trashFrontsAndNodes(){
		 list<FireNode*> toBeTrashed;
		 list<FireNode*>::iterator node;
		 for ( size_t i = 0; i < atmoNX; i++ ) {
			 for ( size_t j = 0; j < atmoNY; j++ ) {
				 // deleting the nodes in the cell
				 toBeTrashed = cells[i][j].fireNodes;
				 for ( node = toBeTrashed.begin();
					  node != toBeTrashed.end(); ++node ){
					 debugOutput<<getDomainID()<<": FireDomain::trashFrontsAndNodes -> ";
					 addToTrashNodes(*node);
				 }
			 }
		 }
		 list<FireFront*> fronts = domainFront->getInnerFronts();
		 list<FireFront*>::iterator front;
		 for ( front = fronts.begin(); front != fronts.end(); ++front ){
			 addToTrashFronts(*front);
		 }
	 }
 
	 // Accessors
	 int FireDomain::getReferenceYear(){
		 return refYear;
	 }
	 int FireDomain::getReferenceDay(){
		 return refDay;
	 }
	 double FireDomain::getSecondsFromReferenceTime(
													const int& year, const int& month
													, const int& day, const int& t){
		 if ( year != refYear ) cout<<"WARNING: change of year "<< year <<" to "<< refYear <<" is not handled in ForeFire data at YY/MM/DD "<<year<<"/"<<month<<"/"<<day<<endl;
		 int iday = getDayNumber(year, month, day);
		 return (iday-refDay)*86400. + t;
	 }
	 list<FireFront*> FireDomain::getMainFronts(){
		 return domainFront->getInnerFronts();
	 }
	 FireFront* FireDomain::getDomainFront(){
		 return domainFront;
	 }
	 TimeTable* FireDomain::getTimeTable(){
		 return schedule;
	 }
	 double FireDomain::getSimulationTime(){
		 if ( schedule->getTime() > this->getTime() ) {
			 return double(int(schedule->getTime()));
		 } else {
			 return this->getTime();
		 }
	 }
 
	 FFPoint& FireDomain::getSWCorner(){
		 return SWCorner;
	 }
	 FFPoint& FireDomain::getNECorner(){
		 return NECorner;
	 }
	 double& FireDomain::SWCornerX(){
		 return SWCorner.getX();
	 }
	 double& FireDomain::SWCornerY(){
		 return SWCorner.getY();
	 }
	 double& FireDomain::NWCornerX(){
		 return NWCorner.getX();
	 }
	 double& FireDomain::NWCornerY(){
		 return NWCorner.getY();
	 }
	 double& FireDomain::NECornerX(){
		 return NECorner.getX();
	 }
	 double& FireDomain::NECornerY(){
		 return NECorner.getY();
	 }
	 double& FireDomain::SECornerX(){
		 return SECorner.getX();
	 }
	 double& FireDomain::SECornerY(){
		 return SECorner.getY();
	 }
 
	 double& FireDomain::getPerimeterResolution(){
		 return perimeterResolution;
	 }
 
	 double& FireDomain::getSpatialIncrement(){
		 return spatialIncrement;
	 }
 
	 double& FireDomain::getMaxTimeStep(){
		 return dtMax;
	 }
 
	 int FireDomain::getNumIterationAtmoModel(){
		 return numIterationAtmoModel;
	 }
	 void FireDomain::increaseNumIterationAtmoModel(){
		 numIterationAtmoModel++;
	 }
 
	 int FireDomain::getDayNumber(const int& year
								  , const int& month, const int& day){
		 bool leap = ( ( year%4 == 0 and year%100 != 0 ) or year%400 == 0 );
		 int febNumDays;
		 leap == true ? febNumDays = 29 : febNumDays = 28;
		 int numDays[11] = {31, febNumDays, 31, 30, 31, 30, 31, 31, 30, 31, 30};
		 int numDay = 0;
		 for ( int imon = 1; imon < month; imon++ ){
			 numDay += numDays[imon-1];
		 }
		 numDay += day;
		 return numDay;
	 }
 
	 // Mutators
	 void FireDomain::setSafeTopologyMode(bool safeMode){
		 safeTopologyMode = safeMode;
	 }
 
	 void FireDomain::setBoundariesFront(FireFront* ff){
		 domainFront = ff;
	 }
	 void FireDomain::setTimeTable(TimeTable* tt){
		 schedule = tt;
	 }
 
	 // input function
	 void FireDomain::input(){
	 }
 
	 // update function
	 void FireDomain::update(){
		 setTime(getUpdateTime());
	 }
 
	 // advance in time function
	 void FireDomain::timeAdvance(){
		 setUpdateTime(numeric_limits<double>::infinity());
	 }
 
	 // Output function
	 void FireDomain::output(){
	 }
 
	 // Adding a newly created Atom to the simulation
	 void FireDomain::addNewAtomToSimulation(ForeFireAtom* atom){
		 schedule->insertBefore(new FFEvent(atom));
	 }
 
	 // Deleting an Atom of the simulation
	 void FireDomain::deleteAtomOfSimulation(ForeFireAtom* atom){
		 schedule->dropAtomEvents(atom);
	 }
 
	 // Searching for an atom in the entire domain
	 FireNode* FireDomain::getFireNodeByID(const long& sid){
		 FireNode* tmpAtom;
		 for ( size_t i=0; i<atmoNX; i++ ){
			 for ( size_t j=0; j<atmoNY; j++ ){
				 tmpAtom = cells[i][j].getFirenodeByID(sid);
				 if ( tmpAtom ) return tmpAtom;
			 }
		 }
		 return 0;
	 }
 
	 FireNode* FireDomain::getFireNodeByID(const double sid){
		 long lsid = getIDfromDouble(sid);
		 return getFireNodeByID(lsid);
	 }
	 
	 FireNode* FireDomain::getFireNodeByIDNeighborCells(
														const double id, FDCell* cell, int range){
		 FireNode* tmpAtom;
		 /* Looking in the current cell */
		 tmpAtom = cell->getFirenodeByID(id);
		 if ( tmpAtom ) return tmpAtom;
		 /* searching in neighboring cells */
		 list<FDCell*> proxCells = getProxCells(cell, range);
		 list<FDCell*>::iterator curCell;
		 for ( curCell = proxCells.begin();
			  curCell != proxCells.end(); ++curCell){
			 tmpAtom = (*curCell)->getFirenodeByID(id);
			 if ( tmpAtom ) return tmpAtom;
		 }
		 // No firenode was found with this ID
		 debugOutput<<getDomainID()<<": "
			 <<"could not find a firenode with Id "
			 <<getDomainID(id)<<", "<<getShortID(id)<<endl;
		 return 0;
	 }
 
	  
 
	 int FireDomain::registerPropagationModelInstantiator(string modelname
														  , PropagationModelInstantiator func){
		 PropModelMap::iterator pmodel = prop_instantiatorMap().find(modelname);
		 if ( pmodel != prop_instantiatorMap().end() ){
			 cout<<"WARNING: propagation model "<<modelname
			 <<" has the same name as a previously"
			 <<" registered propagation model, check your settings"<<endl;
		 }
		 prop_instantiatorMap().insert(make_pair(modelname,func));
		 return 1;
	 }
 
	 void FireDomain::updateFuelTable( string key, double value){
		for ( size_t i = 0; i < NUM_MAX_PROPMODELS; i++ ){
			if(propModelsTable[i] != 0){
 
 
				getDataBroker()->updateFuelValues(propModelsTable[i],key,value);
			}
		}
	 }
	 int FireDomain::registerFluxModelInstantiator(string modelname
												   , FluxModelInstantiator func){
		 FluxModelMap::iterator fmodel = flux_instantiatorMap().find(modelname);
		 if ( fmodel != flux_instantiatorMap().end() ){
			 cout<<"WARNING: flux model "<<modelname
			 <<" has the same name as a previously"
			 <<" registered flux model, check your settings"<<endl;
		 }
		 flux_instantiatorMap().insert(make_pair(modelname,func));
		 return 1;
	 }
 
	 PropagationModel* FireDomain::propModelInstanciation(const int& index, string modelname){
		 PropModelMap::iterator pmodel = prop_instantiatorMap().find(modelname);
		 if ( pmodel != prop_instantiatorMap().end() ){
			 return (pmodel->second)(index, dataBroker);
		 } else {
			 cout<<"ERROR: Propagation model "<<modelname
				 <<" is not recognized, check your configuration !!"<< endl;
			 return 0;
		 }
	 }
 
 
	 FluxModel* FireDomain::fluxModelInstanciation(const int& index, string modelname){
		 FluxModelMap::iterator fmodel;
 
			 
		 fmodel = flux_instantiatorMap().find(modelname);
		 
 
		 if ( fmodel != flux_instantiatorMap().end() ){
 
			 return (fmodel->second)(index, dataBroker);
		 } else {
			 cout<<"ERROR: Flux model "<<modelname
			 <<" is not recognized, check your configuration !!"<< endl;
			 return 0;
		 }
	 }
 
	 void FireDomain::registerPropagationModel(const int& index, PropagationModel* model){
			 // Join all strings in the vector into a single string separated by semicolons
	 std::string properties = "";
	 for(const std::string& prop : model->wantedProperties) {
		 if(properties.empty()) {
			 properties = prop; // First element, add without semicolon
		 } else {
			 properties += ";" + prop; // Subsequent elements, add with semicolon
		 }
	 }
	 
	 // Set the parameter with the joined string
		 params->setParameter(model->getName() + ".keys", properties);
		 //cout<< "loading"<< model->getName()<<index<<endl;
		 propModelsTable[index] = model;
	 }
 
	 void FireDomain::registerFluxModel(const int& index, FluxModel* model){
		  
		 fluxModelsTable[index] = model;
	 }
 
	 bool FireDomain::addPropagativeLayer(string mname){
		 /* searching if there exists a propagation model with associated name */
		 /* affecting it to free index */
		 size_t mindex = getFreePropModelIndex();
		 PropagationModel* model = propModelInstanciation(mindex, mname);
		 if ( model == 0 ) return false;
		 /* Instantiating a flux layer related to this model */
		 propagativeLayer = new PropagativeLayer<double>("ROSlayer", mindex);
		 return true;
	 }
 
	 bool FireDomain::addLayer(string type, string layername, string keyname){
		double timespan = 0;
		size_t  nnx = 1;
		size_t nny = 1;
		size_t nnz = 1;
		size_t nnk = 1;

		double spanx = NECorner.x -SWCorner.x;
		double spany = NECorner.y -SWCorner.y;
 
 
		 if ( type == "BRatio" ){
			 BurningRatioLayer<double>* brlayer = new BurningRatioLayer<double>(layername, SWCorner, NECorner, atmoNX, atmoNY, cells);
			 dataBroker->registerLayer(layername, brlayer);
			 return true;
		 }
		 if ( type == "MaxRos" ){
			 RosLayer<double>* mrlayer = new RosLayer<double>(layername, SWCorner, NECorner, atmoNX, atmoNY, cells);
				 dataBroker->registerLayer(layername, mrlayer);
				 return true;
			 }
 
 
		 if ( !params->isValued(keyname) ){
			 cout << "Unknown parameter "<<keyname<<" for "<<layername << " layer of type "<< type <<", please set parameter"<<endl;
			 return false;
		 }
 

 
		 if ( type == "data" ){
			 double* values = new double[1];
			 values[0] = params->getDouble(keyname);
			 return addScalarLayer( type,  layername,SWCorner.getX(), SWCorner.getY(), getTime(), spanx , spany, timespan,  nnx,	 nny, nnz,  nnk,  values);
		 }else{
			 int* values = new int[1];
			 values[0] = params->getInt(keyname);
			 return addIndexLayer( type,  layername,SWCorner.getX(), SWCorner.getY(), getTime(), spanx , spany, timespan,  nnx,	 nny, nnz,  nnk,  values);
		 }
 
	 }
 
	 bool FireDomain::addScalarLayer(string type, string name, double &x0, double &y0, double& t0, double& width, double& height, double& timespan, size_t& nnx,	size_t& nny, size_t& nnz, size_t& nnk, double* values){
		 FFPoint origin = FFPoint(x0, y0,0);
		 FFPoint span = FFPoint(width, height,0);
 
		 XYZTDataLayer<double>* newLayer = new XYZTDataLayer<double>(name, origin,t0, span, timespan, nnx, nny, nnz, nnk, values);
		 dataBroker->registerLayer(name, newLayer);
 
		 //cout<<"adding scalar Layer "<<name<<" of type "<<type<<" position "<<origin.x<<";"<<origin.y<<" size "<<span.x<< "initial date"<< t0<<" time duration "<<timespan<<endl;
		 if ( type == "windScalDir" ){
			 if ( name == "windU" ){
				 
				 dataBroker->PwindULayer = newLayer;
			 }
			 if ( name == "windV" ){
				 dataBroker->PwindVLayer = newLayer;
			 }
			 }
 
 
		 return true;
	 }
	 bool FireDomain::addIndexLayer(string type, string name, double &x0, double &y0, double& t0, double& width, double& height, double& timespan, size_t& nnx,	size_t& nny, size_t& nnz, size_t& nnk, int* values){
			 FFPoint origin = FFPoint(x0, y0,0);
			 FFPoint span = FFPoint(width, height,0);
 
 
			 if ( type == "flux" ){
				 size_t mindex = values[0];
				 string fmname = name;
				 
				 FluxModel* model = fluxModelInstanciation(mindex, fmname);
				 if ( model != 0 ){
					 /* Instantiating a flux layer related to this model */
					 FluxLayer<double>* newlayer = new FluxLayer<double>(name,SWCorner, NECorner,atmoNX, atmoNY, getCells(), values, origin, t0, span,	timespan, nnx, nny, nnz, nnk);
					 dataBroker->registerFluxLayer(name, newlayer);
					 return true;
				 }
 
			 }
			 if ( type == "propagation" ){
				 size_t mindex = getFreePropModelIndex();
				 PropagationModel* model = propModelInstanciation(mindex, name);
				 if ( model == 0 ) return false;
				 /* Instantiating a prop layer related to this model */
				 propagativeLayer = new PropagativeLayer<double>(name, mindex);
				 return true;
			 }
			 if ( type == "table" ){
 
					 FuelDataLayer<double>* newlayer = new FuelDataLayer<double>(name,	origin, t0, span, timespan, nnx, nny, nnz, nnk,				values);
				 
					 dataBroker->registerLayer(name, newlayer);
			 }
 
			 return false;
		 }
 
	 bool FireDomain::addFluxLayer(string lname){
		 /* searching if there exists a flux model with associated name */
		 /* affecting it to free index */
		
		 // Burning ratio is special as it is derived from the heat flux layer
		 if ( lname == "BRatio" or lname == "Bratio" or lname == "bratio" ){
			 /* Instantiating a burning ratio layer */
			 BurningRatioLayer<double>* brlayer =
			 new BurningRatioLayer<double>(lname,  SWCorner, NECorner,atmoNX, atmoNY, cells);
			 dataBroker->registerLayer(lname, brlayer);
			 return true;
		 }
 
		 // Otherwise, searching for the model in the available ones
		 size_t mindex = getFreeFluxModelIndex();
		 string fmname = "";
		 if ( lname == "heatFlux" ) fmname = "heatFluxBasic";
		 if ( lname == "vaporFlux" ) fmname = "vaporFluxBasic";
		 
		 string defaultModelKey = lname+"DefaultModel";

		 if (params->isValued(defaultModelKey)){
			fmname =  params->getParameter(lname+"DefaultModel");
		 }

		 FluxModel* model = fluxModelInstanciation(mindex, fmname);

 
		 if ( model != 0 ){
			 /* Instantiating a flux layer related to this model */
			 FluxLayer<double>* newlayer = new FluxLayer<double>(lname, SWCorner, NECorner, atmoNX, atmoNY, cells, mindex);
			 registerFluxModel(model->index, model);
			 dataBroker->registerFluxLayer(lname, newlayer);
			 return true;
		 }else{
			 cout<<"ERROR: Flux model "<<fmname
			 <<" is not recognized, check your configuration !!"<< endl;
			 return false;
		 }
		 return false;
	 }
 
	 size_t FireDomain::getFreePropModelIndex(){
		 size_t mindex = NUM_MAX_PROPMODELS - 1;
		 while ( propModelsTable[mindex] != 0 ) mindex--;
		 return mindex;
	 }
 
	 size_t FireDomain::getFreeFluxModelIndex(){
		 size_t mindex = 0;
		 while ( fluxModelsTable[mindex] != 0 ){
			 if(mindex >= NUM_MAX_FLUXMODELS -1){
				 cout<<"ERROR No mor flx models allowed, max:"<< NUM_MAX_FLUXMODELS<<endl;
				 return mindex;
			 }
			 mindex++;
		 }
		 return mindex;
	 }
 
	 double FireDomain::getArrivalTime(FFPoint& loc){
		 if ( !striclyWithinDomain(loc) ) return numeric_limits<double>::infinity();
		 size_t ii = (size_t) ((loc.getX()-SWCornerX())/burningMatrixResX);
		 size_t jj = (size_t) ((loc.getY()-SWCornerY())/burningMatrixResY);
		 return getArrivalTime(ii, jj);
	 }
	 double FireDomain::getMaxSpeed(FFPoint& loc){
		 if ( !striclyWithinDomain(loc) ) return numeric_limits<double>::infinity();
		 size_t ii = (size_t) ((loc.getX()-SWCornerX())/burningMatrixResX);
		 size_t jj = (size_t) ((loc.getY()-SWCornerY())/burningMatrixResY);
		 return getMaxSpeed(ii, jj);
	 }
 
	 double FireDomain::getArrivalTime(const size_t& ii, const size_t& jj){
		 if ( ii > globalBMapSizeX-1 ) return numeric_limits<double>::infinity();
		 if ( jj > globalBMapSizeY-1 ) return numeric_limits<double>::infinity();
		 size_t i = ii/localBMapSizeX;
		 size_t j = jj/localBMapSizeY;
		 return cells[i][j].getArrivalTime(ii%localBMapSizeX,jj%localBMapSizeY);
	 }
 
	 double FireDomain::getMaxSpeed(const size_t& ii, const size_t& jj) {
			 if (ii >= globalBMapSizeX || jj >= globalBMapSizeY) return std::numeric_limits<double>::infinity();
			 double current_time = getArrivalTime(ii, jj);
			 if (current_time == std::numeric_limits<double>::infinity()) return std::numeric_limits<double>::infinity();
 
			 double resolution = params->getDouble("bmapResolution");
			 double grad_x = 0.0, grad_y = 0.0;
 
			 if (ii > 0 && ii < globalBMapSizeX - 1) {
				 double time_left = getArrivalTime(ii - 1, jj);
				 double time_right = getArrivalTime(ii + 1, jj);
				 if (time_left != std::numeric_limits<double>::infinity() && time_right != std::numeric_limits<double>::infinity()) {
					 grad_x = (time_right - time_left);
				 }
			 }
			 if (jj > 0 && jj < globalBMapSizeY - 1) {
				 double time_down = getArrivalTime(ii, jj - 1);
				 double time_up = getArrivalTime(ii, jj + 1);
				 if (time_down != std::numeric_limits<double>::infinity() && time_up != std::numeric_limits<double>::infinity()) {
					 grad_y = (time_up - time_down);
				 }
			 }
 
			 double grad_squared = grad_x * grad_x + grad_y * grad_y;
			 if (grad_squared == 0.0) return std::numeric_limits<double>::infinity(); // Safe division check
			 double grad_mag = (2 * resolution) / sqrt(grad_squared);
			 
			 return grad_mag;
		 }
 
	 void FireDomain::setArrivalTime(const size_t& ii, const size_t& jj, const double& t){
		 // getting the cell where to pixel to set to burning lies
 
		 if ( ii > globalBMapSizeX-1 ) return ;
		 if ( jj > globalBMapSizeY-1 ) return ;
		 if ( ii > globalBMapSizeX-1 ) return ;
		 if ( jj > globalBMapSizeY-1 ) return ;
		 size_t i = ii/localBMapSizeX;
		 size_t j = jj/localBMapSizeY;
		 // setting the arrival time of this pixel
		 cells[i][j].setArrivalTime(ii%localBMapSizeX,jj%localBMapSizeY,t);
	 }
 
	 bool FireDomain::burnCheck(const size_t& ii, const size_t& jj, const double& t){
		 if ( getArrivalTime(ii-1,jj-1) > t ) return true;
		 if ( getArrivalTime(ii-1,jj) > t ) return true;
		 if ( getArrivalTime(ii,jj-1) > t ) return true;
		 if ( getArrivalTime(ii,jj) > t ) return true;
		 return false;
	 }
 
	 void FireDomain::firenodeBurningScan(FireNode* fn){
 
		 /* Local scan of the domain around a firenode. */
		 /* First a bounding box for scanning is computed.
		  * Then a local optimized polygon is defined locally
		  * to optimize the "point in polygon" algorithm,
		  * and lastly the scan is made */
 
		 // I/ Compute the region to be scanned
		 FFPoint swc, nec;
		 if ( fn->getPrev() == 0 or fn->getNext() == 0 ){
			 debugOutput<<getDomainID()<<": WARNING, launching a "
				 <<"FireDomain::firenodeBurningScan for a firenode "
				 <<"with no previous or next"<<endl;
			 swc.setX(fn->getX() - 2.*getPerimeterResolution());
			 swc.setY(fn->getY() - 2.*getPerimeterResolution());
			 nec.setX(fn->getX() + 2.*getPerimeterResolution());
			 nec.setY(fn->getY() + 2.*getPerimeterResolution());
		 } else {
			 computeBoundingBox(fn->getPrev(), fn, fn->getNext(), swc, nec);
		 }
		 // II/ Construct the vertices of a local optimized polygon
		 vector<double> nodesx, nodesy;
		 constructLocalSurroundingPolygon(fn, swc, nec, nodesx, nodesy);
 
		 // III/ Scanning the region
		 size_t nvert = nodesx.size();
		 double* vertx = 0;
		 double* verty = 0;
		 if ( nvert > 0 ){
			 vertx = new double[nvert];
			 verty = new double[nvert];
			 for ( size_t i = 0; i < nvert; i++ ){
				 vertx[i] = nodesx[i];
				 verty[i] = nodesy[i];
			 }
		 } else {
			 size_t nvert = fn->getFront()->getNumFN();
			 vertx = new double[nvert];
			 verty = new double[nvert];
			 fn->getFront()->storeVertices(vertx, verty, nvert);
		 }
		 singlePolygonAreaBurningScan(swc, nec, fn->getTime()
									  , fn->getFront()->isExpanding(), nvert, vertx, verty);
		 delete [] vertx;
		 delete [] verty;
 
	 }
 
	 void FireDomain::constructLocalSurroundingPolygon(FireNode* fn,
													   FFPoint& swc, FFPoint& nec, vector<double>& px, vector<double>& py){
 
		 /* Constructing a local polygon around the desired bounding box
		  * and taking into account the position of the firenodes in the bounding box */