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

/*Headers*/
/*{{{1*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../include/include.h"
#include "../../shared/shared.h"
/*}}}*/

/*Element macros*/
#define NUMVERTICES 2

/*Penpair constructors and destructor*/
/*FUNCTION Penpair::constructor {{{1*/
Penpair::Penpair(){

	this->hnodes=NULL;
	this->nodes=NULL;
	this->parameters=NULL;
	return;
}
/*}}}1*/
/*FUNCTION Penpair::creation {{{1*/
Penpair::Penpair(int penpair_id, int* penpair_node_ids,int in_analysis_type){
	
	this->id=penpair_id;
	this->analysis_type=in_analysis_type;
	this->hnodes=new Hook(penpair_node_ids,2);
	this->parameters=NULL;
	this->nodes=NULL;
	
	return;
}
/*}}}1*/
/*FUNCTION Penpair::destructor {{{1*/
Penpair::~Penpair(){
	delete hnodes;
	return;
}
/*}}}1*/

/*Object virtual functions definitions:*/
/*FUNCTION Penpair::Echo {{{1*/
void Penpair::Echo(void){

	int i;

	printf("Penpair:\n");
	printf("   id: %i\n",id);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	hnodes->Echo();
	
	return;
}
/*}}}1*/
/*FUNCTION Penpair::DeepEcho {{{1*/
void Penpair::DeepEcho(void){

	printf("Penpair:\n");
	printf("   id: %i\n",id);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	hnodes->DeepEcho();

	return;
}		
/*}}}1*/
/*FUNCTION Penpair::Id {{{1*/
int    Penpair::Id(void){ return id; }
/*}}}1*/
/*FUNCTION Penpair::MyRank {{{1*/
int    Penpair::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}1*/
#ifdef _SERIAL_
/*FUNCTION Penpair::Marshall {{{1*/
void  Penpair::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum type of Penpair: */
	enum_type=PenpairEnum;
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	
	/*marshall Penpair data: */
	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);

	/*Marshall hooks*/
	hnodes->Marshall(&marshalled_dataset);

	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}1*/
/*FUNCTION Penpair::MarshallSize {{{1*/
int   Penpair::MarshallSize(){

	return sizeof(id)+
		+sizeof(analysis_type)
		+hnodes->MarshallSize()
		+sizeof(int); //sizeof(int) for enum type
}
/*}}}1*/
/*FUNCTION Penpair::Demarshall {{{1*/
void  Penpair::Demarshall(char** pmarshalled_dataset){

	int i;
	char* marshalled_dataset=NULL;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/
	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);

	/*demarshall hooks: */
	hnodes=new Hook(); hnodes->Demarshall(&marshalled_dataset);

	/*pointers are garbabe, until configuration is carried out: */
	nodes=NULL;
	parameters=NULL;

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}1*/
#endif
/*FUNCTION Penpair::ObjectEnum{{{1*/
int Penpair::ObjectEnum(void){

	return PenpairEnum;
}
/*}}}1*/
/*FUNCTION Penpair::copy {{{1*/
Object* Penpair::copy() {
	
	Penpair* penpair=NULL;

	penpair=new Penpair();

	/*copy fields: */
	penpair->id=this->id;
	penpair->analysis_type=this->analysis_type;

	/*now deal with hooks and objects: */
	penpair->hnodes=(Hook*)this->hnodes->copy();
	penpair->nodes =(Node**)penpair->hnodes->deliverp();

	/*point parameters: */
	penpair->parameters=this->parameters;

	return penpair;

}
/*}}}*/
		
/*Load virtual functions definitions:*/
/*FUNCTION Penpair::Configure {{{1*/
void  Penpair::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){

	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
	 * datasets, using internal ids and offsets hidden in hooks: */
	hnodes->configure(nodesin);

	/*Initialize hooked fields*/
	this->nodes  =(Node**)hnodes->deliverp();

	/*point parameters to real dataset: */
	this->parameters=parametersin;

}
/*}}}1*/
/*FUNCTION Penpair::SetCurrentConfiguration {{{1*/
void  Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){

}
/*}}}1*/
/*FUNCTION Penpair::CreateKMatrix {{{1*/
void  Penpair::CreateKMatrix(Mat Kff, Mat Kfs){
	/*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
	/*No loads applied, do nothing: */
	return;

}
/*}}}1*/
/*FUNCTION Penpair::CreatePVector {{{1*/
void  Penpair::CreatePVector(Vec pf){

	/*No loads applied, do nothing: */
	return;

}
/*}}}1*/
/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
void  Penpair::CreateJacobianMatrix(Mat Jff){
	this->CreateKMatrix(Jff,NULL);
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/
void  Penpair::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){

	/*Retrieve parameters: */
	ElementMatrix* Ke=NULL;
	int analysis_type;
	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	switch(analysis_type){
		case DiagnosticHorizAnalysisEnum:
			Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax);
			break;
		case PrognosticAnalysisEnum:
			Ke=PenaltyCreateKMatrixPrognostic(kmax);
			break;
		default:
			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}

	/*Add to global Vector*/
	if(Ke){
		Ke->AddToGlobal(Kff,Kfs);
		delete Ke;
	}
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreatePVector {{{1*/
void  Penpair::PenaltyCreatePVector(Vec pf,double kmax){
	/*No loads applied, do nothing: */
	return;
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreateJacobianMatrix{{{1*/
void  Penpair::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
}
/*}}}1*/
/*FUNCTION Penpair::InAnalysis{{{1*/
bool Penpair::InAnalysis(int in_analysis_type){
	if (in_analysis_type==this->analysis_type)return true;
	else return false;
}
/*}}}*/

/*Update virtual functions definitions:*/
/*FUNCTION Penpair::InputUpdateFromConstant(double constant, int name) {{{1*/
void  Penpair::InputUpdateFromConstant(double constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Penpair::InputUpdateFromConstant(int constant, int name) {{{1*/
void  Penpair::InputUpdateFromConstant(int constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Penpair::InputUpdateFromConstant(bool constant, int name) {{{1*/
void  Penpair::InputUpdateFromConstant(bool constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Penpair::InputUpdateFromVector(double* vector, int name, int type) {{{1*/
void  Penpair::InputUpdateFromVector(double* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Penpair::InputUpdateFromVector(int* vector, int name, int type) {{{1*/
void  Penpair::InputUpdateFromVector(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Penpair::InputUpdateFromVector(bool* vector, int name, int type) {{{1*/
void  Penpair::InputUpdateFromVector(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/

/*Penpair management:*/
/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticHoriz{{{1*/
ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticHoriz(double kmax){

	int    approximation0=nodes[0]->GetApproximation();
	int    approximation1=nodes[1]->GetApproximation();

	switch(approximation0){
		case MacAyealApproximationEnum:
			switch(approximation1){
				case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); 
				case PattynApproximationEnum:   return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); 
				default: _error_("not supported yet");
			}
		case PattynApproximationEnum:
			switch(approximation1){
				case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); 
				case PattynApproximationEnum:   return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); 
				default: _error_("not supported yet");
			}
		case StokesApproximationEnum:
			switch(approximation1){
				case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); 
				case NoneApproximationEnum: return   PenaltyCreateKMatrixDiagnosticStokes(kmax); 
				default: _error_("not supported yet");
			}
		case NoneApproximationEnum:
			switch(approximation1){
				case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); 
				case NoneApproximationEnum: return   PenaltyCreateKMatrixDiagnosticStokes(kmax); 
				default: _error_("not supported yet");
			}
		default: _error_("not supported yet");
	}
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn {{{1*/
ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn(double kmax){
	
	const int numdof=NUMVERTICES*NDOF2;
	double penalty_offset;

	/*Initialize Element vector and return if necessary*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);

	/*recover parameters: */
	parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);

	//Create elementary matrix: add penalty to 
	Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[0*numdof+2]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[2*numdof+0]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset);

	Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[1*numdof+3]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[3*numdof+1]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset);

	/*Clean up and return*/
	return Ke;
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticStokes(double kmax){
	
	const int numdof=NUMVERTICES*NDOF4;
	double penalty_offset;

	/*Initialize Element vector and return if necessary*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);

	/*recover parameters: */
	parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);

	//Create elementary matrix: add penalty to 
	Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[0*numdof+4]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[4*numdof+0]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[4*numdof+4]=+kmax*pow((double)10.0,penalty_offset);

	Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[1*numdof+5]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[5*numdof+1]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[5*numdof+5]=+kmax*pow((double)10.0,penalty_offset);
	
	Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[2*numdof+6]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[6*numdof+2]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[6*numdof+6]=+kmax*pow((double)10.0,penalty_offset);

	Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset);
	Ke->values[3*numdof+7]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[7*numdof+3]=-kmax*pow((double)10.0,penalty_offset);
	Ke->values[7*numdof+7]=+kmax*pow((double)10.0,penalty_offset);

	/*Clean up and return*/
	return Ke;
}
/*}}}1*/
/*FUNCTION Penpair::PenaltyCreateKMatrixPrognostic {{{1*/
ElementMatrix* Penpair::PenaltyCreateKMatrixPrognostic(double kmax){

	const int numdof=NUMVERTICES*NDOF1;
	double penalty_factor;

	/*Initialize Element vector and return if necessary*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);

	/*recover parameters: */
	parameters->FindParam(&penalty_factor,PrognosticPenaltyFactorEnum);

	//Create elementary matrix: add penalty to 
	Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_factor);
	Ke->values[0*numdof+1]=-kmax*pow((double)10.0,penalty_factor);
	Ke->values[1*numdof+0]=-kmax*pow((double)10.0,penalty_factor);
	Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_factor);

	/*Clean up and return*/
	return Ke;
}
/*}}}1*/