/*!\file FemModel.c
 * \brief: implementation of the FemModel object
 */


#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "stdio.h"
#include "../Container/Container.h"
#include "../modules/ModelProcessorx/ModelProcessorx.h"
#include "./objects.h"
#include "../include/include.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../modules/modules.h"

/*Object constructors and destructor*/
/*FUNCTION FemModel::constructor {{{1*/
FemModel::FemModel(ConstDataHandle IOMODEL,const int in_solution_type,const int* analyses,const int nummodels){

	/*intermediary*/
	int i;
	IoModel* iomodel=NULL;
	int analysis_type;

	/*Initialize internal data: */
	this->nummodels=nummodels;
	this->solution_type=in_solution_type;
	this->analysis_counter=nummodels-1; //point to last analysis_type carried out.
	this->results=new DataSet(); //not initialized by CreateDataSets
	this->partition=NULL;
	this->tpartition=NULL;
	
	/*Dynamically allocate whatever is a list of length nummodels: */
	analysis_type_list=(int*)xmalloc(nummodels*sizeof(int));
	m_Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat));
	m_Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat));
	m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
	m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
	m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec));

	/*Initialize: */
	for(i=0;i<nummodels;i++)analysis_type_list[i]=analyses[i];
	for(i=0;i<nummodels;i++)m_Rmg[i]=NULL;
	for(i=0;i<nummodels;i++)m_Gmn[i]=NULL;
	for(i=0;i<nummodels;i++)m_nodesets[i]=NULL;
	for(i=0;i<nummodels;i++)m_yg[i]=NULL;
	for(i=0;i<nummodels;i++)m_ys[i]=NULL;

	/*create datasets for all analyses: */
	ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,IOMODEL,this->solution_type,nummodels,analyses);

	/*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
	for(i=0;i<nummodels;i++){

		_printf_("   processing finite element model of analysis %s:\n",EnumAsString(analysis_type_list[i]));
		analysis_type=analysis_type_list[i];
		this->SetCurrentAnalysis(analysis_type);
	
		_printf_("      create degrees of freedom\n");
		VerticesDofx(&partition,&tpartition,vertices,parameters);
		NodesDofx(nodes,parameters,analysis_type);

		_printf_("      create single point constraints\n");
		SpcNodesx( &m_yg[i], nodes,constraints,analysis_type); 

		_printf_("      create rigid body constraints\n");
		MpcNodesx( &m_Rmg[i], nodes,constraints,analysis_type); 

		_printf_("      create node sets\n");
		BuildNodeSetsx(&m_nodesets[i], nodes,analysis_type);

		_printf_("      reducing single point constraints vector\n");
		Reducevectorgtosx(&m_ys[i], m_yg[i],m_nodesets[i]);

		_printf_("      normalizing rigid body constraints matrix\n");
		NormalizeConstraintsx(&m_Gmn[i], m_Rmg[i],m_nodesets[i]);

		_printf_("      configuring element and loads\n");
		ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters);
	}
}
/*}}}1*/
/*FUNCTION FemModel::destructor {{{1*/
FemModel::~FemModel(){

	/*Intermediary*/
	int i;

	/*Delete all the datasets: */
	xfree((void**)&analysis_type_list);
	delete elements;
	delete nodes;
	delete vertices;
	delete constraints;
	delete loads;
	delete materials;
	delete parameters;
	delete results;
	VecFree(&partition);
	VecFree(&tpartition);

	for(i=0;i<nummodels;i++){
		Mat temp_Rmg=m_Rmg[i];
		MatFree(&temp_Rmg);
		Mat temp_Gmn=m_Gmn[i];
		MatFree(&temp_Gmn);
		NodeSets* temp_nodesets=m_nodesets[i];
		delete temp_nodesets;
		Vec temp_yg=m_yg[i];
		VecFree(&temp_yg);
		Vec temp_ys=m_ys[i];
		VecFree(&temp_ys);
	}

	/*Delete dynamically allocated arrays: */
	xfree((void**)&m_Rmg);
	xfree((void**)&m_Gmn);
	xfree((void**)&m_nodesets);
	xfree((void**)&m_yg);
	xfree((void**)&m_ys);

}
/*}}}1*/

/*Object management*/
/*FUNCTION FemModel::Echo {{{1*/
void FemModel::Echo(void){

	printf("FemModel echo: \n");
	printf("   number of fem models: %i\n",nummodels);
	printf("   analysis_type_list: \n");
	for(int i=0;i<nummodels;i++)printf("     %i: %s\n",i,EnumAsString(analysis_type_list[i]));
	printf("   current analysis_type: \n");
	printf("     %i: %s\n",analysis_counter,EnumAsString(analysis_type_list[analysis_counter]));


}
/*}}}*/

/*Numerics: */
/*FUNCTION FemModel::GetCurrentAnalysis {{{1*/
int FemModel::GetCurrentAnalysis(){
	return analysis_type_list[analysis_counter];
}
/*}}}1*/
/*FUNCTION FemModel::SetCurrentAnalysis {{{1*/
void FemModel::SetCurrentAnalysis(int analysis_type){

	int found=-1;
	for(int i=0;i<nummodels;i++){
		if (analysis_type_list[i]==analysis_type){
			found=i;
			break;
		}
	}
	if(found!=-1) analysis_counter=found;
	else ISSMERROR("Could not find analysis_type %s in list of FemModel analyses",EnumAsString(analysis_type));

	/*activate matrices/vectors: */
	Rmg=m_Rmg[analysis_counter];
	Gmn=m_Gmn[analysis_counter];
	nodesets=m_nodesets[analysis_counter];
	yg=m_yg[analysis_counter];
	ys=m_ys[analysis_counter];

	/*Now, plug analysis_counter and analysis_type inside the parameters: */
	this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
	this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
}
/*}}}1*/
/*FUNCTION FemModel::SetCurrentAnalysisAlias {{{1*/
void FemModel::SetCurrentAnalysisAlias(int base_analysis_type,int real_analysis_type){

	/*Use base_analysis_type to setup the analysis counter, but the analysis type of the FemModel will remain 
	 * real_analysis_type. This means we are using the base_analysis_type settings to run a similar, compatible 
	 * analysis called real_analysis_type: */

	int found=-1;
	for(int i=0;i<nummodels;i++){
		if (analysis_type_list[i]==base_analysis_type){
			found=i;
			break;
		}
	}
	if(found!=-1) analysis_counter=found;
	else ISSMERROR("Could not find alias for analysis_type %s in list of FemModel analyses",EnumAsString(base_analysis_type));

	/*activate matrices/vectors: */
	Rmg=m_Rmg[analysis_counter];
	Gmn=m_Gmn[analysis_counter];
	nodesets=m_nodesets[analysis_counter];
	yg=m_yg[analysis_counter];
	ys=m_ys[analysis_counter];

	/*Now, plug analysis_counter and real_analysis_type inside the parameters: */
	this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
	this->parameters->SetParam(real_analysis_type,AnalysisTypeEnum);
}
/*}}}1*/