/*
 * CreateElementsNodesAndMaterialsDiagnosticHoriz.c:
 */

#include "../../Container/Container.h"
#include "../../toolkits/toolkits.h"
#include "../../io/io.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../classes/objects/objects.h"
#include "../../shared/shared.h"
#include "../MeshPartitionx/MeshPartitionx.h"
#include "../../include/include.h"
#include "./ModelProcessorx.h"

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

	/*Intermediary*/
	int i;
	int dim,materials_type;
	int numberofelements;
	int numberofvertices;
	bool control_analysis;

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

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	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();

	/*First, partition elements and vertices. Nodes will partitioned on a per analysis_type basis. If partitining already done, ignore: */
	ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel);

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

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

			/*Create and add tria element to elements dataset: */
			if(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<numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
			break;
		case MatdamageiceEnum:
			iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum);
			for (i=0;i<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(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum,
				MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum,
				InversionMaxParametersEnum);

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

	/*Create vertices: */

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

	for (i=0;i<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,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum);

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