/*!\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

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

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

	int i,n;
	Result* result=NULL;
	Result* newresult=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     sx_g=NULL;
	double* sx_g_serial=NULL;
	double* slopex=NULL;

	Vec     sy_g=NULL;
	double* sy_g_serial=NULL;
	double* slopey=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;

	Vec     grad_g=NULL;
	double* grad_g_serial=NULL;
	double* gradient=NULL;

	Vec     sigma_zz=NULL;
	double* sigma_zz_serial=NULL;

	Vec     v_g=NULL;
	double* v_g_serial=NULL;

	int numberofnodes,numberofvertices,numberofelements;

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

	/*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);
	if(analysis_type==PrognosticAnalysisEnum){
		fem_p=model->GetFormulation(PrognosticAnalysisEnum);
	}
	if(analysis_type==Prognostic2AnalysisEnum){
		fem_p=model->GetFormulation(Prognostic2AnalysisEnum);
	}
	if(analysis_type==TransientAnalysisEnum){
		fem_p=model->GetFormulation(PrognosticAnalysisEnum);
	}
	if(analysis_type==BalancedthicknessAnalysisEnum){
		fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
	}
	if(analysis_type==Balancedthickness2AnalysisEnum){
		fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
	}
	if(analysis_type==BalancedvelocitiesAnalysisEnum){
		fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum);
	}
	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(&numberofnodes,NumberOfNodesEnum);
					VecToMPISerial(&partition,fem_dh->partition->vector);
					fem_dh->parameters->FindParam(&yts,YtsEnum);
				}
				else{
					fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
					VecToMPISerial(&partition,fem_dhu->partition->vector);
					fem_dhu->parameters->FindParam(&yts,YtsEnum);
				}
				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(&numberofnodes,NumberOfNodesEnum);
						VecToMPISerial(&partition,fem_dh->partition->vector);
						fem_dh->parameters->FindParam(&yts,YtsEnum);
					}
					else{
						fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
						VecToMPISerial(&partition,fem_dhu->partition->vector);
						fem_dhu->parameters->FindParam(&yts,YtsEnum);
					}
					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(&numberofnodes,NumberOfNodesEnum);
					VecToMPISerial(&partition,fem_ds->partition->vector);
					fem_ds->parameters->FindParam(&yts,YtsEnum);
					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);
			xfree((void**)&vx);
			xfree((void**)&vy);
			xfree((void**)&vz);
			xfree((void**)&vel);
			VecFree(&u_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
					VecToMPISerial(&partition,fem_dh->partition->vector);
				}
				else{
					fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
					VecToMPISerial(&partition,fem_dhu->partition->vector);
				}
			}
			else{
				fem_ds->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
				VecToMPISerial(&partition,fem_ds->partition->vector);
			}

			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);
			xfree((void**)&pressure);
			VecFree(&p_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_t->partition->vector);

			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);
			xfree((void**)&temperature);
			VecFree(&t_g);
		}
		else if(strcmp(result->GetFieldName(),"grad_g")==0){
			
			/*easy, grad_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&grad_g);
			VecToMPISerial(&grad_g_serial,grad_g);
			fem_c->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_c->partition->vector);

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

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

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

			/*do some cleanup: */
			xfree((void**)&grad_g_serial);
			xfree((void**)&partition);
			xfree((void**)&gradient);
			VecFree(&grad_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
			fem_t->parameters->FindParam(&yts,YtsEnum);
			VecToMPISerial(&partition,fem_t->partition->vector);

			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);
			xfree((void**)&melting);
			VecFree(&m_g);
		}
		else if(strcmp(result->GetFieldName(),"h_g")==0){
			/*easy, h_g is of size numberofvertices, on 1 dof, just repartition: */
			result->GetField(&h_g);
			VecToMPISerial(&h_g_serial,h_g);
			fem_p->parameters->FindParam(&numberofvertices,NumberOfVerticesEnum);
			VecToMPISerial(&partition,fem_p->partition->vector);

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

			for(i=0;i<numberofvertices;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,numberofvertices);
			newresults->AddObject(newresult);

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

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

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

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

			/*do some cleanup: */
			xfree((void**)&v_g_serial);
			xfree((void**)&vel);
			xfree((void**)&partition);
			VecFree(&v_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_p->partition->vector);

			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);
			xfree((void**)&surface);
			VecFree(&s_g);
		}
		else if(strcmp(result->GetFieldName(),"sx_g")==0){
			/*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&sx_g);
			VecToMPISerial(&sx_g_serial,sx_g);
			fem_sl->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_sl->partition->vector);

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

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

			/*do some cleanup: */
			xfree((void**)&sx_g_serial);
			xfree((void**)&partition);
			xfree((void**)&slopex);
			VecFree(&sx_g);
		}
		else if(strcmp(result->GetFieldName(),"sy_g")==0){
			/*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
			result->GetField(&sy_g);
			VecToMPISerial(&sy_g_serial,sy_g);
			fem_sl->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_sl->partition->vector);

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

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

			/*do some cleanup: */
			xfree((void**)&sy_g_serial);
			xfree((void**)&partition);
			xfree((void**)&slopey);
			VecFree(&sy_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_p->partition->vector);

			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);
			xfree((void**)&bed);
			VecFree(&b_g);
		}
		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(&numberofnodes,NumberOfNodesEnum);
			VecToMPISerial(&partition,fem_dh->partition->vector);

			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);
			xfree((void**)&param_g);
			xfree((void**)&parameter);
		}
		else if(strcmp(result->GetFieldName(),"riftproperties")==0){
			result->GetField(&riftproperties);
			fem_dh->parameters->FindParam(&numrifts,NumRiftsEnum);
			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);
			xfree((void**)&riftproperties);

		}
		else if(strcmp(result->GetFieldName(),"sigma_zz")==0){
			/*easy, param_g is of size numberofelements, on 1 dof, just repartition: */
			fem_ds->parameters->FindParam(&numberofelements,NumberOfElementsEnum);
			result->GetField(&sigma_zz);
			VecToMPISerial(&sigma_zz_serial,sigma_zz);

			/*Ok, add parameter to newresults: */
			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"sigma_zz",sigma_zz_serial,numberofelements);
			newresults->AddObject(newresult);

			/*do some cleanup: */
			xfree((void**)&sigma_zz_serial);

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

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