/* \file WriteData.c:
 * \brief: general interface for writing data
 */

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

#include "../../include/include.h"
#include "../../shared/shared.h"
#include "./matlabio.h"

#include <mex.h>


/*Primitive data types*/
/*FUNCTION WriteData(mxArray** pdataref,double* matrix, int M,int N){{{1*/
void WriteData(mxArray** pdataref,double* matrix, int M,int N){

	mxArray *dataref  = NULL;
	double  *tmatrix  = NULL;
		
	if(matrix){
		/*create the matlab matrixwith Matlab's memory manager */   
		tmatrix=(double*)mxMalloc(M*N*sizeof(double));
		for(int i=0;i<M;i++){
			for(int j=0;j<N;j++){
				tmatrix[i*N+j]=matrix[j*M+i];
			}
		}
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(dataref,(mwSize)M);
		mxSetN(dataref,(mwSize)N);
		mxSetPr(dataref,(double*)tmatrix);
	}
	else{
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
	}
	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,int* matrix, int M,int N){{{1*/
void WriteData(mxArray** pdataref,int* matrix, int M,int N){

	mxArray* dataref = NULL;
	double*  tmatrix = NULL;

	if(matrix){

		/*convert to double matrix using Matlab's memory manager*/
		double* tmatrix=(double*)mxMalloc(M*N*sizeof(double));
		for(int i=0;i<M;i++){
			for(int j=0;j<N;j++){
				tmatrix[i*N+j]=(double)matrix[j*M+i];
			}
		}
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(dataref,(mwSize)M);
		mxSetN(dataref,(mwSize)N);
		mxSetPr(dataref,(double*)tmatrix);

	}
	else{
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
	}
	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,double* vector, int M){{{1*/
void WriteData(mxArray** pdataref,double* vector, int M){
	
	mxArray* dataref       = NULL;
	double*  vector_matlab = NULL;

	if(vector){

		/*create the matlab vector with Matlab's memory manager */
		vector_matlab=(double*)mxMalloc(M*sizeof(double));
		for(int i=0;i<M;i++) vector_matlab[i]=vector[i];
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(dataref,(mwSize)M);
		mxSetN(dataref,(mwSize)1);
		mxSetPr(dataref,vector_matlab);
	}
	else{
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
	}

	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,double scalar){{{1*/
void WriteData(mxArray** pdataref,double scalar){

	*pdataref=mxCreateDoubleScalar(scalar);
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,int integer){{{1*/
void WriteData(mxArray** pdataref,int integer){

		*pdataref=mxCreateDoubleScalar((double)integer);

}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,int boolean){{{1*/
void WriteData(mxArray** pdataref,bool boolean){

	*pdataref=mxCreateDoubleScalar((double)boolean);

}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,char* string){{{1*/
void WriteData(mxArray** pdataref,char* string){

		*pdataref=mxCreateString(string);
}
/*}}}*/

/*ISSM objects*/
/*FUNCTION WriteData(mxArray** pdataref,BamgGeom* bamggeom){{{1*/
void WriteData(mxArray** pdataref,BamgGeom* bamggeom){

	/*Intermediary*/
	int         i;
	mxArray    *dataref           = NULL;
	const int   numfields         = 7;
	const char *fnames[numfields];
	mwSize      ndim              = 2;
	mwSize      dimensions[2]     = {1,1};

	/*Initialize field names*/
	i=0;
	fnames[i++] = "Vertices";
	fnames[i++] = "Edges";
	fnames[i++] = "TangentAtEdges";
	fnames[i++] = "RequiredVertices";
	fnames[i++] = "RequiredEdges";
	fnames[i++] = "CrackedEdges";
	fnames[i++] = "SubDomains";
	_assert_(i==numfields);

	/*Initialize Matlab structure*/
	dataref=mxCreateStructArray(ndim,dimensions,numfields,fnames);

	/*set each matlab each field*/
	i=0;
	i++; SetStructureField(dataref,"Vertices",        bamggeom->VerticesSize[0],        bamggeom->VerticesSize[1],        bamggeom->Vertices);
	i++; SetStructureField(dataref,"Edges",           bamggeom->EdgesSize[0],           bamggeom->EdgesSize[1],           bamggeom->Edges);
	i++; SetStructureField(dataref,"TangentAtEdges",  bamggeom->TangentAtEdgesSize[0],  bamggeom->TangentAtEdgesSize[1],  bamggeom->TangentAtEdges);
	i++; SetStructureField(dataref,"RequiredVertices",bamggeom->RequiredVerticesSize[0],bamggeom->RequiredVerticesSize[1],bamggeom->RequiredVertices);
	i++; SetStructureField(dataref,"RequiredEdges",   bamggeom->RequiredEdgesSize[0],   bamggeom->RequiredEdgesSize[1],   bamggeom->RequiredEdges);
	i++; SetStructureField(dataref,"CrackedEdges",    bamggeom->CrackedEdgesSize[0],    bamggeom->CrackedEdgesSize[1],    bamggeom->CrackedEdges);
	i++; SetStructureField(dataref,"SubDomains",      bamggeom->SubDomainsSize[0],      bamggeom->SubDomainsSize[1],      bamggeom->SubDomains);
	_assert_(i==numfields);

	/*Assign output*/
	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,BamgMesh* bamgmesh){{{1*/
void WriteData(mxArray** pdataref,BamgMesh* bamgmesh){

	/*Intermediary*/
	int         i;
	mxArray    *dataref           = NULL;
	const int   numfields         = 16;
	const char *fnames[numfields];
	mwSize      ndim              = 2;
	mwSize      dimensions[2]     = {1,1};

	/*Initialize field names*/
	i=0;
	fnames[i++] = "Triangles";
	fnames[i++] = "Vertices";
	fnames[i++] = "Edges";
	fnames[i++] = "IssmSegments";
	fnames[i++] = "IssmEdges";
	fnames[i++] = "Quadrilaterals";
	fnames[i++] = "VerticesOnGeomVertex";
	fnames[i++] = "VerticesOnGeomEdge";
	fnames[i++] = "EdgesOnGeomEdge";
	fnames[i++] = "SubDomains";
	fnames[i++] = "SubDomainsFromGeom";
	fnames[i++] = "ElementConnectivity";
	fnames[i++] = "NodalConnectivity";
	fnames[i++] = "NodalElementConnectivity";
	fnames[i++] = "CrackedVertices";
	fnames[i++] = "CrackedEdges";
	_assert_(i==numfields);

	/*Initialize Matlab structure*/
	dataref=mxCreateStructArray(ndim,dimensions,numfields,fnames);

	/*set each matlab each field*/
	i=0;
	i++; SetStructureField(dataref,"Triangles", bamgmesh->TrianglesSize[0],bamgmesh->TrianglesSize[1], bamgmesh->Triangles);
	i++; SetStructureField(dataref,"Vertices",bamgmesh->VerticesSize[0], bamgmesh->VerticesSize[1],bamgmesh->Vertices);
	i++; SetStructureField(dataref,"Edges", bamgmesh->EdgesSize[0],bamgmesh->EdgesSize[1], bamgmesh->Edges);
	i++; SetStructureField(dataref,"IssmSegments",bamgmesh->IssmSegmentsSize[0], bamgmesh->IssmSegmentsSize[1],bamgmesh->IssmSegments);
	i++; SetStructureField(dataref,"IssmEdges", bamgmesh->IssmEdgesSize[0],bamgmesh->IssmEdgesSize[1], bamgmesh->IssmEdges);
	i++; SetStructureField(dataref,"Quadrilaterals",bamgmesh->QuadrilateralsSize[0], bamgmesh->QuadrilateralsSize[1],bamgmesh->Quadrilaterals);
	i++; SetStructureField(dataref,"VerticesOnGeomVertex",bamgmesh->VerticesOnGeomVertexSize[0],bamgmesh->VerticesOnGeomVertexSize[1], bamgmesh->VerticesOnGeomVertex);
	i++; SetStructureField(dataref,"VerticesOnGeomEdge",bamgmesh->VerticesOnGeomEdgeSize[0],bamgmesh->VerticesOnGeomEdgeSize[1], bamgmesh->VerticesOnGeomEdge);
	i++; SetStructureField(dataref,"EdgesOnGeomEdge", bamgmesh->EdgesOnGeomEdgeSize[0], bamgmesh->EdgesOnGeomEdgeSize[1],bamgmesh->EdgesOnGeomEdge);
	i++; SetStructureField(dataref,"SubDomains",bamgmesh->SubDomainsSize[0], bamgmesh->SubDomainsSize[1],bamgmesh->SubDomains);
	i++; SetStructureField(dataref,"SubDomainsFromGeom", bamgmesh->SubDomainsFromGeomSize[0], bamgmesh->SubDomainsFromGeomSize[1],bamgmesh->SubDomainsFromGeom);
	i++; SetStructureField(dataref,"ElementConnectivity",bamgmesh->ElementConnectivitySize[0],bamgmesh->ElementConnectivitySize[1], bamgmesh->ElementConnectivity);
	i++; SetStructureField(dataref,"NodalConnectivity",bamgmesh->NodalConnectivitySize[0],bamgmesh->NodalConnectivitySize[1], bamgmesh->NodalConnectivity);
	i++; SetStructureField(dataref,"NodalElementConnectivity", bamgmesh->NodalElementConnectivitySize[0], bamgmesh->NodalElementConnectivitySize[1],bamgmesh->NodalElementConnectivity);
	i++; SetStructureField(dataref,"CrackedVertices", bamgmesh->CrackedVerticesSize[0],bamgmesh->CrackedVerticesSize[1], bamgmesh->CrackedVertices);
	i++; SetStructureField(dataref,"CrackedEdges",bamgmesh->CrackedEdgesSize[0], bamgmesh->CrackedEdgesSize[1],bamgmesh->CrackedEdges);
	_assert_(i==numfields);

	/*Assign output*/
	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,Matrix* matrix){{{1*/
void WriteData(mxArray** pdataref,Matrix* matrix){
		
	int      i,j;
	int      rows,cols;
	mxArray *dataref     = NULL;
	double  *matrix_ptr  = NULL;
	double  *tmatrix_ptr = NULL;
	
	if(matrix){
		
		#ifdef _HAVE_PETSC_
		PetscMatrixToDoubleMatrix(&matrix_ptr,&rows,&cols,matrix->matrix);
		#else
		matrix_ptr=matrix->matrix->ToSerial();
		matrix->matrix->GetSize(&rows,cols);
		#endif

		/*Now transpose the matrix and allocate with Matlab's memory manager: */
		tmatrix_ptr=(double*)mxMalloc(rows*cols*sizeof(double));
		for(i=0;i<cols;i++){
			for(j=0;j<rows;j++){
				tmatrix_ptr[i*rows+j]=matrix_ptr[j*cols+i];
			}
		}
		
		/*create matlab matrix: */
		dataref=mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(dataref,rows); 
		mxSetN(dataref,cols);
		mxSetPr(dataref,tmatrix_ptr);

		/*Free ressources:*/
		xfree((void**)&matrix_ptr);

	}
	else{
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
	}

	*pdataref=dataref;
}
/*}}}*/
/*FUNCTION WriteData(mxArray** pdataref,Vector* vector){{{1*/
void WriteData(mxArray** pdataref,Vector* vector){
	
	mxArray* dataref=NULL;
	double*  vector_ptr=NULL;
	double*  vector_matlab=NULL;
	int      rows;
	
	if(vector){
		/*call toolkit routine: */
		#ifdef _HAVE_PETSC_
		PetscVectorToDoubleVector(&vector_ptr,&rows,vector->vector);
		#else
		vector_ptr=vector->vector->ToMPISerial();
		vector->vector->GetSize(&rows);
		#endif
		
		/*now create the matlab vector with Matlab's memory manager */
		vector_matlab=(double*)mxMalloc(rows*sizeof(double));
		for(int i=0;i<rows;i++) vector_matlab[i]=vector_ptr[i];

		dataref = mxCreateDoubleMatrix(0,0,mxREAL);                         
		mxSetM(dataref,rows);
		mxSetN(dataref,1);                                                                                          
		mxSetPr(dataref,vector_matlab);           
	}
	else{
		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
	}

	/*Clean-up and return*/
	xfree((void**)&vector_ptr);
	*pdataref=dataref;
}
/*}}}*/

/*Toolkit*/
/*FUNCTION SetStructureField{{{1*/
void SetStructureField(mxArray* dataref,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer){

	/*Intermediary*/
	int      i1,i2;
	mxArray *pfield  = NULL;

	/*Copy field using Matlab's API and transpose*/
	double*  fieldcopy=NULL;
	if (fieldrows*fieldcols){
		fieldcopy=(double*)mxMalloc(fieldrows*fieldcols*sizeof(double));
		for(i1=0;i1<fieldrows;i1++){
			for(i2=0;i2<fieldcols;i2++){
				fieldcopy[fieldcols*i1+i2]=fieldpointer[fieldrows*i2+i1];
			}
		}
	}

	/*Set matlab field*/
	pfield=mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pfield,fieldcols);
	mxSetN(pfield,fieldrows);
	mxSetPr(pfield,fieldcopy);
	mxSetField(dataref,0,fieldname,pfield);
}
/*}}}*/