/*
 * CreateElementsNodesAndMaterialsStressbalanceHoriz.c:
 */

#include "../../toolkits/toolkits.h"
#include "../../classes/classes.h"
#include "../../shared/shared.h"
#include "./ModelProcessorx.h"

void	CreateElementsVerticesAndMaterials(Elements** pelements,Vertices** pvertices,Materials** pmaterials, IoModel* iomodel,const int nummodels){

	/*Intermediary*/
	int i;
	int materials_type;
	bool control_analysis;
	bool dakota_analysis;

	/*DataSets: */
	Elements  *elements  = NULL;
	Vertices  *vertices  = NULL;
	Materials *materials = NULL;

	/*Fetch parameters: */
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
	iomodel->Constant(&materials_type,MaterialsEnum);

	/*Did we already create the elements? : */
	if(*pelements)return;

	/*First create the elements, vertices, nodes and material properties, if they don't already exist */
	elements  = new Elements();
	vertices  = new Vertices();
	materials = new Materials();

	#ifdef _HAVE_3D_
	if(iomodel->dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum);
	#endif
	if(control_analysis)iomodel->FetchData(3,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum);

	/*Create elements*/
	for(i=0;i<iomodel->numberofelements;i++){
		if(iomodel->my_elements[i]){

			/*Create and add tria element to elements dataset: */
			if(iomodel->dim==2) elements->AddObject(new Tria(i+1,i,i,iomodel,nummodels));
	        #ifdef _HAVE_3D_
			else       elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels));
	        #endif
		}
	}

	/*Create materials*/
	switch(materials_type){
		case MaticeEnum:
			iomodel->FetchData(2,MaterialsRheologyBEnum,MaterialsRheologyNEnum);
			for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
			if(dakota_analysis){
				if(iomodel->dim==2) materials->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum); 
           #ifdef _HAVE_3D_
				else       materials->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); 
           #endif
			}
			break;
		case MatdamageiceEnum:
			iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum);
			for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matdamageice(i+1,i,iomodel));
			break;
		default:
			_error_("Materials "<<EnumToStringx(materials_type)<<" not supported");
	}

	/*Free data: */
	iomodel->DeleteData(8,MeshUpperelementsEnum,MeshLowerelementsEnum,
				MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum,
				InversionMaxParametersEnum);

	/*Add new constant material property to materials, at the end: */
	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials

	/*Create vertices: */

	/*Fetch data:*/
	iomodel->FetchData(6,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum,MaskIceLevelsetEnum);
	CreateNumberNodeToElementConnectivity(iomodel);

	for(i=0;i<iomodel->numberofvertices;i++){

		/*vertices and nodes (same number, as we are running continuous galerkin formulation): */
		if(iomodel->my_vertices[i]){

			/*Add vertex to vertices dataset: */
			vertices->AddObject(new Vertex(i+1,i,i,iomodel));
		}
	}

	/*Free data: */
	iomodel->DeleteData(6,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum,MaskIceLevelsetEnum);

	/*Assign output pointer: */
	*pelements=elements;
	*pvertices=vertices;
	*pmaterials=materials;
}
