Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 12042)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 12043)
@@ -748,19 +748,13 @@
 				    ./matlab/io/MatlabMatrixToMatrix.cpp\
 				    ./matlab/io/MatlabVectorToVector.cpp\
-					./matlab/io/MatlabVectorToDoubleVector.cpp\
-				    ./matlab/io/MatlabMatrixToDoubleMatrix.cpp\
-					./matlab/Container/Options.cpp\
-					./matlab/objects/Bamg/BamgGeom.h\
-					./matlab/objects/Bamg/BamgGeom.cpp\
-					./matlab/objects/Bamg/BamgMesh.h\
-					./matlab/objects/Bamg/BamgMesh.cpp\
-					./matlab/objects/Bamg/BamgOpts.h\
-					./matlab/objects/Bamg/BamgOpts.cpp\
-					./matlab/io/MatlabMatrixToPetscMatrix.cpp\
-					./matlab/io/MatlabVectorToPetscVector.cpp\
-					./matlab/io/PetscMatrixToDoubleMatrix.cpp\
-					./matlab/io/PetscVectorToDoubleVector.cpp\
-					./matlab/io/MatlabMatrixToSeqMat.cpp\
-					./matlab/io/MatlabVectorToSeqVec.cpp
+					 ./matlab/io/MatlabVectorToDoubleVector.cpp\
+					 ./matlab/io/MatlabMatrixToDoubleMatrix.cpp\
+					 ./matlab/Container/Options.cpp\
+					 ./matlab/io/MatlabMatrixToPetscMatrix.cpp\
+					 ./matlab/io/MatlabVectorToPetscVector.cpp\
+					 ./matlab/io/PetscMatrixToDoubleMatrix.cpp\
+					 ./matlab/io/PetscVectorToDoubleVector.cpp\
+					 ./matlab/io/MatlabMatrixToSeqMat.cpp\
+					 ./matlab/io/MatlabVectorToSeqVec.cpp
 #}}}
 #Modules sources{{{1
Index: /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp
===================================================================
--- /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp	(revision 12042)
+++ /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp	(revision 12043)
@@ -14,4 +14,5 @@
 #include "./matlabio.h"
 
+/*Primitive data types*/
 /*FUNCTION FetchData(double** pmatrix,int* pM,int *pN,const mxArray* dataref){{{1*/
 void FetchData(double** pmatrix,int* pM,int *pN,const mxArray* dataref){
@@ -242,25 +243,4 @@
 }
 /*}}}*/
-/*FUNCTION FetchData(Matrix** pmatrix,const mxArray* dataref){{{1*/
-void FetchData(Matrix** pmatrix,const mxArray* dataref){
-	
-	Matrix* outmatrix=NULL;
-	int dummy=0;
-
-	if (mxIsClass(dataref,"double") ){
-			
-		/*Convert matlab matrix to matrix: */
-		outmatrix=MatlabMatrixToMatrix(dataref);
-
-	}
-	else{
-		/*This is an error: we don't have the correct input!: */
-		_error_("Input parameter of class %s not supported yet",mxGetClassName(dataref));
-	}
-
-	/*Assign output pointers:*/
-	*pmatrix=outmatrix;
-}
-/*}}}*/
 /*FUNCTION FetchData(double** pvector,int* pM,const mxArray* dataref){{{1*/
 void FetchData(double** pvector,int* pM,const mxArray* dataref){
@@ -384,28 +364,4 @@
 	*pvector=outvector;
 	if (pM)*pM=outvector_rows;
-}
-/*}}}*/
-/*FUNCTION FetchData(Vector** pvector,const mxArray* dataref){{{1*/
-void FetchData(Vector** pvector,const mxArray* dataref){
-	
-	Vector* vector=NULL;
-	int dummy;
-
-	if(mxIsEmpty(dataref)){
-		/*Nothing to pick up. Just initialize matrix pointer to NULL: */
-		vector=new Vector(0);
-	}
-	else if (mxIsClass(dataref,"double") ){
-
-		/*Convert matlab vector to petsc vector: */
-		vector=MatlabVectorToVector(dataref);
-	}
-	else{
-		/*This is an error: we don't have the correct input!: */
-		_error_("Input parameter of class %s not supported yet",mxGetClassName(dataref));
-	}
-
-	/*Assign output pointers:*/
-	*pvector=vector;
 }
 /*}}}*/
@@ -526,2 +482,132 @@
 }
 /*}}}*/
+
+/*ISSM objects*/
+/*FUNCTION FetchData(Matrix** pmatrix,const mxArray* dataref){{{1*/
+void FetchData(Matrix** pmatrix,const mxArray* dataref){
+
+	Matrix* outmatrix=NULL;
+	int dummy=0;
+
+	if (mxIsClass(dataref,"double") ){
+
+		/*Convert matlab matrix to matrix: */
+		outmatrix=MatlabMatrixToMatrix(dataref);
+
+	}
+	else{
+		/*This is an error: we don't have the correct input!: */
+		_error_("Input parameter of class %s not supported yet",mxGetClassName(dataref));
+	}
+
+	/*Assign output pointers:*/
+	*pmatrix=outmatrix;
+}
+/*}}}*/
+/*FUNCTION FetchData(Vector** pvector,const mxArray* dataref){{{1*/
+void FetchData(Vector** pvector,const mxArray* dataref){
+
+	Vector* vector=NULL;
+	int dummy;
+
+	if(mxIsEmpty(dataref)){
+		/*Nothing to pick up. Just initialize matrix pointer to NULL: */
+		vector=new Vector(0);
+	}
+	else if (mxIsClass(dataref,"double") ){
+
+		/*Convert matlab vector to petsc vector: */
+		vector=MatlabVectorToVector(dataref);
+	}
+	else{
+		/*This is an error: we don't have the correct input!: */
+		_error_("Input parameter of class %s not supported yet",mxGetClassName(dataref));
+	}
+
+	/*Assign output pointers:*/
+	*pvector=vector;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgGeom** pbamggeom,const mxArray* dataref){{{1*/
+void FetchData(BamgGeom** pbamggeom,const mxArray* dataref){
+
+	/*Initialize output*/
+	BamgGeom* bamggeom=new BamgGeom();
+
+	/*Fetch all fields*/
+	FetchData(&bamggeom->Vertices,&bamggeom->VerticesSize[0],&bamggeom->VerticesSize[1],mxGetAssignedField(dataref,0,"Vertices"));
+	FetchData(&bamggeom->Edges, &bamggeom->EdgesSize[0], &bamggeom->EdgesSize[1], mxGetAssignedField(dataref,0,"Edges"));
+	FetchData(&bamggeom->Corners, &bamggeom->CornersSize[0], &bamggeom->CornersSize[1], mxGetAssignedField(dataref,0,"Corners"));
+	FetchData(&bamggeom->RequiredVertices,&bamggeom->RequiredVerticesSize[0],&bamggeom->RequiredVerticesSize[1],mxGetAssignedField(dataref,0,"RequiredVertices"));
+	FetchData(&bamggeom->RequiredEdges, &bamggeom->RequiredEdgesSize[0], &bamggeom->RequiredEdgesSize[1], mxGetAssignedField(dataref,0,"RequiredEdges"));
+	FetchData(&bamggeom->CrackedEdges,&bamggeom->CrackedEdgesSize[0],&bamggeom->CrackedEdgesSize[1],mxGetAssignedField(dataref,0,"CrackedEdges"));
+	FetchData(&bamggeom->SubDomains,&bamggeom->SubDomainsSize[0],&bamggeom->SubDomainsSize[1],mxGetAssignedField(dataref,0,"SubDomains"));
+
+	/*Assign output pointers:*/
+	*pbamggeom=bamggeom;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgMesh** pbamgmesh,const mxArray* dataref){{{1*/
+void FetchData(BamgMesh** pbamgmesh,const mxArray* dataref){
+
+	/*Initialize output*/
+	BamgMesh* bamgmesh=new BamgMesh();
+
+	/*Fetch all fields*/
+	FetchData(&bamgmesh->Vertices,&bamgmesh->VerticesSize[0],&bamgmesh->VerticesSize[1],mxGetAssignedField(dataref,0,"Vertices"));
+	FetchData(&bamgmesh->Edges, &bamgmesh->EdgesSize[0], &bamgmesh->EdgesSize[1], mxGetAssignedField(dataref,0,"Edges"));
+	FetchData(&bamgmesh->Triangles, &bamgmesh->TrianglesSize[0], &bamgmesh->TrianglesSize[1], mxGetAssignedField(dataref,0,"Triangles"));
+	FetchData(&bamgmesh->CrackedEdges,&bamgmesh->CrackedEdgesSize[0],&bamgmesh->CrackedEdgesSize[1],mxGetAssignedField(dataref,0,"CrackedEdges"));
+	FetchData(&bamgmesh->VerticesOnGeomEdge,&bamgmesh->VerticesOnGeomEdgeSize[0],&bamgmesh->VerticesOnGeomEdgeSize[1],mxGetAssignedField(dataref,0,"VerticesOnGeomEdge"));
+	FetchData(&bamgmesh->VerticesOnGeomVertex,&bamgmesh->VerticesOnGeomVertexSize[0],&bamgmesh->VerticesOnGeomVertexSize[1],mxGetAssignedField(dataref,0,"VerticesOnGeomVertex"));
+	FetchData(&bamgmesh->EdgesOnGeomEdge, &bamgmesh->EdgesOnGeomEdgeSize[0], &bamgmesh->EdgesOnGeomEdgeSize[1], mxGetAssignedField(dataref,0,"EdgesOnGeomEdge"));
+	FetchData(&bamgmesh->IssmSegments,&bamgmesh->IssmSegmentsSize[0],&bamgmesh->IssmSegmentsSize[1],mxGetAssignedField(dataref,0,"IssmSegments"));
+
+	/*Assign output pointers:*/
+	*pbamgmesh=bamgmesh;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgOpts** pbamgopts,const mxArray* dataref){{{1*/
+void FetchData(BamgOpts** pbamgopts,const mxArray* dataref){
+
+	/*Initialize output*/
+	BamgOpts* bamgopts=new BamgOpts();
+
+	/*Fetch all fields*/
+	FetchData(&bamgopts->anisomax,mxGetField(dataref,0,"anisomax"));
+	FetchData(&bamgopts->cutoff,mxGetField(dataref,0,"cutoff"));
+	FetchData(&bamgopts->coeff,mxGetField(dataref,0,"coeff"));
+	FetchData(&bamgopts->errg,mxGetField(dataref,0,"errg"));
+	FetchData(&bamgopts->gradation,mxGetField(dataref,0,"gradation"));
+	FetchData(&bamgopts->Hessiantype,mxGetField(dataref,0,"Hessiantype"));
+	FetchData(&bamgopts->MaxCornerAngle,mxGetField(dataref,0,"MaxCornerAngle"));
+	FetchData(&bamgopts->maxnbv,mxGetField(dataref,0,"maxnbv"));
+	FetchData(&bamgopts->maxsubdiv,mxGetField(dataref,0,"maxsubdiv"));
+	FetchData(&bamgopts->Metrictype,mxGetField(dataref,0,"Metrictype"));
+	FetchData(&bamgopts->nbjacobi,mxGetField(dataref,0,"nbjacobi"));
+	FetchData(&bamgopts->nbsmooth,mxGetField(dataref,0,"nbsmooth"));
+	FetchData(&bamgopts->omega,mxGetField(dataref,0,"omega"));
+	FetchData(&bamgopts->power,mxGetField(dataref,0,"power"));
+	FetchData(&bamgopts->verbose,mxGetField(dataref,0,"verbose"));
+
+	FetchData(&bamgopts->Crack,mxGetField(dataref,0,"Crack"));
+	FetchData(&bamgopts->geometricalmetric,mxGetField(dataref,0,"geometricalmetric"));
+	FetchData(&bamgopts->KeepVertices,mxGetField(dataref,0,"KeepVertices"));
+	FetchData(&bamgopts->splitcorners,mxGetField(dataref,0,"splitcorners"));
+
+	FetchData(&bamgopts->hmin,mxGetField(dataref,0,"hmin"));
+	FetchData(&bamgopts->hmax,mxGetField(dataref,0,"hmax"));
+	FetchData(&bamgopts->hminVertices,&bamgopts->hminVerticesSize[0],&bamgopts->hminVerticesSize[1],mxGetField(dataref,0,"hminVertices"));
+	FetchData(&bamgopts->hmaxVertices,&bamgopts->hmaxVerticesSize[0],&bamgopts->hmaxVerticesSize[1],mxGetField(dataref,0,"hmaxVertices"));
+	FetchData(&bamgopts->hVertices,&bamgopts->hVerticesSize[0],&bamgopts->hVerticesSize[1],mxGetField(dataref,0,"hVertices"));
+	FetchData(&bamgopts->metric,&bamgopts->metricSize[0],&bamgopts->metricSize[1],mxGetField(dataref,0,"metric"));
+	FetchData(&bamgopts->field,&bamgopts->fieldSize[0],&bamgopts->fieldSize[1],mxGetField(dataref,0,"field"));
+	FetchData(&bamgopts->err,&bamgopts->errSize[0],&bamgopts->errSize[1],mxGetField(dataref,0,"err"));
+
+	/*Additional checks*/
+	bamgopts->Check();
+
+	/*Assign output pointers:*/
+	*pbamgopts=bamgopts;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp
===================================================================
--- /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp	(revision 12042)
+++ /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp	(revision 12043)
@@ -15,12 +15,215 @@
 #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;
-	mxArray* dataref=NULL;
-	double*  matrix_ptr=NULL;
+	int      i,j;
 	int      rows,cols;
-	double*  tmatrix_ptr=NULL;
+	mxArray *dataref     = NULL;
+	double  *matrix_ptr  = NULL;
+	double  *tmatrix_ptr = NULL;
 	
 	if(matrix){
@@ -93,103 +296,30 @@
 }
 /*}}}*/
-/*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];
+
+/*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];
 			}
 		}
-		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);
-}
-/*}}}*/
+	}
+
+	/*Set matlab field*/
+	pfield=mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pfield,fieldcols);
+	mxSetN(pfield,fieldrows);
+	mxSetPr(pfield,fieldcopy);
+	mxSetField(dataref,0,fieldname,pfield);
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/matlab/io/matlabio.h
===================================================================
--- /issm/trunk-jpl/src/c/matlab/io/matlabio.h	(revision 12042)
+++ /issm/trunk-jpl/src/c/matlab/io/matlabio.h	(revision 12043)
@@ -27,4 +27,6 @@
 void WriteData(mxArray** pdataref,double scalar);
 void WriteData(mxArray** pdataref,char* string);
+void WriteData(mxArray** pdataref,BamgGeom* bamggeom);
+void WriteData(mxArray** pdataref,BamgMesh* bamgmesh);
 
 void FetchData(double** pmatrix,int* pM,int *pN,const mxArray* dataref);
@@ -44,4 +46,7 @@
 void FetchData(int* pinteger,const mxArray* dataref);
 void FetchData(bool* pbool,const mxArray* dataref);
+void FetchData(BamgGeom** bamggeom,const mxArray* dataref);
+void FetchData(BamgMesh** bamgmesh,const mxArray* dataref);
+void FetchData(BamgOpts** bamgopts,const mxArray* dataref);
 
 Option* OptionParse(char* name, const mxArray* prhs[]);
@@ -53,4 +58,5 @@
 
 mxArray* mxGetAssignedField(const mxArray* pmxa_array,int number, const char* field);
+void SetStructureField(mxArray* dataref,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer);
 int CheckNumMatlabArguments(int nlhs,int NLHS, int nrhs,int NRHS, const char* THISFUNCTION, void (*function)( void ));
 
Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgGeom.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/BamgGeom.h	(revision 12042)
+++ /issm/trunk-jpl/src/c/objects/Bamg/BamgGeom.h	(revision 12043)
@@ -26,9 +26,5 @@
 
 		BamgGeom();
-		BamgGeom(void* module_struct);
 		~BamgGeom();
-
-		void SetStructureFields(void* matlab_struct);
-		void SetStructureField(void* matlab_struct,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer);
 };
 
Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.cpp	(revision 12042)
+++ /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.cpp	(revision 12043)
@@ -18,7 +18,7 @@
 	this->CrackedEdgesSize[0]=0,              this->CrackedEdgesSize[1]=0;             this->CrackedEdges=NULL;
 
-	this->VerticesOnGeomVertexSize[0]=0, this->VerticesOnGeomVertexSize[1]=0;this->VerticesOnGeomVertex=NULL;
-	this->VerticesOnGeomEdgeSize[0]=0,   this->VerticesOnGeomEdgeSize[1]=0;  this->VerticesOnGeomEdge=NULL;
-	this->EdgesOnGeomEdgeSize[0]=0,      this->EdgesOnGeomEdgeSize[1]=0;     this->EdgesOnGeomEdge=NULL;
+	this->VerticesOnGeomVertexSize[0]=0,      this->VerticesOnGeomVertexSize[1]=0;     this->VerticesOnGeomVertex=NULL;
+	this->VerticesOnGeomEdgeSize[0]=0,        this->VerticesOnGeomEdgeSize[1]=0;       this->VerticesOnGeomEdge=NULL;
+	this->EdgesOnGeomEdgeSize[0]=0,           this->EdgesOnGeomEdgeSize[1]=0;          this->EdgesOnGeomEdge=NULL;
 
 	this->IssmEdgesSize[0]=0,                 this->IssmEdgesSize[1]=0;                this->IssmEdges=NULL;
@@ -28,6 +28,4 @@
 	this->NodalConnectivitySize[0]=0,         this->NodalConnectivitySize[1]=0;        this->NodalConnectivity=NULL;
 	this->NodalElementConnectivitySize[0]=0,  this->NodalElementConnectivitySize[1]=0; this->NodalElementConnectivity=NULL;
-
-
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.h	(revision 12042)
+++ /issm/trunk-jpl/src/c/objects/Bamg/BamgMesh.h	(revision 12043)
@@ -47,10 +47,5 @@
 
 		BamgMesh();
-		BamgMesh(void* module_struct);
 		~BamgMesh();
-
-		void SetStructureFields(void* module_struct);
-		void SetStructureField(void* module_struct,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer);
-
 };
 
Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgOpts.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/BamgOpts.h	(revision 12042)
+++ /issm/trunk-jpl/src/c/objects/Bamg/BamgOpts.h	(revision 12043)
@@ -50,5 +50,4 @@
 
 		BamgOpts();
-		BamgOpts(void* module_struct);
 		~BamgOpts();