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

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

#include "../classes.h"
#include "shared/shared.h"
/*}}}*/	

/*Load macros*/
#define NUMVERTICESSEG 2
#define NUMVERTICESQUA 4

/*Icefront constructors and destructor*/
/*FUNCTION Icefront::Icefront() {{{*/
Icefront::Icefront(){

	this->inputs=NULL;
	this->parameters=NULL;

	this->hnodes=NULL;
	this->nodes= NULL;
	this->hvertices=NULL;
	this->vertices= NULL;
	this->helement=NULL;
	this->element= NULL;
	this->hmatpar=NULL;
	this->matpar= NULL;
}
/*}}}*/
/*FUNCTION Icefront::Icefront(int id, int i, IoModel* iomodel,int analysis_type) {{{*/
Icefront::Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int in_analysis_type){

	int segment_width;
	int element;
	int num_nodes; 
	int dim;
	int numberofelements;

	/*icefront constructor data: */
	int  icefront_eid;
	int  icefront_mparid;
	int  icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
	int  icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
	int  icefront_fill;

	/*find parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);

	/*First, retrieve element index and element type: */
	if (dim==2){
		segment_width=4;
	}
	else{
		segment_width=6;
	}
	_assert_(iomodel->Data(DiagnosticIcefrontEnum));
	element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)

	/*Build ids for hook constructors: */
	icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
	icefront_mparid=numberofelements+1; //matlab indexing

	if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
		icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
		icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
		icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
		icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
	}
	else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
		icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
		icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
		icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
		icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
		icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
		icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
		icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
		icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
	}
	else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");

	if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum)
	 num_nodes=4;
	else 
	 num_nodes=2;

	/*Fill*/
	icefront_fill=reCast<int>(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-1]);

	/*Ok, we have everything to build the object: */
	this->id=icefront_id;
	this->analysis_type=in_analysis_type;

	/*Hooks: */
	this->hnodes=new Hook(icefront_node_ids,num_nodes);
	this->hvertices=new Hook(icefront_vertex_ids,num_nodes);
	this->helement=new Hook(&icefront_eid,1);
	this->hmatpar=new Hook(&icefront_mparid,1);

	//intialize  and add as many inputs per element as requested: 
	this->inputs=new Inputs();
	this->inputs->AddInput(new IntInput(FillEnum,icefront_fill));
	this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type));

	//parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
	this->parameters = NULL;
	this->nodes      = NULL;
	this->vertices   = NULL;
	this->element    = NULL;
	this->matpar     = NULL;
}

/*}}}*/
/*FUNCTION Icefront::~Icefront() {{{*/
Icefront::~Icefront(){
	delete inputs;
	this->parameters=NULL;
	delete hnodes;
	delete hvertices;
	delete helement;
	delete hmatpar;
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION Icefront::Echo {{{*/
void Icefront::Echo(void){
	_printString_("Icefront:" << "\n");
	_printString_("   id: " << id << "\n");
	_printString_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
	hnodes->Echo();
	hvertices->Echo();
	helement->Echo();
	hmatpar->Echo();
	_printString_("   parameters: " << parameters << "\n");
	if(parameters)parameters->Echo();
	_printString_("   inputs: " << inputs << "\n");
	if(inputs)inputs->Echo();
}
/*}}}*/
/*FUNCTION Icefront::DeepEcho{{{*/
void Icefront::DeepEcho(void){

	_printString_("Icefront:" << "\n");
	_printString_("   id: " << id << "\n");
	_printString_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
	hnodes->DeepEcho();
	hvertices->DeepEcho();
	helement->DeepEcho();
	hmatpar->DeepEcho();
	_printString_("   parameters: " << parameters << "\n");
	if(parameters)parameters->DeepEcho();
	_printString_("   inputs: " << inputs << "\n");
	if(inputs)inputs->DeepEcho();
}
/*}}}*/
/*FUNCTION Icefront::Id {{{*/
int    Icefront::Id(void){ return id; }
/*}}}*/
/*FUNCTION Icefront::ObjectEnum{{{*/
int Icefront::ObjectEnum(void){

	return IcefrontEnum;

}
/*}}}*/
/*FUNCTION Icefront::copy {{{*/
Object* Icefront::copy() {

	Icefront* icefront=NULL;

	icefront=new Icefront();

	/*copy fields: */
	icefront->id=this->id;
	icefront->analysis_type=this->analysis_type;
	if(this->inputs){
		icefront->inputs=(Inputs*)this->inputs->Copy();
	}
	else{
		icefront->inputs=new Inputs();
	}
	/*point parameters: */
	icefront->parameters=this->parameters;

	/*now deal with hooks and objects: */
	icefront->hnodes    = (Hook*)this->hnodes->copy();
	icefront->hvertices = (Hook*)this->hvertices->copy();
	icefront->helement  = (Hook*)this->helement->copy();
	icefront->hmatpar   = (Hook*)this->hmatpar->copy();

	/*corresponding fields*/
	icefront->nodes    = (Node**)icefront->hnodes->deliverp();
	icefront->vertices = (Vertex**)icefront->hvertices->deliverp();
	icefront->element  = (Element*)icefront->helement->delivers();
	icefront->matpar   = (Matpar*)icefront->hmatpar->delivers();

	return icefront;

}
/*}}}*/

/*Load virtual functions definitions:*/
/*FUNCTION Icefront::Configure {{{*/
void  Icefront::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((DataSet*)nodesin);
	hvertices->configure((DataSet*)verticesin);
	helement->configure((DataSet*)elementsin);
	hmatpar->configure((DataSet*)materialsin);

	/*Initialize hooked fields*/
	this->nodes    = (Node**)hnodes->deliverp();
	this->vertices = (Vertex**)hvertices->deliverp();
	this->element  = (Element*)helement->delivers();
	this->matpar   = (Matpar*)hmatpar->delivers();

	/*point parameters to real dataset: */
	this->parameters=parametersin;
}
/*}}}*/
/*FUNCTION Icefront::SetCurrentConfiguration {{{*/
void  Icefront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
}
/*}}}*/
/*FUNCTION Icefront::CreateKMatrix {{{*/
void  Icefront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){

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

}
/*}}}*/
/*FUNCTION Icefront::CreatePVector {{{*/
void  Icefront::CreatePVector(Vector<IssmDouble>* pf){

	/*Checks in debugging mode*/
	/*{{{*/
	_assert_(nodes);
	_assert_(element);
	_assert_(matpar);
	/*}}}*/

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

	/*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
		#ifdef _HAVE_DIAGNOSTIC_
		case DiagnosticHorizAnalysisEnum:
			pe=CreatePVectorDiagnosticHoriz();
			break;
		#endif
		#ifdef _HAVE_CONTROL_
		case AdjointHorizAnalysisEnum:
			pe=CreatePVectorAdjointHoriz();
			break;
		#endif
		default:
			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
	}

	/*Add to global Vector*/
	if(pe){
		pe->AddToGlobal(pf);
		delete pe;
	}
}
/*}}}*/
/*FUNCTION Icefront::CreateJacobianMatrix{{{*/
void  Icefront::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
	this->CreateKMatrix(Jff,NULL);
}
/*}}}*/
/*FUNCTION Icefront::GetNodesSidList{{{*/
void Icefront::GetNodesSidList(int* sidlist){

	int type;
	inputs->GetInputValue(&type,IceFrontTypeEnum);
	_assert_(sidlist);
	_assert_(nodes);

	switch(type){
		case MacAyeal2dIceFrontEnum:
		case MacAyeal3dIceFrontEnum:
			for(int i=0;i<NUMVERTICESSEG;i++) sidlist[i]=nodes[i]->Sid();
			return;
#ifdef _HAVE_3D_
		case PattynIceFrontEnum:
		case StokesIceFrontEnum:
			for(int i=0;i<NUMVERTICESQUA;i++) sidlist[i]=nodes[i]->Sid();
			return;
#endif
		default:
			_error_("Icefront type " << EnumToStringx(type) << " not supported yet");
	}
}
/*}}}*/
/*FUNCTION Icefront::GetNumberOfNodes{{{*/
int Icefront::GetNumberOfNodes(void){

	int type;
	inputs->GetInputValue(&type,IceFrontTypeEnum);

	switch(type){
		case MacAyeal2dIceFrontEnum:
			return NUMVERTICESSEG;
#ifdef _HAVE_3D_
		case MacAyeal3dIceFrontEnum:
			return NUMVERTICESSEG;
		case PattynIceFrontEnum:
			return NUMVERTICESQUA;
		case StokesIceFrontEnum:
			return NUMVERTICESQUA;
#endif
		default:
			_error_("Icefront type " << EnumToStringx(type) << " not supported yet");
	}

}
/*}}}*/
/*FUNCTION Icefront::IsPenalty{{{*/
bool Icefront::IsPenalty(void){
	return false;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreateKMatrix {{{*/
void  Icefront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs, IssmDouble kmax){
	/*do nothing: */
	return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreatePVector{{{*/
void  Icefront::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){
	/*do nothing: */
	return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{*/
void  Icefront::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){
	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
}
/*}}}*/
/*FUNCTION Icefront::SetwiseNodeConnectivity{{{*/
void Icefront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){

	/*Output */
	int d_nz = 0;
	int o_nz = 0;

	/*Loop over all nodes*/
	for(int i=0;i<this->GetNumberOfNodes();i++){

		if(!flags[this->nodes[i]->Sid()]){

			/*flag current node so that no other element processes it*/
			flags[this->nodes[i]->Sid()]=true;

			/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
			switch(set2_enum){
				case FsetEnum:
					if(nodes[i]->indexing.fsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case GsetEnum:
					if(nodes[i]->indexing.gsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case SsetEnum:
					if(nodes[i]->indexing.ssize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				default: _error_("not supported");
			}
		}
	}

	/*Assign output pointers: */
	*pd_nz=d_nz;
	*po_nz=o_nz;
}
/*}}}*/
/*FUNCTION Icefront::InAnalysis{{{*/
bool Icefront::InAnalysis(int in_analysis_type){
	if (in_analysis_type==this->analysis_type)return true;
	else return false;
}
/*}}}*/

/*Update virtual functions definitions:*/
/*FUNCTION Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(int* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(bool* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type) {{{*/
void  Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(IssmDouble constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(int constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(int constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(bool constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(bool constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromSolution{{{*/
void  Icefront::InputUpdateFromSolution(IssmDouble* solution){
	/*Nothing updated yet*/
}
/*}}}*/

/*Icefront numerics: */
#ifdef _HAVE_DIAGNOSTIC_
/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{*/
ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){

	int type;
	inputs->GetInputValue(&type,IceFrontTypeEnum);

	switch(type){
		case MacAyeal2dIceFrontEnum:
			return CreatePVectorDiagnosticMacAyeal2d();
        #ifdef _HAVE_3D_
		case MacAyeal3dIceFrontEnum:
			return CreatePVectorDiagnosticMacAyeal3d();
		case PattynIceFrontEnum:
			return CreatePVectorDiagnosticPattyn();
		case StokesIceFrontEnum:
			return CreatePVectorDiagnosticStokes();
	    #endif
		default:
			_error_("Icefront type " << EnumToStringx(type) << " not supported yet");
	}
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){

	/*Constants*/
	const int numnodes= NUMVERTICESSEG;
	const int numdofs = numnodes *NDOF2;

	/*Intermediary*/
	int        ig,index1,index2,fill;
	IssmDouble     Jdet;
	IssmDouble     thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
	IssmDouble     water_pressure,air_pressure,surface_under_water,base_under_water;
	IssmDouble     xyz_list[numnodes][3];
	IssmDouble     normal[2];
	IssmDouble     L[2];
	GaussTria *gauss;

	Tria* tria=((Tria*)element);

	/*Initialize Element vector and return if necessary*/
	if(tria->IsOnWater()) return NULL;
	ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESSEG);
	Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
	Input* bed_input      =tria->inputs->GetInput(BedEnum);       _assert_(bed_input);
	inputs->GetInputValue(&fill,FillEnum);
	rho_water=matpar->GetRhoWater();
	rho_ice  =matpar->GetRhoIce();
	gravity  =matpar->GetG();
	GetSegmentNormal(&normal[0],xyz_list);

	/*Start looping on Gaussian points*/
	index1=tria->GetNodeIndex(nodes[0]);
	index2=tria->GetNodeIndex(nodes[1]);
	gauss=new GaussTria(index1,index2,3);

	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		thickness_input->GetInputValue(&thickness,gauss);
		bed_input->GetInputValue(&bed,gauss);

		switch(fill){
			case WaterEnum:
				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));
				break;
			case AirEnum:
				water_pressure=0;
				break;
			case IceEnum:
				water_pressure=-1.0/2.0*gravity*rho_ice*pow(thickness,2); // we are facing a wall of ice. use water_pressure to cancel the lithostatic pressure.
				break;
			default:
				_error_("fill type " << EnumToStringx(fill) << " not supported yet");
		}
		ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
		air_pressure=0;
		pressure = ice_pressure + water_pressure + air_pressure;

		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);

		for (int i=0;i<numnodes;i++){
			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i];
			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i];
		}
	}

	/*Transform load vector*/
	TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,XYEnum);

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
#endif

#ifdef _HAVE_CONTROL_
/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{*/
ElementVector* Icefront::CreatePVectorAdjointHoriz(void){

	/*No load vector applied to the adjoint*/
	return NULL;
}
/*}}}*/
#endif
#ifdef _HAVE_3D_
/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){

	Icefront *icefront = NULL;
	Penta    *penta    = NULL;
	Tria     *tria     = NULL;

	/*Cast element onto Penta*/
	penta   =(Penta*)this->element;

	/*Return if not on bed*/
	if(!penta->IsOnBed() || penta->IsOnWater()) return NULL;

	/*Spawn Tria and call MacAyeal2d*/
	tria    =(Tria*)penta->SpawnTria(0,1,2);
	icefront=(Icefront*)this->copy();
	icefront->element=tria;
	icefront->inputs->AddInput(new IntInput(IceFrontTypeEnum,MacAyeal2dIceFrontEnum));
	ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();

	/*clean-up and return*/
	delete tria->material;
	delete tria;
	delete icefront;
	return pe;
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){

	/*Constants*/
	const int numdofs = NUMVERTICESQUA *NDOF2;

	/*Intermediaries*/
	int         i,j,ig,index1,index2,index3,index4;
	int         fill;
	IssmDouble      surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
	IssmDouble      water_pressure,air_pressure;
	IssmDouble      Jdet,z_g;
	IssmDouble      xyz_list[NUMVERTICESQUA][3];
	IssmDouble      normal[3];
	IssmDouble      l1l4[4];
	GaussPenta *gauss = NULL;

	Penta* penta=(Penta*)element;

	/*Initialize Element vector and return if necessary*/
	if(penta->IsOnWater()) return NULL;
	ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,PattynApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
	Input* surface_input  =penta->inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
	inputs->GetInputValue(&fill,FillEnum);
	rho_water=matpar->GetRhoWater();
	rho_ice  =matpar->GetRhoIce();
	gravity  =matpar->GetG();
	GetQuadNormal(&normal[0],xyz_list);

	/*Identify which nodes are in the quad: */
	index1=element->GetNodeIndex(nodes[0]);
	index2=element->GetNodeIndex(nodes[1]);
	index3=element->GetNodeIndex(nodes[2]);
	index4=element->GetNodeIndex(nodes[3]);

	/* Start  looping on the number of gaussian points: */
	IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
	IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
	if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity
	else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
		penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
		z_g=penta->GetZcoord(gauss);
		surface_input->GetInputValue(&surface,gauss);

		switch(fill){
			case WaterEnum:
				water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
				break;
			case AirEnum:
				water_pressure=0;
				break;
			default:
				_error_("fill type " << EnumToStringx(fill) << " not supported yet");
		}
		ice_pressure=rho_ice*gravity*(surface-z_g);
		air_pressure=0;
		pressure = ice_pressure + water_pressure + air_pressure;

		for(i=0;i<NUMVERTICESQUA;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
	}

	/*Transform load vector*/
	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum);

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){

	/*Constants*/
	const int numdofs = NUMVERTICESQUA *NDOF4;

	/*Intermediaries*/
	int         i,j,ig,index1,index2,index3,index4;
	int         fill;
	IssmDouble      pressure,rho_water,gravity;
	IssmDouble      water_pressure,air_pressure;
	IssmDouble      Jdet,z_g;
	IssmDouble      xyz_list[NUMVERTICESQUA][3];
	IssmDouble      normal[3];
	IssmDouble      l1l4[4];
	GaussPenta *gauss = NULL;

	Penta* penta=(Penta*)element;

	/*Initialize Element vector and return if necessary*/
	if(penta->IsOnWater()) return NULL;
	ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,StokesApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
	inputs->GetInputValue(&fill,FillEnum);
	rho_water=matpar->GetRhoWater();
	gravity  =matpar->GetG();
	GetQuadNormal(&normal[0],xyz_list);

	/*Identify which nodes are in the quad: */
	index1=element->GetNodeIndex(nodes[0]);
	index2=element->GetNodeIndex(nodes[1]);
	index3=element->GetNodeIndex(nodes[2]);
	index4=element->GetNodeIndex(nodes[3]);

	/* Start  looping on the number of gaussian points: */
	IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
	IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
	if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity
	else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
		penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
		z_g=penta->GetZcoord(gauss);

		switch(fill){
			case WaterEnum:
				water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
				break;
			case AirEnum:
				water_pressure=0;
				break;
			default:
				_error_("fill type " << EnumToStringx(fill) << " not supported yet");
		}
		air_pressure=0;
		pressure = water_pressure + air_pressure; //no ice pressure fore Stokes

		for(i=0;i<NUMVERTICESQUA;i++){
			for(j=0;j<NDOF4;j++){
				if(j<3)  pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
				else     pe->values[i*NDOF4+j]+=0; //pressure term
			}
		}
	}

	/*Transform load vector*/
	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum);

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
#endif
/*FUNCTION Icefront::GetDofList {{{*/
void  Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){

	int numberofdofs=0;
	int count=0;
	int type;
	int numberofnodes=2;

	/*output: */
	int* doflist=NULL;

	/*recover type: */
	inputs->GetInputValue(&type,IceFrontTypeEnum);

	/*Some checks for debugging*/
	_assert_(nodes);

	/*How many nodes? :*/
	if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
	 numberofnodes=2;
	else 
	 numberofnodes=4;

	/*Figure out size of doflist: */
	for(int i=0;i<numberofnodes;i++){
		numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Allocate: */
	doflist=xNew<int>(numberofdofs);

	/*Populate: */
	count=0;
	for(int i=0;i<numberofnodes;i++){
		nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
/*FUNCTION Icefront::GetSegmentNormal {{{*/
void Icefront:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){

	/*Build unit outward pointing vector*/
	const int numnodes=NUMVERTICESSEG;
	IssmDouble vector[2];
	IssmDouble norm;

	vector[0]=xyz_list[1][0] - xyz_list[0][0];
	vector[1]=xyz_list[1][1] - xyz_list[0][1];

	norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));

	normal[0]= + vector[1]/norm;
	normal[1]= - vector[0]/norm;
}
/*}}}*/
/*FUNCTION Icefront::GetQuadNormal {{{*/
void Icefront:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){

	/*Build unit outward pointing vector*/
	IssmDouble AB[3];
	IssmDouble AC[3];
	IssmDouble norm;

	AB[0]=xyz_list[1][0] - xyz_list[0][0];
	AB[1]=xyz_list[1][1] - xyz_list[0][1];
	AB[2]=xyz_list[1][2] - xyz_list[0][2];
	AC[0]=xyz_list[2][0] - xyz_list[0][0];
	AC[1]=xyz_list[2][1] - xyz_list[0][1];
	AC[2]=xyz_list[2][2] - xyz_list[0][2];

	cross(normal,AB,AC);
	norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));

	for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
}
/*}}}*/