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

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

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

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

	this->hnodes=NULL;
	this->nodes= 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) {{{1*/
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_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=(int)(*(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=(int) *(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+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0);
		icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1);
	}
	else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
		icefront_node_ids[0]=iomodel->nodecounter+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0);
		icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1);
		icefront_node_ids[2]=iomodel->nodecounter+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2);
		icefront_node_ids[3]=iomodel->nodecounter+(int)*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3);
	}
	else _error_("in_icefront_type %s not supported yet!",EnumToStringx(in_icefront_type));

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

	/*Fill*/
	icefront_fill=(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->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(TypeEnum,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->element= NULL;
	this->matpar= NULL;
}


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

/*Object virtual functions definitions:*/
/*FUNCTION Icefront::Echo {{{1*/
void Icefront::Echo(void){
	printf("Icefront:\n");
	printf("   id: %i\n",id);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	hnodes->Echo();
	helement->Echo();
	hmatpar->Echo();
	printf("   parameters: %p\n",parameters);
	if(parameters)parameters->Echo();
	printf("   inputs: %p\n",inputs);
	if(inputs)inputs->Echo();
}
/*}}}*/
/*FUNCTION Icefront::DeepEcho{{{1*/
void Icefront::DeepEcho(void){

	printf("Icefront:\n");
	printf("   id: %i\n",id);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	hnodes->DeepEcho();
	helement->DeepEcho();
	hmatpar->DeepEcho();
	printf("   parameters: %p\n",parameters);
	if(parameters)parameters->DeepEcho();
	printf("   inputs: %p\n",inputs);
	if(inputs)inputs->DeepEcho();
}
/*}}}*/
/*FUNCTION Icefront::Id {{{1*/
int    Icefront::Id(void){ return id; }
/*}}}*/
/*FUNCTION Icefront::MyRank {{{1*/
int    Icefront::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
#ifdef _SERIAL_
/*FUNCTION Icefront::Marshall {{{1*/
void  Icefront::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;
	char* marshalled_inputs=NULL;
	int   marshalled_inputs_size;

	/*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,&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);
	helement->Marshall(&marshalled_dataset);
	hmatpar->Marshall(&marshalled_dataset);

	/*Marshall inputs: */
	marshalled_inputs_size=inputs->MarshallSize();
	marshalled_inputs=inputs->Marshall();
	memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
	marshalled_dataset+=marshalled_inputs_size;

	/*parameters: don't do anything about it. parameters are marshalled somewhere else!*/

	xfree((void**)&marshalled_inputs);

	*pmarshalled_dataset=marshalled_dataset;
}
/*}}}*/
/*FUNCTION Icefront::MarshallSize {{{1*/
int   Icefront::MarshallSize(){
	
	return sizeof(id)
		+sizeof(analysis_type)
		+hnodes->MarshallSize()
		+helement->MarshallSize()
		+hmatpar->MarshallSize()
		+inputs->MarshallSize()
		+sizeof(int); //sizeof(int) for enum type
}
/*}}}*/
/*FUNCTION Icefront::Demarshall {{{1*/
void  Icefront::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   i;

	/*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);
	helement=new Hook(); helement->Demarshall(&marshalled_dataset);
	hmatpar=new Hook(); hmatpar->Demarshall(&marshalled_dataset);
	
	/*pointers are garbabe, until configuration is carried out: */
	nodes=NULL;
	element=NULL;
	matpar=NULL;

	/*demarshall inputs: */
	inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 

	/*parameters: may not exist even yet, so let Configure handle it: */
	this->parameters=NULL;

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

	return IcefrontEnum;

}
/*}}}*/
/*FUNCTION Icefront::copy {{{1*/
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->helement=(Hook*)this->helement->copy();
	icefront->hmatpar=(Hook*)this->hmatpar->copy();

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

	return icefront;

}
/*}}}*/

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

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

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

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

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

	/*Checks in debugging mode*/
	/*{{{2*/
	_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 %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}

	/*Add to global Vector*/
	if(pe){
		pe->AddToGlobal(pf);
		delete pe;
	}
}
/*}}}*/
/*FUNCTION Icefront::CreateJacobianMatrix{{{1*/
void  Icefront::CreateJacobianMatrix(Mat Jff){
	this->CreateKMatrix(Jff,NULL);
}
/*}}}1*/
/*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/
void  Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){
	/*do nothing: */
	return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreatePVector{{{1*/
void  Icefront::PenaltyCreatePVector(Vec pf,double kmax){
	/*do nothing: */
	return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/
void  Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
}
/*}}}1*/
/*FUNCTION Icefront::InAnalysis{{{1*/
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(double* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVector(double* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(int* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVector(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(bool* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVector(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type) {{{1*/
void  Icefront::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(double* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVectorDakota(double* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{1*/
void  Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(double constant, int name) {{{1*/
void  Icefront::InputUpdateFromConstant(double constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(int constant, int name) {{{1*/
void  Icefront::InputUpdateFromConstant(int constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(bool constant, int name) {{{1*/
void  Icefront::InputUpdateFromConstant(bool constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromSolution{{{1*/
void  Icefront::InputUpdateFromSolution(double* solution){
	/*Nothing updated yet*/
}
/*}}}*/

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

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

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

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

	/*Intermediary*/
	int        ig,index1,index2,fill;
	double     Jdet;
	double     thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
	double     water_pressure,air_pressure,surface_under_water,base_under_water;
	double     xyz_list[numnodes][3];
	double     normal[2];
	double     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],nodes,numnodes);
	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 %s not supported yet",EnumToStringx(fill));
		}
		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 {{{1*/
ElementVector* Icefront::CreatePVectorAdjointHoriz(void){

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

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

	/*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(TypeEnum,MacAyeal2dIceFrontEnum));
	ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();

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

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

	/*Intermediaries*/
	int         i,j,ig,index1,index2,index3,index4;
	int         fill;
	double      surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
	double      water_pressure,air_pressure;
	double      Jdet,z_g;
	double      xyz_list[NUMVERTICESQUA][3];
	double      normal[3];
	double      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], nodes, 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: */
	double zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
	double 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 %s not supported yet",EnumToStringx(fill));
		}
		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{{{1*/
ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){

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

	/*Intermediaries*/
	int         i,j,ig,index1,index2,index3,index4;
	int         fill;
	double      pressure,rho_water,gravity;
	double      water_pressure,air_pressure;
	double      Jdet,z_g;
	double      xyz_list[NUMVERTICESQUA][3];
	double      normal[3];
	double      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], nodes, 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: */
	double zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
	double 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 %s not supported yet",EnumToStringx(fill));
		}
		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 {{{1*/
void  Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){

	int i,j;
	int numberofdofs=0;
	int count=0;
	int type;
	int numberofnodes=2;

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

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

	/*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(i=0;i<numberofnodes;i++){
		numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Allocate: */
	doflist=(int*)xmalloc(numberofdofs*sizeof(int));

	/*Populate: */
	count=0;
	for(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 {{{1*/
void Icefront:: GetSegmentNormal(double* normal,double xyz_list[4][3]){

	/*Build unit outward pointing vector*/
	const int numnodes=NUMVERTICESSEG;
	double vector[2];
	double 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 {{{1*/
void Icefront:: GetQuadNormal(double* normal,double xyz_list[4][3]){

	/*Build unit outward pointing vector*/
	double AB[3];
	double AC[3];
	double 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;
}
/*}}}*/
