/*!\file Pengrid.c
 * \brief: implementation of the Pengrid 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 "./Pengrid.h"
#include <string.h>
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../shared/shared.h"
#include "../include/typedefs.h"

		
Pengrid::Pengrid(){
	return;
}

Pengrid::Pengrid(int	pengrid_id, int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){
	
	id=pengrid_id;
	dof=pengrid_dof;
	active=pengrid_active;
	penalty_offset =pengrid_penalty_offset;
	thermal_steadystate=pengrid_thermal_steadystate;

	node_id=pengrid_node_id;
	node_offset=UNDEF;
	node=NULL;

	return;
}

Pengrid::~Pengrid(){
	return;
}
		
void Pengrid::Echo(void){

	printf("Pengrid:\n");
	printf("   id: %i\n",id);
	printf("   dof: %i\n",dof);
	printf("   active: %i\n",active);
	printf("   penalty_offset: %g\n",penalty_offset);
	printf("   thermal_steadystate: %i\n",thermal_steadystate);
	printf("   node_id: [%i]\n",node_id);
	printf("   node_offset: [%i]\n",node_offset);
	
	if(node)node->Echo();
	return;
}
		
void  Pengrid::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;

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

	/*get enum type of Pengrid: */
	enum_type=PengridEnum();
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	
	/*marshall Pengrid data: */
	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
	memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
	memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
	memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
	memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
	memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);

	*pmarshalled_dataset=marshalled_dataset;
	return;
}
		
int   Pengrid::MarshallSize(){

	return sizeof(id)+
		sizeof(dof)+
		sizeof(active)+
		sizeof(penalty_offset)+
		sizeof(thermal_steadystate)+
		sizeof(node_id)+
		sizeof(node_offset)+
		sizeof(int); //sizeof(int) for enum type
}

char* Pengrid::GetName(void){
	return "pengrid";
}
		

void  Pengrid::Demarshall(char** pmarshalled_dataset){

	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(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
	memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
	memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
	memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
	memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
	memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);

	node=NULL;

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
int Pengrid::Enum(void){

	return PengridEnum();
}

int    Pengrid::GetId(void){ return id; }

int    Pengrid::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;}

#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::Configure"

void  Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){

	DataSet* nodesin=NULL;

	/*Recover pointers :*/
	nodesin=(DataSet*)pnodesin;

	/*Link this load with its nodes: */
	ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);

}


#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::CreateKMatrix"

void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){

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

}

#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::CreatePVector"
void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type){

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

}
#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::UpdateFromInputs"
void  Pengrid::UpdateFromInputs(void* inputs){
	
}

#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){

	if ((analysis_type==DiagnosticStokesAnalysisEnum())){

		PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type);

	}
	else{
		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
	}

}

#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){
	
	const int numgrids=1;
	const int NDOF3=3;
	const int numdof=numgrids*NDOF3;
	int       doflist[numdof];
	int       numberofdofspernode;

	int dofs1=0;
	int dofs2=1;
	double slope[2];
	int found=0;
	double Ke[4][4]={0.0};
	
	ParameterInputs* inputs=NULL;

	/*recover pointers: */
	inputs=(ParameterInputs*)vinputs;

	/*Get dof list: */
	GetDofList(&doflist[0],&numberofdofspernode);
	
	/*recover slope: */
	found=inputs->Recover("bedslopex",&slope[0],1,&dofs1,numgrids,(void**)&node);
	if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
	found=inputs->Recover("bedslopey",&slope[1],1,&dofs2,numgrids,(void**)&node);
	if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!");

	//Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
	Ke[0][0]=-slope[0]*kmax*pow(10.0,penalty_offset);
	Ke[1][1]=-slope[1]*kmax*pow(10.0,penalty_offset);
	Ke[2][2]=kmax*pow(10,penalty_offset);
	
	/*Add Ke to global matrix Kgg: */
	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
}

#undef __FUNCT__ 
#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
void  Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){

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

}

Object* Pengrid::copy() {
	return new Pengrid(*this); 
}


void  Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){

	int j;
	int doflist_per_node[MAXDOFSPERNODE];
	int numberofdofspernode;
	
	node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
	for(j=0;j<numberofdofspernode;j++){
		doflist[j]=doflist_per_node[j];
	}

	/*Assign output pointers:*/
	*pnumberofdofspernode=numberofdofspernode;
}