/*
 * 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,j,k,n;
	int dim;
	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);

	/*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);
	
	/*Fetch data needed: */
	iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum);
	#ifdef _HAVE_3D_
	if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum);
	#endif
	if(control_analysis)iomodel->FetchData(3,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum);
	
	/*Create elements and materials: */
	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 and add material property to materials dataset: */
			materials->AddObject(new Matice(i+1,i,iomodel));
		}
	}
	
	/*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;
}
