/*
 * CreateElementsNodesAndMaterialsDiagnosticHoriz.c:
 */

#undef __FUNCT__ 
#define __FUNCT__ "CreateElementsNodesAndMaterialsDiagnosticHoriz"

#include "../DataSet/DataSet.h"
#include "../toolkits/toolkits.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../objects/objects.h"
#include "../shared/shared.h"
#include "./ModelProcessorx.h"
#include "./MeshPartitionx/MeshPartitionx.h"


int	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){


	/*output: int* epart, int* my_grids, double* my_bordergrids*/


	int noerr=1;
	int i,j,k,n;
	extern int my_rank;
	extern int num_procs;

	/*DataSets: */
	DataSet*    elements  = NULL;
	DataSet*    nodes = NULL;
	DataSet*    materials = NULL;
	
	/*Objects: */
	Node*       node   = NULL;
	Tria*       tria = NULL;
	Penta*      penta = NULL;
	Matice*     matice  = NULL;
	Matpar*     matpar  = NULL;

	int         analysis_type;
	
	/*output: */
	int* epart=NULL; //element partitioning.
	int* npart=NULL; //node partitioning.
	int* my_grids=NULL;
	double* my_bordergrids=NULL;


	/*intermediary: */
	int elements_width; //size of elements
	double B_avg;
			
	/*tria constructor input: */
	int tria_id;
	int tria_mid;
	int tria_mparid;
	int tria_g[3];
	double tria_h[3];
	double tria_s[3];
	double tria_b[3];
	double tria_k[3];
	int    tria_friction_type;
	double tria_p;
	double tria_q;
	int    tria_shelf;
	double tria_fit;
	double tria_meanvel;/*!scaling ratio for velocities*/
	double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
	int    tria_acceleration;

	/*matice constructor input: */
	int    matice_mid;
	double matice_B;
	double matice_n;
	
	/*penta constructor input: */

	int penta_id;
	int penta_mid;
	int penta_mparid;
	int penta_g[6];
	double penta_h[6];
	double penta_s[6];
	double penta_b[6];
	double penta_k[6];
	int penta_friction_type;
	double penta_p;
	double penta_q;
	int penta_shelf;
	int penta_onbed;
	int penta_onsurface;
	double penta_fit;
	double penta_meanvel;/*!scaling ratio for velocities*/
	double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
	int penta_acceleration;
	int penta_collapse;
	double penta_melting[6];
	double penta_accumulation[6];
	double penta_geothermalflux[6];
	int penta_artdiff;
	int penta_thermal_steadystate;

	/*matpar constructor input: */
	int	matpar_mid;
	double matpar_rho_ice; 
	double matpar_rho_water;
	double matpar_heatcapacity;
	double matpar_thermalconductivity;
	double matpar_latentheat;
	double matpar_beta;
	double matpar_meltingpoint;
	double matpar_mixed_layer_capacity;
	double matpar_thermal_exchange_velocity;
	double matpar_g;
	double matpar_viscosity_overshoot;

	/* node constructor input: */
	int node_id;
	int node_partitionborder=0;
	double node_x[3];
	int node_onbed;
	int node_onsurface;
	int node_numdofs;

		
	#ifdef _PARALLEL_
	/*Metis partitioning: */
	int  range;
	Vec  gridborder;
	int  my_numgrids;
	int* all_numgrids=NULL;
	int  gridcount;
	int  count;
	#endif
	int  first_grid_index;

	/*Rifts:*/
	int* riftsnumpenaltypairs;
	double** riftspenaltypairs;
	int* riftsfill;
	double* riftsfriction;
	double* riftpenaltypairs=NULL;
	int el1,el2;
	 
	/*Penalty partitioning: */
	int num_grids3d_collapsed;
	double* double_penalties_grids3d_collapsed=NULL;
	double* double_penalties_grids3d_noncollapsed=NULL;
	int     grid_ids[6];
	int     num_grid_ids;
	int     grid_id;


	/*First create the elements, nodes and material properties: */
	elements  = new DataSet(ElementsEnum());
	nodes     = new DataSet(NodesEnum());
	materials = new DataSet(MaterialsEnum());


	/*Width of elements: */
	if(strcmp(model->meshtype,"2d")==0){
		elements_width=3; //tria elements
	}
	else{
		elements_width=6; //penta elements
	}


	#ifdef _PARALLEL_
	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
	if(strcmp(model->meshtype,"2d")==0){
		/*load elements: */
		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
	}
	else{
		/*load elements2d: */
		ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
	}

	MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);

	/*Free elements and elements2d: */
	xfree((void**)&model->elements);
	xfree((void**)&model->elements2d);
		

	/*Deal with rifts, they have to be included into one partition only, not several: */
	FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
	
	for(i=0;i<model->numrifts;i++){
		riftpenaltypairs=model->riftspenaltypairs[i];
		for(j=0;j<model->riftsnumpenaltypairs[i];j++){
			el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
			el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
			epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
		}
	}
	/*Free rifts: */
	xfree((void**)&riftsnumpenaltypairs); 
	for(i=0;i<model->numrifts;i++){
		double* temp=riftspenaltypairs[i];
		xfree((void**)&temp);
	}
	xfree((void**)&riftspenaltypairs);
	xfree((void**)&riftsfill);
	xfree((void**)&riftsfriction);

	/*Used later on: */
	my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
	#endif



	/*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */

	/*2d mesh: */
	if (strcmp(model->meshtype,"2d")==0){

		/*Fetch data needed: */
		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
		ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
		ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
		ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
		ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
		ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
		ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
		ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
		ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat");
		ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
		ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
	
		
		for (i=0;i<model->numberofelements;i++){

		#ifdef _PARALLEL_
		/*!All elements have been partitioned above, only create elements for this CPU: */
		if(my_rank==epart[i]){ 
		#endif
			
			
			/*ids: */
			tria_id=i+1; //matlab indexing.
			tria_mid=i+1; //refers to the corresponding material property card
			tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card

			/*vertices offsets: */
			tria_g[0]=(int)*(model->elements+elements_width*i+0);
			tria_g[1]=(int)*(model->elements+elements_width*i+1);
			tria_g[2]=(int)*(model->elements+elements_width*i+2);

			/*thickness,surface and bed:*/
			tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
			tria_h[1]=*(model->thickness+  ((int)*(model->elements+elements_width*i+1)-1)); 
			tria_h[2]=*(model->thickness+  ((int)*(model->elements+elements_width*i+2)-1)) ;

			tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1)); 
			tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1)); 
			tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1)); 

			tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1)); 
			tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1)); 
			tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1)); 

			/*basal drag:*/
			tria_friction_type=(int)model->drag_type;

			tria_k[0]=*(model->drag+        ((int)*(model->elements+elements_width*i+0)-1)); 
			tria_k[1]=*(model->drag+        ((int)*(model->elements+elements_width*i+1)-1)); 
			tria_k[2]=*(model->drag+        ((int)*(model->elements+elements_width*i+2)-1)); 
			
			tria_p=model->p[i];
			tria_q=model->q[i];

			/*element on iceshelf?:*/
			tria_shelf=(int)*(model->elementoniceshelf+i);

			/*type of fitting? : */
			tria_fit=model->fit[0];
			tria_meanvel=model->meanvel;
			tria_epsvel=model->epsvel;

			/*diverse:*/
			tria_acceleration=0;

			/*Create tria element using its constructor:*/
			tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_fit, tria_meanvel, tria_epsvel, tria_acceleration);

			/*Add tria element to elements dataset: */
			elements->AddObject(tria);

			/*Deal with material property card: */
			matice_mid=i+1; //same as the material id from the geom2 elements.
			
			/*Average B over 3 grid elements: */
			B_avg=0;
			for(j=0;j<3;j++){
				B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
			}
			B_avg=B_avg/3;
			matice_B=B_avg;
			matice_n=(double)*(model->n+i);
			
			/*Create matice using its constructor:*/
			matice= new Matice(matice_mid,matice_B,matice_n);
	
			/*Add matice element to materials dataset: */
			materials->AddObject(matice);
	
			#ifdef _PARALLEL_
			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
			 will hold which grids belong to this partition*/
			my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
			#endif

		#ifdef _PARALLEL_
		}//if(my_rank==epart[i])
		#endif

		}//for (i=0;i<numberofelements;i++)

	
		/*Free data : */
		/*xfree((void**)&model->elements);
		xfree((void**)&model->thickness);
		xfree((void**)&model->surface);
		xfree((void**)&model->bed);
		xfree((void**)&model->drag);
		xfree((void**)&model->p);
		xfree((void**)&model->q);
		xfree((void**)&model->elementoniceshelf);
		xfree((void**)&model->fit);
		xfree((void**)&model->B);
		xfree((void**)&model->n);*/
	}
	else{ //	if (strcmp(meshtype,"2d")==0)

		/*Fetch data needed: */
		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
		ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
		ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
		ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
		ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
		ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
		ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
		ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
		ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat");
		ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
		ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
		ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
		ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
		ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
		
		for (i=0;i<model->numberofelements;i++){
		#ifdef _PARALLEL_
		/*We are using our element partition to decide which elements will be created on this node: */
		if(my_rank==epart[i]){
		#endif

			
			/*name and id: */
			penta_id=i+1; //matlab indexing.
			penta_mid=i+1; //refers to the corresponding material property card
			penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card

			/*vertices,thickness,surface,bed and drag: */
			for(j=0;j<6;j++){
				penta_g[j]=(int)*(model->elements+elements_width*i+j);
				penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1)); 
				penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1)); 
				penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1)); 
				penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1)); 
			}

			/*basal drag:*/
			penta_friction_type=(int)model->drag_type;
	
			penta_p=model->p[i];
			penta_q=model->q[i];

			/*diverse: */
			penta_shelf=(int)*(model->elementoniceshelf+i);
			penta_fit=(int)model->fit[0];
			penta_onbed=(int)*(model->elementonbed+i);
			penta_onsurface=(int)*(model->elementonsurface+i);
			penta_meanvel=model->meanvel;
			penta_epsvel=model->epsvel;
			penta_acceleration=0;

			if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
				penta_collapse=1;
			}
			else{
				penta_collapse=0;
			}


			/*Create Penta using its constructor:*/
			penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_fit,penta_meanvel,penta_epsvel,
					penta_acceleration,penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
					penta_thermal_steadystate); 

			/*Add penta element to elements dataset: */
			elements->AddObject(penta);
	

			/*Deal with material:*/
			matice_mid=i+1; //same as the material id from the geom2 elements.
			/*Average B over 6 element grids: */
			B_avg=0;
			for(j=0;j<6;j++){
				B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
			}
			B_avg=B_avg/6;
			matice_B= B_avg;
			matice_n=(double)*(model->n+i);
	
			/*Create matice using its constructor:*/
			matice= new Matice(matice_mid,matice_B,matice_n);
	
			/*Add matice element to materials dataset: */
			materials->AddObject(matice);
			
			#ifdef _PARALLEL_
			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
			 will hold which grids belong to this partition*/
			my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
			my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
			#endif

		#ifdef _PARALLEL_
		}//if(my_rank==epart[i])
		#endif

		}//for (i=0;i<numberofelements;i++)

		/*Free data: */
		/*xfree((void**)&model->elements);
		xfree((void**)&model->thickness);
		xfree((void**)&model->surface);
		xfree((void**)&model->bed);
		xfree((void**)&model->drag);
		xfree((void**)&model->p);
		xfree((void**)&model->q);
		xfree((void**)&model->elementoniceshelf);
		xfree((void**)&model->fit);
		xfree((void**)&model->elementonbed);
		xfree((void**)&model->elementonsurface);
		xfree((void**)&model->elements_type);
		xfree((void**)&model->n);
		xfree((void**)&model->B);*/

	} //if (strcmp(meshtype,"2d")==0)

	/*Add one constant material property to materials: */
	matpar_mid=model->numberofelements+1; //put it at the end of the materials
	matpar_g=model->g; 
	matpar_viscosity_overshoot=model->viscosity_overshoot; 
	matpar_rho_ice=model->rho_ice; 
	matpar_rho_water=model->rho_water; 
	matpar_thermalconductivity=model->thermalconductivity; 
	matpar_heatcapacity=model->heatcapacity; 
	matpar_latentheat=model->latentheat; 
	matpar_beta=model->beta; 
	matpar_meltingpoint=model->meltingpoint; 
	matpar_mixed_layer_capacity=model->mixed_layer_capacity; 
	matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; 

	/*Create matpar object using its constructor: */
	matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
			matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
			matpar_thermal_exchange_velocity,matpar_g,matpar_viscosity_overshoot);
		
	/*Add to materials datset: */
	materials->AddObject(matpar);
	
	#ifdef _PARALLEL_
		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
		 *element partition, and which are on its border with other nodes:*/
		gridborder=NewVec(model->numberofnodes);

		for (i=0;i<model->numberofnodes;i++){
			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
		}
		VecAssemblyBegin(gridborder);
		VecAssemblyEnd(gridborder);

		#ifdef _DEBUG_
		VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
		#endif
		
		VecToMPISerial(&my_bordergrids,gridborder);

		#ifdef _DEBUG_
		if(my_rank==0){
			for (i=0;i<model->numberofnodes;i++){
				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
			}
		}
		#endif
	#endif

	/*Partition penalties in 3d: */
	if(strcmp(model->meshtype,"3d")==0){
	
		/*Get penalties: */
		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");

		if(model->numpenalties){

			model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
			for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;

			for(i=0;i<model->numpenalties;i++){
				first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
					/*All grids that are being penalised belong to this node's internal grid partition.:*/
					model->penaltypartitioning[i]=1;
				}
				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
					model->penaltypartitioning[i]=0;
				}
			}
		}

		/*Free penalties: */
		xfree((void**)&model->penalties);
	}

	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
	 We can therefore determine  which grids are internal to this node's partition 
	 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
	 that, go and create the grids*/

	/*Create nodes from x,y,z, as well as the spc values on those grids: */
		
	/*First fetch data: */
	if (strcmp(model->meshtype,"3d")==0){
		ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
	}
	ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
	ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
	ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");

	/*Get analysis_type: */
	analysis_type=AnalysisTypeAsEnum(model->analysis_type);

	/*Get number of dofs per node: */
	DistributeNumDofs(&node_numdofs,analysis_type);

	for (i=0;i<model->numberofnodes;i++){
	#ifdef _PARALLEL_
	/*keep only this partition's nodes:*/
	if((my_grids[i]==1)){
	#endif

		node_id=i+1; //matlab indexing
			
		
		
		#ifdef _PARALLEL_
		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
			node_partitionborder=1;
		}
		else{
			node_partitionborder=0;
		}
		#else
			node_partitionborder=0;
		#endif


		node_x[0]=model->x[i];
		node_x[1]=model->y[i];
		node_x[2]=model->z[i];
		
		node_onbed=(int)model->gridonbed[i];
		node_onsurface=(int)model->gridonsurface[i];

		/*Create node using its constructor: */
		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface);

		/*set single point constraints.: */
		if (strcmp(model->meshtype,"3d")==0){
			/*On a 3d mesh, we may have collapsed grids. Spc all their dofs: */
			if (model->deadgrids[i]){
				for(k=0;k<node_numdofs;k++){
					node->DofInSSet(k);
				}
			}
		}
		/*Add node to nodes dataset: */
		nodes->AddObject(node);

	#ifdef _PARALLEL_
	} //if((my_grids[i]==1))
	#endif
	}

	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
	 * datasets, it will not be redone: */
	elements->Presort();
	nodes->Presort();
	materials->Presort();

	/*Clean fetched data: */
	xfree((void**)&model->deadgrids);
	xfree((void**)&model->x);
	xfree((void**)&model->y);
	xfree((void**)&model->z);
	xfree((void**)&model->gridonbed);
	xfree((void**)&model->gridonsurface);
		
	cleanup_and_return:

	/*Assign output pointer: */
	*pelements=elements;
	*pnodes=nodes;
	*pmaterials=materials;

	/*Keep partitioning information into model*/
	model->epart=epart;
	model->my_grids=my_grids;
	model->my_bordergrids=my_bordergrids;

	/*Free ressources:*/
	#ifdef _PARALLEL_
	xfree((void**)&all_numgrids);
	VecFree(&gridborder);
	#endif

	return noerr;

}
