Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 8900)
+++ /issm/trunk/src/c/Makefile.am	(revision 8901)
@@ -398,4 +398,6 @@
 					./toolkits/metis/metisincludes.h\
 					./toolkits/triangle/triangleincludes.h\
+					./toolkits/matlab/matlabincludes.h\
+					./toolkits/matlab/MatlabNArrayToNArray.cpp\
 					./toolkits.h\
 					./io/io.h\
@@ -407,5 +409,4 @@
 					./io/FetchParams.cpp\
 					./io/OptionParse.cpp\
-					./io/MatlabNArrayToNArray.cpp\
 					./io/pfopen.cpp\
 					./io/pfclose.cpp\
Index: sm/trunk/src/c/io/MatlabNArrayToNArray.cpp
===================================================================
--- /issm/trunk/src/c/io/MatlabNArrayToNArray.cpp	(revision 8900)
+++ 	(revision )
@@ -1,275 +1,0 @@
-/* \file MatlabNArrayToNArray.cpp
- * \brief: convert a sparse or dense matlab n-dimensional array to cpp n-dimensional array
- */
-
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./io.h"
-#include "../shared/shared.h"
-#include "../include/include.h"
-
-#ifdef _SERIAL_
-#include <mex.h>
-
-/*FUNCTION MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
-int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int  i,j,rows,cols;
-	int  numel,ndims;
-	int *size,*dims;
-	double* mxmatrix_ptr=NULL;
-	const mwSize* ipt=NULL;
-
-	/*output: */
-	double* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-	double* pr=NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*get Matlab matrix information: */
-	numel=mxGetNumberOfElements(mxmatrix);
-	ndims=mxGetNumberOfDimensions(mxmatrix);
-	ipt  =mxGetDimensions(mxmatrix);
-	size =(int *) xcalloc(ndims,sizeof(int));
-	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-		nnz=mxGetNzmax(mxmatrix);
-		nz=(int)((double)nnz/(double)rows);
-
-		matrix=(double*)xcalloc(rows*cols,sizeof(double));
-
-		/*Now, get ir,jc and pr: */
-		ir=mxGetIr(mxmatrix);
-		jc=mxGetJc(mxmatrix);
-		pr=mxGetPr(mxmatrix);
-
-		/*Now, start inserting data into double* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				*(matrix+rows*ir[count]+i)=pr[count];
-				count++;
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
-		
-		/*Create serial matrix: */
-		matrix=(double*)xcalloc(numel,sizeof(double));
-
-		dims=(int *) xcalloc(ndims,sizeof(int));
-		for(i=0;i<numel;i++){
-			ColumnWiseDimsFromIndex(dims,i,size,ndims);
-			j=IndexFromRowWiseDims(dims,size,ndims);
-			*(matrix+j)=*(mxmatrix_ptr+i);
-		}
-		xfree((void**)&dims);
-		
-	}
-
-	/*Assign output pointer: */
-	*pmatrix=matrix;
-	*pmatrix_numel=numel;
-	*pmatrix_ndims=ndims;
-	*pmatrix_size=size;
-
-	return 1;
-}
-/*}}}*/
-/*FUNCTION MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
-int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int  i,j,rows,cols;
-	int  numel,ndims;
-	int *size,*dims;
-	bool* mxmatrix_ptr=NULL;
-	const mwSize* ipt=NULL;
-
-	/*output: */
-	bool* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-	bool*   pm=NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*get Matlab matrix information: */
-	numel=mxGetNumberOfElements(mxmatrix);
-	ndims=mxGetNumberOfDimensions(mxmatrix);
-	ipt  =mxGetDimensions(mxmatrix);
-	size =(int *) xcalloc(ndims,sizeof(int));
-	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-		nnz=mxGetNzmax(mxmatrix);
-		nz=(int)((double)nnz/(double)rows);
-
-		matrix=(bool*)xcalloc(rows*cols,sizeof(bool));
-
-		/*Now, get ir,jc and pm: */
-		ir=mxGetIr(mxmatrix);
-		jc=mxGetJc(mxmatrix);
-		pm=(bool*)mxGetData(mxmatrix);
-
-		/*Now, start inserting data into bool* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				*(matrix+rows*ir[count]+i)=pm[count];
-				count++;
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=(bool*)mxGetData(mxmatrix);
-		
-		/*Create serial matrix: */
-		matrix=(bool*)xcalloc(numel,sizeof(bool));
-
-		dims=(int *) xcalloc(ndims,sizeof(int));
-		for(i=0;i<numel;i++){
-			ColumnWiseDimsFromIndex(dims,i,size,ndims);
-			j=IndexFromRowWiseDims(dims,size,ndims);
-			*(matrix+j)=(bool)*(mxmatrix_ptr+i);
-		}
-		xfree((void**)&dims);
-		
-	}
-
-	/*Assign output pointer: */
-	*pmatrix=matrix;
-	*pmatrix_numel=numel;
-	*pmatrix_ndims=ndims;
-	*pmatrix_size=size;
-
-	return 1;
-}
-/*}}}*/
-/*FUNCTION MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
-int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int  i,j,rows,cols;
-	int  numel,ndims;
-	int *size,*dims;
-	mxChar* mxmatrix_ptr=NULL;
-	const mwSize* ipt=NULL;
-
-	/*output: */
-	char* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-	char*   pm=NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*get Matlab matrix information: */
-	numel=mxGetNumberOfElements(mxmatrix);
-	ndims=mxGetNumberOfDimensions(mxmatrix);
-	ipt  =mxGetDimensions(mxmatrix);
-	size =(int *) xcalloc(ndims,sizeof(int));
-	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-		nnz=mxGetNzmax(mxmatrix);
-		nz=(int)((double)nnz/(double)rows);
-
-		matrix=(char*)xcalloc(rows*cols,sizeof(double));
-
-		/*Now, get ir,jc and pm: */
-		ir=mxGetIr(mxmatrix);
-		jc=mxGetJc(mxmatrix);
-		pm=(char*)mxGetData(mxmatrix);
-
-		/*Now, start inserting data into char* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				*(matrix+rows*ir[count]+i)=(char)pm[count];
-				count++;
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=mxGetChars(mxmatrix);
-		
-		/*Create serial matrix: */
-		matrix=(char*)xcalloc(numel+1,sizeof(mxChar));
-
-		/*looping code adapted from Matlab example explore.c: */
-		int elements_per_page = size[0] * size[1];
-		/* total_number_of_pages = size[2] x size[3] x ... x size[N-1] */
-		int total_number_of_pages = 1;
-		for (i=2; i<ndims; i++) {
-			total_number_of_pages *= size[i];
-		}
-
-		i=0;
-		for (int page=0; page < total_number_of_pages; page++) {
-			int row;
-			/* On each page, walk through each row. */
-			for (row=0; row<size[0]; row++)  {
-				int column;
-				j = (page * elements_per_page) + row;
-
-				/* Walk along each column in the current row. */
-				for (column=0; column<size[1]; column++) {
-					*(matrix+i++)=(char)*(mxmatrix_ptr+j);
-					j += size[0];
-				}
-			}
-		}
-
-	}
-
-	/*Assign output pointer: */
-	*pmatrix=matrix;
-	*pmatrix_numel=numel;
-	*pmatrix_ndims=ndims;
-	*pmatrix_size=size;
-
-	return 1;
-}
-/*}}}*/
-#endif
Index: /issm/trunk/src/c/io/io.h
===================================================================
--- /issm/trunk/src/c/io/io.h	(revision 8900)
+++ /issm/trunk/src/c/io/io.h	(revision 8901)
@@ -80,17 +80,9 @@
 OptionCell*     OptionCellParse( char* name, const mxArray* prhs[]);
 
-int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
-int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
-int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
-
 #endif
 
 #ifdef _PARALLEL_
-
-
 void WriteData(int* pdummy,void* data,char* data_type);
 void WriteDataToDisk(void* data,int* pM,int* pN,char* datatype,FILE* fid);
-
-
 #endif
 
Index: sm/trunk/src/c/objects/NodeSets.cpp
===================================================================
--- /issm/trunk/src/c/objects/NodeSets.cpp	(revision 8900)
+++ 	(revision )
@@ -1,86 +1,0 @@
-/*!\file: NodeSet.cpp
- * \brief routines for handling node sets (m,n,f and s sets)
- */ 
-
-#include "../include/include.h"
-#include "../shared/shared.h"
-#include "./objects.h"
-
-/*constructors and destructors*/
-/*FUNCTION NodeSets::NodeSets{{{1*/
-NodeSets::NodeSets( double* nodesets_pv_f,double* nodesets_pv_s, int nodesets_gsize,int nodesets_fsize,int nodesets_ssize){
-
-	gsize=nodesets_gsize;
-	fsize=nodesets_fsize;
-	ssize=nodesets_ssize;
-
-	pv_f=nodesets_pv_f;
-	pv_s=nodesets_pv_s;
-
-}/*}}}*/
-/*FUNCTION NodeSets::~NodeSets{{{1*/
-NodeSets::~NodeSets(void){
-
-	xfree((void**)&pv_f);
-	xfree((void**)&pv_s);
-
-}/*}}}*/
-
-/*Object management*/
-/*FUNCTION NodeSets::Echo{{{1*/
-void NodeSets::Echo(void){
-	
-	int i;
-
-	printf("\nNodeSets echo:\n");
-	printf(" gsize: %i\n",gsize);
-	printf(" fsize: %i\n",fsize);
-	printf(" ssize: %i\n",ssize);
-	printf(" pv_f: %p\n",pv_f);
-	printf(" pv_s: %p\n",pv_s);
-
-}/*}}}*/
-/*FUNCTION NodeSets::DeepEcho{{{1*/
-void NodeSets::DeepEcho(void){
-	
-	int i;
-
-	printf("\nNodeSets echo:\n");
-	printf(" gsize: %i\n",gsize);
-	printf(" fsize: %i\n",fsize);
-	printf(" ssize: %i\n",ssize);
-
-	if(pv_f){
-		printf("f set:\n");
-		for(i=0;i<fsize;i++){
-			printf("%g\n",pv_f[i]);
-		}
-	}
-	
-	if(pv_s){
-		printf("s set:\n");
-		for(i=0;i<ssize;i++){
-			printf("%g\n",pv_s[i]);
-		}
-	}	
-}/*}}}*/
-/*FUNCTION NodeSets::GetGSize{{{1*/
-int NodeSets::GetGSize(){
-	return gsize;
-}/*}}}*/
-/*FUNCTION NodeSets::GetFSize{{{1*/
-int NodeSets::GetFSize(){
-	return fsize;
-}/*}}}*/
-/*FUNCTION NodeSets::GetSSize{{{1*/
-int NodeSets::GetSSize(){
-	return ssize;
-}/*}}}*/
-/*FUNCTION NodeSets::GetPV_F{{{1*/
-double* NodeSets::GetPV_F(){
-	return pv_f;
-}/*}}}*/
-/*FUNCTION NodeSets::GetPV_S{{{1*/
-double* NodeSets::GetPV_S(){
-	return pv_s;
-}/*}}}*/
Index: sm/trunk/src/c/objects/NodeSets.h
===================================================================
--- /issm/trunk/src/c/objects/NodeSets.h	(revision 8900)
+++ 	(revision )
@@ -1,37 +1,0 @@
-/*!\file:  NodeSets.h
- * \brief header file for  a strcture holding all the nodesets.
- */ 
-
-#ifndef _NODESETS_H
-#define _NODESETS_H
-
-/*Headers:*/
-/*{{{1*/
-/*}}}*/
-
-class NodeSets {
-
-	private:
-		int gsize;
-		int fsize;
-		int ssize;
-
-		double* pv_f;
-		double* pv_s;
-
-	public:
-
-	NodeSets( double* pv_f,double* pv_s,int gsize,int fsize,int ssize);
-	~NodeSets();
-
-	void Echo();
-	void DeepEcho(void);
-	int GetGSize();
-	int GetFSize();
-	int GetSSize();
-	double* GetPV_F();
-	double* GetPV_S();
-
-};
-#endif  /* _NODESETS_H */
-
Index: /issm/trunk/src/c/toolkits/matlab/MatlabNArrayToNArray.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/matlab/MatlabNArrayToNArray.cpp	(revision 8901)
+++ /issm/trunk/src/c/toolkits/matlab/MatlabNArrayToNArray.cpp	(revision 8901)
@@ -0,0 +1,274 @@
+/* \file MatlabNArrayToNArray.cpp
+ * \brief: convert a sparse or dense matlab n-dimensional array to cpp n-dimensional array
+ */
+
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+
+#ifdef _SERIAL_
+#include <mex.h>
+
+/*FUNCTION MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
+int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int  i,j,rows,cols;
+	int  numel,ndims;
+	int *size,*dims;
+	double* mxmatrix_ptr=NULL;
+	const mwSize* ipt=NULL;
+
+	/*output: */
+	double* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex*    ir=NULL;
+	mwIndex*    jc=NULL;
+	double* pr=NULL;
+	int     count;
+	int     nnz;
+	int     nz;
+
+	/*get Matlab matrix information: */
+	numel=mxGetNumberOfElements(mxmatrix);
+	ndims=mxGetNumberOfDimensions(mxmatrix);
+	ipt  =mxGetDimensions(mxmatrix);
+	size =(int *) xcalloc(ndims,sizeof(int));
+	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+		nnz=mxGetNzmax(mxmatrix);
+		nz=(int)((double)nnz/(double)rows);
+
+		matrix=(double*)xcalloc(rows*cols,sizeof(double));
+
+		/*Now, get ir,jc and pr: */
+		ir=mxGetIr(mxmatrix);
+		jc=mxGetJc(mxmatrix);
+		pr=mxGetPr(mxmatrix);
+
+		/*Now, start inserting data into double* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				*(matrix+rows*ir[count]+i)=pr[count];
+				count++;
+			}
+		}
+
+	}
+	else{
+
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
+		
+		/*Create serial matrix: */
+		matrix=(double*)xcalloc(numel,sizeof(double));
+
+		dims=(int *) xcalloc(ndims,sizeof(int));
+		for(i=0;i<numel;i++){
+			ColumnWiseDimsFromIndex(dims,i,size,ndims);
+			j=IndexFromRowWiseDims(dims,size,ndims);
+			*(matrix+j)=*(mxmatrix_ptr+i);
+		}
+		xfree((void**)&dims);
+		
+	}
+
+	/*Assign output pointer: */
+	*pmatrix=matrix;
+	*pmatrix_numel=numel;
+	*pmatrix_ndims=ndims;
+	*pmatrix_size=size;
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
+int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int  i,j,rows,cols;
+	int  numel,ndims;
+	int *size,*dims;
+	bool* mxmatrix_ptr=NULL;
+	const mwSize* ipt=NULL;
+
+	/*output: */
+	bool* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex*    ir=NULL;
+	mwIndex*    jc=NULL;
+	bool*   pm=NULL;
+	int     count;
+	int     nnz;
+	int     nz;
+
+	/*get Matlab matrix information: */
+	numel=mxGetNumberOfElements(mxmatrix);
+	ndims=mxGetNumberOfDimensions(mxmatrix);
+	ipt  =mxGetDimensions(mxmatrix);
+	size =(int *) xcalloc(ndims,sizeof(int));
+	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+		nnz=mxGetNzmax(mxmatrix);
+		nz=(int)((double)nnz/(double)rows);
+
+		matrix=(bool*)xcalloc(rows*cols,sizeof(bool));
+
+		/*Now, get ir,jc and pm: */
+		ir=mxGetIr(mxmatrix);
+		jc=mxGetJc(mxmatrix);
+		pm=(bool*)mxGetData(mxmatrix);
+
+		/*Now, start inserting data into bool* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				*(matrix+rows*ir[count]+i)=pm[count];
+				count++;
+			}
+		}
+
+	}
+	else{
+
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=(bool*)mxGetData(mxmatrix);
+		
+		/*Create serial matrix: */
+		matrix=(bool*)xcalloc(numel,sizeof(bool));
+
+		dims=(int *) xcalloc(ndims,sizeof(int));
+		for(i=0;i<numel;i++){
+			ColumnWiseDimsFromIndex(dims,i,size,ndims);
+			j=IndexFromRowWiseDims(dims,size,ndims);
+			*(matrix+j)=(bool)*(mxmatrix_ptr+i);
+		}
+		xfree((void**)&dims);
+		
+	}
+
+	/*Assign output pointer: */
+	*pmatrix=matrix;
+	*pmatrix_numel=numel;
+	*pmatrix_ndims=ndims;
+	*pmatrix_size=size;
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{1*/
+int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int  i,j,rows,cols;
+	int  numel,ndims;
+	int *size,*dims;
+	mxChar* mxmatrix_ptr=NULL;
+	const mwSize* ipt=NULL;
+
+	/*output: */
+	char* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex*    ir=NULL;
+	mwIndex*    jc=NULL;
+	char*   pm=NULL;
+	int     count;
+	int     nnz;
+	int     nz;
+
+	/*get Matlab matrix information: */
+	numel=mxGetNumberOfElements(mxmatrix);
+	ndims=mxGetNumberOfDimensions(mxmatrix);
+	ipt  =mxGetDimensions(mxmatrix);
+	size =(int *) xcalloc(ndims,sizeof(int));
+	for (i=0; i<ndims; i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+		nnz=mxGetNzmax(mxmatrix);
+		nz=(int)((double)nnz/(double)rows);
+
+		matrix=(char*)xcalloc(rows*cols,sizeof(double));
+
+		/*Now, get ir,jc and pm: */
+		ir=mxGetIr(mxmatrix);
+		jc=mxGetJc(mxmatrix);
+		pm=(char*)mxGetData(mxmatrix);
+
+		/*Now, start inserting data into char* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				*(matrix+rows*ir[count]+i)=(char)pm[count];
+				count++;
+			}
+		}
+
+	}
+	else{
+
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=mxGetChars(mxmatrix);
+		
+		/*Create serial matrix: */
+		matrix=(char*)xcalloc(numel+1,sizeof(mxChar));
+
+		/*looping code adapted from Matlab example explore.c: */
+		int elements_per_page = size[0] * size[1];
+		/* total_number_of_pages = size[2] x size[3] x ... x size[N-1] */
+		int total_number_of_pages = 1;
+		for (i=2; i<ndims; i++) {
+			total_number_of_pages *= size[i];
+		}
+
+		i=0;
+		for (int page=0; page < total_number_of_pages; page++) {
+			int row;
+			/* On each page, walk through each row. */
+			for (row=0; row<size[0]; row++)  {
+				int column;
+				j = (page * elements_per_page) + row;
+
+				/* Walk along each column in the current row. */
+				for (column=0; column<size[1]; column++) {
+					*(matrix+i++)=(char)*(mxmatrix_ptr+j);
+					j += size[0];
+				}
+			}
+		}
+
+	}
+
+	/*Assign output pointer: */
+	*pmatrix=matrix;
+	*pmatrix_numel=numel;
+	*pmatrix_ndims=ndims;
+	*pmatrix_size=size;
+
+	return 1;
+}
+/*}}}*/
+#endif
Index: /issm/trunk/src/c/toolkits/matlab/matlabincludes.h
===================================================================
--- /issm/trunk/src/c/toolkits/matlab/matlabincludes.h	(revision 8901)
+++ /issm/trunk/src/c/toolkits/matlab/matlabincludes.h	(revision 8901)
@@ -0,0 +1,15 @@
+/* \file matlabincludes.h
+ * \brief all includes for Matlab patches
+ */
+
+#ifndef _MATLAB_INCLUDES_H_
+#define _MATLAB_INCLUDES_H_
+
+#ifdef _SERIAL_
+#include <mex.h>
+int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
+int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
+int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
+#endif
+
+#endif
Index: /issm/trunk/src/c/toolkits/toolkits.h
===================================================================
--- /issm/trunk/src/c/toolkits/toolkits.h	(revision 8900)
+++ /issm/trunk/src/c/toolkits/toolkits.h	(revision 8901)
@@ -10,4 +10,5 @@
 #include "./metis/metisincludes.h"
 #include "./triangle/triangleincludes.h"
+#include "./matlab/matlabincludes.h"
 
 #endif