/*!\file:  ProcessResults.cpp
 * \brief: go through results dataset, and for each result, process it for easier retrieval 
 * by the Matlab side. This usually means splitting the velocities from the g-size nodeset 
 * to the grid set (ug->vx,vy,vz), same for pressure (p_g->pressure), etc ... It also implies 
 * departitioning of the results.
 * This phase is necessary prior to outputting the results on disk.
 */ 

#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#undef __FUNCT__ 
#define __FUNCT__ "ProcessResults"

#include "../DataSet/DataSet.h"
#include "../objects/objects.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../shared/shared.h"

void ProcessResults(DataSet** presults,Model* model,int analysis_type){

	int i,n;
	Result* result=NULL;
	Result* newresult=NULL;
	
	/*input: */
	DataSet* results=NULL;

	/*output: */
	DataSet* newresults=NULL;

	/*fem diagnostic models: */
	FemModel* fem_dh=NULL;
	FemModel* fem_dv=NULL;
	FemModel* fem_dhu=NULL;
	FemModel* fem_ds=NULL;
	FemModel* fem_sl=NULL;

	/*fem thermal models: */
	FemModel* fem_t=NULL;

	/*fem prognostic models: */
	FemModel* fem_p=NULL;

	/*fem control models: */
	FemModel* fem_c=NULL;

	/*some parameters*/
	int ishutter;
	int ismacayealpattyn;
	int isstokes;
	int dim;

	/*intermediary: */
	Vec     u_g=NULL;
	double* u_g_serial=NULL;
	double* vx=NULL;
	double* vy=NULL;
	double* vz=NULL;
	double* vel=NULL;
	Vec     p_g=NULL;
	double* p_g_serial=NULL;
	double* pressure=NULL;
	double* partition=NULL;
	double  yts;

	double* param_g=NULL;
	double* parameter=NULL;
	Vec     riftproperties=NULL;
	double* riftproperties_serial=NULL;
	int     numrifts=0;

	Vec     s_g=NULL;
	double* s_g_serial=NULL;
	double* surface=NULL;

	Vec     b_g=NULL;
	double* b_g_serial=NULL;
	double* bed=NULL;

	Vec     h_g=NULL;
	double* h_g_serial=NULL;
	double* thickness=NULL;

	Vec     t_g=NULL;
	double* t_g_serial=NULL;
	double* temperature=NULL;

	Vec     m_g=NULL;
	double* m_g_serial=NULL;
	double* melting=NULL;

	int numberofnodes;

	/*recover input results: */
	results=*presults;

	/*Initialize new results: */
	newresults=new DataSet(ResultsEnum());
	
	/*some flags needed: */
	model->FindParam(&dim,"dim");
	model->FindParam(&ishutter,"ishutter");
	model->FindParam(&isstokes,"isstokes");
	model->FindParam(&ismacayealpattyn,"ismacayealpattyn");

	/*Recover femmodels first: */
	fem_dh=model->GetFormulation(DiagnosticAnalysisEnum(),HorizAnalysisEnum());
	fem_c=fem_dh;
	fem_dv=model->GetFormulation(DiagnosticAnalysisEnum(),VertAnalysisEnum());
	fem_ds=model->GetFormulation(DiagnosticAnalysisEnum(),StokesAnalysisEnum());
	fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum());
	fem_sl=model->GetFormulation(SlopeComputeAnalysisEnum());
	fem_p=model->GetFormulation(PrognosticAnalysisEnum());
	fem_t=model->GetFormulation(ThermalAnalysisEnum());

	for(n=0;n<results->Size();n++){
		result=(Result*)results->GetObjectByOffset(n);

		if(strcmp(result->GetFieldName(),"u_g")==0){

			/*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
			 *Stokes on 4 dofs: */
			result->GetField(&u_g);
			VecToMPISerial(&u_g_serial,u_g);

			//2d results -> 2 dofs per node
			if (dim==2){
				/*ok, 2 dofs, on number of nodes: */
				if(ismacayealpattyn){
					fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
					VecToMPISerial(&partition,fem_dh->partition);
					fem_dh->parameters->FindParam((void*)&yts,"yts");
				}
				else{
					fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
					VecToMPISerial(&partition,fem_dhu->partition);
					fem_dhu->parameters->FindParam((void*)&yts,"yts");
				}
				vx=(double*)xmalloc(numberofnodes*sizeof(double));
				vy=(double*)xmalloc(numberofnodes*sizeof(double));
				vz=(double*)xmalloc(numberofnodes*sizeof(double)); 
				vel=(double*)xmalloc(numberofnodes*sizeof(double)); 

				for(i=0;i<numberofnodes;i++){
					vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
					vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
					vz[i]=0;
					vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
				}
			}
			//3d results -> 3 or 4 (stokes) dofs per node
			else{
				if(!isstokes){
					/*ok, 3 dofs, on number of nodes: */
					if(ismacayealpattyn){
						fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
						VecToMPISerial(&partition,fem_dh->partition);
						fem_dh->parameters->FindParam((void*)&yts,"yts");
					}
					else{
						fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
						VecToMPISerial(&partition,fem_dhu->partition);
						fem_dhu->parameters->FindParam((void*)&yts,"yts");
					}
					vx=(double*)xmalloc(numberofnodes*sizeof(double));
					vy=(double*)xmalloc(numberofnodes*sizeof(double));
					vz=(double*)xmalloc(numberofnodes*sizeof(double)); 
					vel=(double*)xmalloc(numberofnodes*sizeof(double)); 

					for(i=0;i<numberofnodes;i++){
						vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
						vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
						vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
						vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
					}
				}
				else{
					/* 4 dofs on number of nodes. discard pressure: */
					fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
					VecToMPISerial(&partition,fem_ds->partition);
					fem_ds->parameters->FindParam((void*)&yts,"yts");
					vx=(double*)xmalloc(numberofnodes*sizeof(double));
					vy=(double*)xmalloc(numberofnodes*sizeof(double));
					vz=(double*)xmalloc(numberofnodes*sizeof(double));
					vel=(double*)xmalloc(numberofnodes*sizeof(double));
					for(i=0;i<numberofnodes;i++){
						vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
						vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
						vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
						vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
					}
				}
			}

			/*Ok, add vx,vy and vz to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
			newresults->AddObject(newresult);

			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
			newresults->AddObject(newresult);

			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
			newresults->AddObject(newresult);

			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&u_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"p_g")==0){
			/*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&p_g);
			VecToMPISerial(&p_g_serial,p_g);

			if(!isstokes){
				if(ismacayealpattyn){
					fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
					VecToMPISerial(&partition,fem_dh->partition);
				}
				else{
					fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
					VecToMPISerial(&partition,fem_dhu->partition);
				}
			}
			else{
				fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
				VecToMPISerial(&partition,fem_ds->partition);
			}

			pressure=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				pressure[i]=p_g_serial[(int)partition[i]];
			}
			
			/*Ok, add pressure,vy and vz to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&p_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"t_g")==0){
			/*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&t_g);
			VecToMPISerial(&t_g_serial,t_g);
			fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			VecToMPISerial(&partition,fem_t->partition);

			temperature=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				temperature[i]=t_g_serial[(int)partition[i]];
			}

			/*Ok, add pressure to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&t_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"m_g")==0){
			/*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&m_g);
			VecToMPISerial(&m_g_serial,m_g);
			fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			fem_t->parameters->FindParam((void*)&yts,"yts");
			VecToMPISerial(&partition,fem_t->partition);

			melting=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				melting[i]=m_g_serial[(int)partition[i]]*yts;
			}

			/*Ok, add pressure to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&m_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"h_g")==0){
			/*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&h_g);
			VecToMPISerial(&h_g_serial,h_g);
			fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			VecToMPISerial(&partition,fem_p->partition);

			thickness=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				thickness[i]=h_g_serial[(int)partition[i]];
			}
			
			/*Ok, add pressure to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&h_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"s_g")==0){
			/*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&s_g);
			VecToMPISerial(&s_g_serial,s_g);
			fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			VecToMPISerial(&partition,fem_p->partition);

			surface=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				surface[i]=s_g_serial[(int)partition[i]];
			}
			
			/*Ok, add pressure to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&s_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"b_g")==0){
			/*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&b_g);
			VecToMPISerial(&b_g_serial,b_g);
			fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			VecToMPISerial(&partition,fem_p->partition);

			bed=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				bed[i]=b_g_serial[(int)partition[i]];
			}
			
			/*Ok, add pressure to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&b_g_serial);
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"param_g")==0){
			/*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&param_g);
			fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
			VecToMPISerial(&partition,fem_dh->partition);

			parameter=(double*)xmalloc(numberofnodes*sizeof(double));

			for(i=0;i<numberofnodes;i++){
				parameter[i]=param_g[(int)partition[i]];
			}
			
			/*Ok, add parameter to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&partition);
		}
		else if(strcmp(result->GetFieldName(),"riftproperties")==0){
			result->GetField(&riftproperties);
			fem_dh->parameters->FindParam((void*)&numrifts,"numrifts");
			VecToMPISerial(&riftproperties_serial,riftproperties);
			
			/*Ok, add parameter to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts);
			newresults->AddObject(newresult);

		}
		else{
			/*Just copy the result into the new results dataset: */
			newresults->AddObject(result->copy());
		}
	}

	/*Delete results: */
	delete results;


	/*Assign output pointers:*/
	*presults=newresults;
}