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

		
Icefront::Icefront(){
	return;
}

Icefront::Icefront(char icefront_type[ICEFRONTSTRING],int icefront_sid, int icefront_mparid, int icefront_eid, int icefront_element_type, 
		int icefront_node_ids[MAX_ICEFRONT_GRIDS],double icefront_h[MAX_ICEFRONT_GRIDS],double	icefront_b[MAX_ICEFRONT_GRIDS]){

	int i;
	
	strcpy(type,icefront_type);
	sid=icefront_sid;
	
	mparid=icefront_mparid;
	
	eid=icefront_eid;
	element_type=icefront_element_type;
	for(i=0;i<MAX_ICEFRONT_GRIDS;i++){
		node_ids[i]=icefront_node_ids[i];
		node_offsets[i]=UNDEF;
		nodes[i]=NULL;
		h[i]=icefront_h[i];
		b[i]=icefront_b[i];
	}
	
	matpar=NULL;
	matpar_offset=UNDEF;

	element=NULL;
	element_offset=UNDEF;

	return;
}

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

	int i;
	
	printf("Icefront:\n");
	printf("   type: %s\n",type);
	printf("   sid: %i\n",sid);
	
	printf("   mparid: %i\n",mparid);
	printf("   matpar_offset: %i\n",matpar_offset);
	printf("   matpar: \n");
	if(matpar)matpar->Echo();
	
	printf("   element_type: %i\n",element_type);
	printf("   eid: %i\n",eid);
	printf("   element_offset: %i\n",eid);
	printf("   element\n");
	if(element)element->Echo();
		
	if (strcmp(type,"segment")==0){
		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
		printf("   h=[%g,%g]\n",h[0],h[1]);
		printf("   b=[%g,%g]\n",b[0],b[1]);
		for(i=0;i<2;i++){
			if(nodes[i])nodes[i]->Echo();
		}
	}
	else{
		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
		printf("   h=[%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3]);
		printf("   b=[%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3]);
		for(i=0;i<4;i++){
			if(nodes[i])nodes[i]->Echo();
		}
	}

	return;
}
		
void  Icefront::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;

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

	/*get enum type of Icefront: */
	enum_type=IcefrontEnum();
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	
	/*marshall Icefront data: */
	memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
	
	memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
	memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
	
	memcpy(marshalled_dataset,&element_type,sizeof(element_type));marshalled_dataset+=sizeof(element_type);
	memcpy(marshalled_dataset,&eid,sizeof(eid));marshalled_dataset+=sizeof(eid);
	memcpy(marshalled_dataset,&element_offset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);
	
	memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
	memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
	
	memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h);
	memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b);

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

	return sizeof(type)+
		sizeof(sid)+
		sizeof(mparid)+
		sizeof(matpar_offset)+
		sizeof(eid)+
		sizeof(element_offset)+
		sizeof(element_type)+
		sizeof(node_ids)+
		sizeof(node_offsets)+
		sizeof(h)+
		sizeof(b)+
		sizeof(int); //sizeof(int) for enum type
}

char* Icefront::GetName(void){
	return "icefront";
}
		

void  Icefront::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(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
	
	memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
	memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);

	memcpy(&element_type,marshalled_dataset,sizeof(element_type));marshalled_dataset+=sizeof(element_type);
	memcpy(&eid,marshalled_dataset,sizeof(eid));marshalled_dataset+=sizeof(eid);
	memcpy(&element_offset,marshalled_dataset,sizeof(element_offset));marshalled_dataset+=sizeof(element_offset);

	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
	
	memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
	memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);

	for(i=0;i<MAX_ICEFRONT_GRIDS;i++)nodes[i]=NULL;
	matpar=NULL;
	element=NULL;

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}

int Icefront::Enum(void){

	return IcefrontEnum();

}

int    Icefront::GetId(void){ return sid; }

int    Icefront::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
		
#undef __FUNCT__ 
#define __FUNCT__ "Icefront::Configure"

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

	DataSet* elementsin=NULL;
	DataSet* nodesin=NULL;
	DataSet* materialsin=NULL;

	/*Recover pointers :*/
	elementsin=(DataSet*)pelementsin;
	nodesin=(DataSet*)pnodesin;
	materialsin=(DataSet*)pmaterialsin;
	
	/*Link this load with its nodes: */
	if (strcmp(type,"segment")==0){
		
		ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
	
	}
	else{
		ResolvePointers((Object**)nodes,node_ids,node_offsets,4,nodesin);
	}

	
	/*Same for materials: */
	ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);


	/*Same for element: */
	ResolvePointers((Object**)&element,&eid,&element_offset,1,elementsin);

}



#undef __FUNCT__ 
#define __FUNCT__ "Icefront::CreateKMatrix"

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

	/*No stiffness loads applied, do nothing: */

	return;

}

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

	/*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
	}
	else{
		throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
	}
}


#undef __FUNCT__ 
#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz"
void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){

	/*Branck on the type of icefront: */
	if (strcmp(type,"segment")==0){
		CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type);
	}
	else{
		CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type);
	}
}	

#undef __FUNCT__ 
#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment"

void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,ParameterInputs* inputs, int analysis_type){

	int i,j;
	
	const int numgrids=2;
	const int NDOF2=2;
	int   numberofdofspernode;
	const int numdofs=numgrids*NDOF2;
	int   doflist[numdofs];
	double xyz_list[numgrids][3];

	double x1,y1,x2,y2;

	/* input parameters: */
	double thickness_list_element[3];
	double thickness_list[2];
	double bed_list_element[3];
	double bed_list[2];
	
	int    grid1,grid2;
	double normal[2];
	double length;
	int    fill;

	/*Objects: */
	double  pe_g[numdofs];
	Matpar* matpar=NULL;
	Node**  element_nodes=NULL;

	/*Recover material and fill parameters: */
	matpar=element->GetMatPar();
	fill=element->GetShelf();

	//check that the element is onbed (collapsed formulation) otherwise:pe=0
	if(element_type==PentaEnum()){
		if  (!element->GetOnBed()){
			return;
		}
	}

	/* Set pe_g to 0: */
	for(i=0;i<numdofs;i++) pe_g[i]=0.0;
		
	/*Identify which grids are comprised in the segment. */
	if(element_type==TriaEnum())element_nodes=(Node**)xmalloc(3*sizeof(Node*));
	if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*));
	element->GetNodes(element_nodes);
		
	/*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */
	for(i=0;i<3;i++){
		if (nodes[0]==element_nodes[i])grid1=i;
		if (nodes[1]==element_nodes[i])grid2=i;
	}
	
	/* Get dof list and node coordinates: */
	GetDofList(&doflist[0],&numberofdofspernode);
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	
	/*Now build thickness_list_element and bed_list_element: */
	element->GetThicknessList(&thickness_list_element[0]);
	element->GetBedList(&thickness_list_element[0]);
	
	/*Build thickness_list and bed_list: */
	thickness_list[0]=thickness_list_element[grid1];
	thickness_list[1]=thickness_list_element[grid2];
	bed_list[0]=bed_list_element[grid1];
	bed_list[1]=bed_list_element[grid2];

	//Recover grid coordinates
	x1=xyz_list[0][0];
	y1=xyz_list[0][1];
	x2=xyz_list[1][0];
	y2=xyz_list[1][1];

	//Compute length and normal of segment
	normal[0]=cos(atan2(x1-x2,y2-y1));
	normal[1]=sin(atan2(x1-x2,y2-y1));
	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));

	//Compute load contribution for this segment:
	SegmentPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list,bed_list,normal,length,fill);
		
	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && this->icefront.eid==ELID){ 
		printf("\nicefront load\n");
		printf("grids %i %i\n",grid1,grid2);
		printf("rho_water %g\n",rho_water);
		printf("rho_ice %g\n",rho_ice);
		printf("gravity %g\n",gravity);
		printf("thickness_list (%g,%g)\n",thickness_list[0],thickness_list[1]);
		printf("bed_list (%g,%g)\n",bed_list[0],bed_list[1]);
		printf("normal (%g,%g)\n",normal[0],normal[1]);
		printf("length %g\n",length);
		printf("fill %i\n",fill);
		printf("pe_g->terms\n");
		for(i=0;i<numgrids*NDOF2;i++){
			printf("%g ",*(pe_g->terms+i));
		}
	}
	#endif

	/*Plug pe_g into vector: */
	VecSetValues(pg,numdofs,doflist,pe_g,INSERT_VALUES);

	/*Free ressources:*/
	xfree((void**)&element_nodes);

}


#undef __FUNCT__ 
#define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad"
void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type){

	throw ErrorException(__FUNCT__," not implemented yet!");
}

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

#undef __FUNCT__ 
#define __FUNCT__ "Icefront::GetDofList"

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

	int i,j;
	int doflist_per_node[MAXDOFSPERNODE];
	int numberofdofspernode;
	
	if (strcmp(type,"segment")==0){
		for(i=0;i<2;i++){
			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
			for(j=0;j<numberofdofspernode;j++){
				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
			}
		}
	}
	else{
		for(i=0;i<4;i++){
			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
			for(j=0;j<numberofdofspernode;j++){
				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
			}
		}
	}

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

}


#undef __FUNCT__ 
#define __FUNCT__ "Icefront::SegmentPressureLoad"

void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill){

	double   nx,ny;
	double   h1,h2,b1,b2;
	int      num_gauss;
	double*  segment_gauss_coord=NULL;
	double*  gauss_weights=NULL;
	int      ig;
	double   Jdet;
	double   thickness,bed;
	double   ice_pressure,water_pressure,air_pressure;
	double   surface_under_water,base_under_water;
	double   pressure;

	nx=normal[0];
	ny=normal[1];

	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
	num_gauss=3;
	GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);

	//recover thickness and bed at two grids
	h1=thickness_list[0];
	h2=thickness_list[1];
	b1=bed_list[0];
	b2=bed_list[1];

	//compute Jacobian of segment
	Jdet=1./2.*length;

	for(ig=0;ig<num_gauss;ig++){

		thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2;
		bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;

		if (fill==1){
			//icefront ends in water: 
			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
			air_pressure=0;

			//Now deal with water pressure
			surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
			base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
			water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
		}
		else if (fill==0){
			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
			air_pressure=0;
			water_pressure=0;
		}
		else{
			throw ErrorException(__FUNCT__," fill type not supported yet");
		}

		pressure = ice_pressure + water_pressure + air_pressure;

		pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx;
		pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny;
		pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx;
		pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny;

	} //for(ig=0;ig<num_gauss;ig++)

	xfree((void**)&segment_gauss_coord);
	xfree((void**)&gauss_weights);
}

#undef __FUNCT__ 
#define __FUNCT__ "Icefront::PenaltyCreateKMatrix"
void  Icefront::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){
	/*do nothing: */
}
		
#undef __FUNCT__ 
#define __FUNCT__ "Icefront::PenaltyCreatePVector"
void  Icefront::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){
	/*do nothing: */
}

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

