#include "./HydrologyDCEfficientAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Model processing*/
int  HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
	return 1;
}/*}}}*/
void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

	int         hydrology_model;
	bool        isefficientlayer;

	/*retrieve some parameters: */
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);

	/*Now, do we really want DC?*/
	if(hydrology_model!=HydrologydcEnum) return;

	/*Do we want an efficient layer*/
	iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
	if(!isefficientlayer) return;

	/*Nothing for now*/
}/*}}}*/
void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

	bool   isefficientlayer;
	int    hydrology_model;

	/*Now, do we really want DC?*/
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
	if(hydrology_model!=HydrologydcEnum) return;

	/*Do we want an efficient layer*/
	iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
	if(!isefficientlayer) return;

	/*Update elements: */
	int counter=0;
	for(int i=0;i<iomodel->numberofelements;i++){
		if(iomodel->my_elements[i]){
			Element* element=(Element*)elements->GetObjectByOffset(counter);
			element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
			counter++;
		}
	}
	iomodel->FetchDataToInput(elements,ThicknessEnum);
	iomodel->FetchDataToInput(elements,BaseEnum);
	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
	iomodel->FetchDataToInput(elements,EplHeadEnum);
	iomodel->FetchDataToInput(elements,SedimentHeadEnum);
	iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum);
	iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
	iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
	if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
	
	//	elements->InputDuplicate(HydrologydcEplInitialThicknessEnum,HydrologydcEplThicknessEnum);

}/*}}}*/
void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/

	/*Now, do we really want DC?*/
	int  hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
	if(hydrology_model!=HydrologydcEnum) return;

	/*Do we want an efficient layer*/
	bool isefficientlayer;
	iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
	if(!isefficientlayer) return;

	if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
	::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum);
	iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
}/*}}}*/
void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/

	/*Do we really want DC?*/
	int  hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
	if(hydrology_model!=HydrologydcEnum) return;

	/*Do we want an efficient layer*/
	bool isefficientlayer;
	iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
	if(!isefficientlayer) return;

	IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);

}/*}}}*/
void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
	/*Nothing for now*/
}/*}}}*/

/*Finite Element Analysis*/
void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/
	_error_("not implemented");
}/*}}}*/
ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
	/*Default, return NULL*/
	return NULL;
}/*}}}*/
ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
_error_("Not implemented");
}/*}}}*/
ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/

	/*Intermediaries*/
	bool     active_element;
	int      meshtype;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&meshtype,MeshTypeEnum);
	switch(meshtype){
		case Mesh2DhorizontalEnum:
			basalelement = element;
			break;
		case Mesh3DEnum:
			if(!element->IsOnBed()) return NULL;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
	}

	Input* active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
	active_element_input->GetInputValue(&active_element);

	/*Check that all nodes are active, else return empty matrix*/
	if(!active_element) {
	if(meshtype!=Mesh2DhorizontalEnum){
			basalelement->DeleteMaterials(); 
			delete basalelement;
		}
		return NULL;
	}
	/* Intermediaries */
	IssmDouble  D_scalar,Jdet,dt;
	IssmDouble  epl_thickness;
	IssmDouble  transfer;
	IssmDouble *xyz_list = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = basalelement->GetNumberOfNodes();

	/*Initialize Element vector*/
	ElementMatrix* Ke     = basalelement->NewElementMatrix();
	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
	IssmDouble     D[2][2]={0.};

	/*Retrieve all inputs and parameters*/
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
	Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum);          _assert_(thickness_input);
	Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
	Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
	Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
	Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);

	IssmDouble epl_specificstoring   = EplSpecificStoring(basalelement);
	IssmDouble epl_conductivity      = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss           ->GaussPoint(ig);
		basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
		thickness_input ->GetInputValue(&epl_thickness,gauss);

		/*Diffusivity*/
		D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;
		if(dt!=0.) D_scalar=D_scalar*dt;
		D[0][0]=D_scalar;
		D[1][1]=D_scalar;
		GetB(B,basalelement,xyz_list,gauss); 
		TripleMultiply(B,2,numnodes,1,
					&D[0][0],2,2,0,
					B,2,numnodes,0,
					&Ke->values[0],1);

		/*Transient*/
		if(dt!=0.){
			basalelement->NodalFunctions(&basis[0],gauss);
			D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;
			TripleMultiply(basis,numnodes,1,0,
						&D_scalar,1,1,0,
						basis,1,numnodes,0,
						&Ke->values[0],1);
			
			/*Transfer EPL part*/
			transfer=GetHydrologyKMatrixTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input);
			D_scalar=transfer*gauss->weight*Jdet*dt;
			TripleMultiply(basis,numnodes,1,0,
										 &D_scalar,1,1,0,
										 basis,1,numnodes,0,
										 &Ke->values[0],1);
		}
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	xDelete<IssmDouble>(B);
	delete gauss;
	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	return Ke;

}/*}}}*/
ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/

	/*Intermediaries*/
	bool     active_element;
	int      meshtype;
	Element* basalelement;

	/*Get basal element*/
	element->FindParam(&meshtype,MeshTypeEnum);
	switch(meshtype){
		case Mesh2DhorizontalEnum:
			basalelement = element;
			break;
		case Mesh3DEnum:
			if(!element->IsOnBed()) return NULL;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
	}

	Input* active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
	active_element_input->GetInputValue(&active_element);

	/*Check that all nodes are active, else return empty matrix*/
	if(!active_element) {
	if(meshtype!=Mesh2DhorizontalEnum){
			basalelement->DeleteMaterials(); 
			delete basalelement;
		}
		return NULL;
	}
	/*Intermediaries */
	IssmDouble dt,scalar,water_head;
	IssmDouble transfer;
	IssmDouble epl_thickness;
	IssmDouble Jdet;
	IssmDouble residual,connectivity;

	IssmDouble *xyz_list     = NULL;
	Input*      old_wh_input = NULL;

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes    = basalelement->GetNumberOfNodes();
	int numvertices = basalelement->GetNumberOfVertices();

	/*Initialize Element vector*/
	ElementVector* pe    = basalelement->NewElementVector();
	IssmDouble*    basis = xNew<IssmDouble>(numnodes);

	/*Retrieve all inputs and parameters*/
	basalelement->GetVerticesCoordinates(&xyz_list);
	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);	

	Input* thickness_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
	Input* sed_head_input  = basalelement->GetInput(SedimentHeadEnum);
	Input* epl_head_input  = basalelement->GetInput(EplHeadEnum);
	Input* sed_trans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
	Input* residual_input  = basalelement->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
	if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum);            _assert_(old_wh_input);}

	IssmDouble epl_specificstoring = EplSpecificStoring(basalelement);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);

		basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
		basalelement ->NodalFunctions(basis,gauss);

		/*Transient and transfer terms*/
		if(dt!=0.){
			old_wh_input    ->GetInputValue(&water_head,gauss);
			thickness_input ->GetInputValue(&epl_thickness,gauss);
			
			/*Dealing with the sediment part of the transfer term*/
			transfer=GetHydrologyPVectorTransfer(basalelement,gauss,thickness_input,sed_head_input,epl_head_input,sed_trans_input,residual_input);
			scalar = Jdet*gauss->weight*((water_head*epl_specificstoring*epl_thickness)+(transfer*dt));
			for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
		}
	}
	delete gauss;

	/*	Add residual if necessary*/
	gauss=basalelement->NewGauss();
	for(int iv=0;iv<numvertices;iv++){
		gauss->GaussVertex(iv);

		connectivity = IssmDouble(basalelement->VertexConnectivity(iv));
		residual_input->GetInputValue(&residual,gauss);
		pe->values[iv]+=residual/connectivity;

	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(basis);
	delete gauss;
	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	return pe;
}/*}}}*/
void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
	element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
}/*}}}*/
void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

	int meshtype,i;
	Element*   basalelement=NULL;

	element->FindParam(&meshtype,MeshTypeEnum);

	if(meshtype!=Mesh2DhorizontalEnum){
		if(!element->IsOnBed()) return;
		basalelement=element->SpawnBasalElement();
	}
	else{
		basalelement = element;
	}
	
	/*Intermediary*/
	int* doflist = NULL;

	/*Fetch number of nodes for this finite element*/
	int numnodes = basalelement->GetNumberOfNodes();

	/*Fetch dof list and allocate solution vector*/
	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
	IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numnodes;i++){
		eplHeads[i]=solution[doflist[i]];
		if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
	}
	/*Add input to the element: */
	element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum);

	/*Free ressources:*/
	xDelete<IssmDouble>(eplHeads);
	xDelete<int>(doflist);
	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
} /*}}}*/
void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
	/*Default, do nothing*/
	return;
}/*}}}*/

/*Intermediaries*/
IssmDouble HydrologyDCEfficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/
	IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
	IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
	IssmDouble epl_porosity          = element->GetMaterialParameter(HydrologydcEplPorosityEnum);
	IssmDouble epl_compressibility   = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum);
	IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
	return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));		 
}/*}}}*/
IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
	IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
	IssmDouble g                        = element->GetMaterialParameter(ConstantsGEnum);
	IssmDouble sediment_porosity        = element->GetMaterialParameter(HydrologydcSedimentPorosityEnum);
	IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
	IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum);
	IssmDouble water_compressibility    = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
	return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));		 
}/*}}}*/
IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/
	
	int transfermethod;
	IssmDouble epl_thickness;
	IssmDouble epl_head,sed_head;
	IssmDouble sediment_transmitivity;
	IssmDouble leakage,residual,transfer;

	IssmDouble sediment_thickness  = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
	IssmDouble sediment_storing    = SedimentStoring(element);
 	IssmDouble epl_specificstoring = EplSpecificStoring(element);		

	element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
	/*Switch between the different transfer methods cases*/
	switch(transfermethod){
	case 0:
		/*Just keepping the transfer to zero*/
		transfer=0.0;
		break;
	case 1:
		_assert_(epl_thick_input); 
		_assert_(sed_head_input); 
		_assert_(epl_head_input); 
		_assert_(sed_trans_input); 
		_assert_(residual_input); 
		/* get input */
		epl_thick_input->GetInputValue(&epl_thickness,gauss);
		sed_head_input->GetInputValue(&sed_head,gauss);
		epl_head_input->GetInputValue(&epl_head,gauss);
		sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
		residual_input->GetInputValue(&residual,gauss);
		element->FindParam(&leakage,HydrologydcLeakageFactorEnum);

		if(epl_head>sed_head){
			if(residual>0.0){	
				transfer=0.0;
			}
			else{
				transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity)/(sediment_thickness*leakage);
			}
		}
		else{
			transfer=(sediment_storing*sediment_transmitivity)/(sediment_thickness*leakage);
		}
		break;
	default:
		_error_("no case higher than 1 for the Transfer method");
	}
	
	return transfer;
}/*}}}*/
IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_thick_input, Input* sed_head_input, Input* epl_head_input, Input* sed_trans_input, Input* residual_input){/*{{{*/

	int transfermethod;
	IssmDouble epl_thickness;
	IssmDouble epl_head,sediment_head;
	IssmDouble sediment_transmitivity;
	IssmDouble leakage,residual,transfer;

	IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
	IssmDouble sediment_storing   = SedimentStoring(element);
 	IssmDouble epl_specificstoring = EplSpecificStoring(element);		

	element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
	/*Switch between the different transfer methods cases*/
	switch(transfermethod){
	case 0:
		/*Just keepping the transfer to zero*/
		transfer=0.0;
		break;
	case 1:
		_assert_(epl_thick_input); 
		_assert_(sed_head_input); 
		_assert_(epl_head_input); 
		_assert_(sed_trans_input); 
		_assert_(residual_input); 
		/* get input */
		epl_thick_input->GetInputValue(&epl_thickness,gauss);
		sed_head_input->GetInputValue(&sediment_head,gauss);
		epl_head_input->GetInputValue(&epl_head,gauss);
		sed_trans_input->GetInputValue(&sediment_transmitivity,gauss);
		residual_input->GetInputValue(&residual,gauss);
		element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
		
		if(epl_head>sediment_head){
			if(residual>0.0){	
				transfer=0.0;
			}
			else{
				transfer=(epl_specificstoring*epl_thickness*sediment_transmitivity*sediment_head)/(sediment_thickness*leakage);
			}
		}
		else{
			transfer=(sediment_storing*sediment_transmitivity*sediment_head)/(sediment_thickness*leakage);
		}
		break;
	default:
		_error_("no case higher than 1 for the Transfer method");
	}
	return transfer;
}/*}}}*/

void HydrologyDCEfficientAnalysis::ComputeEPLThickness(FemModel* femmodel){/*{{{*/

	bool        active_element;
	int         meshtype;
	IssmDouble  dt,A,B;
	IssmDouble  EPLgrad2;
	IssmDouble  EPL_N;

	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);

	for(int j=0;j<femmodel->elements->Size();j++){
		
		Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
		
		switch(meshtype){
		case Mesh2DhorizontalEnum:
			if(!element->IsOnBed()) return;			
			B = element->GetMaterialParameter(MaterialsRheologyBbarEnum);
			break;
		case Mesh3DEnum:
			B = element->GetMaterialParameter(MaterialsRheologyBEnum);
			break;
		default:
		_error_("not Implemented Yet");
		}
			
		int         numnodes      = element->GetNumberOfNodes();
		IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
		IssmDouble* eplhead       = xNew<IssmDouble>(numnodes);
		IssmDouble* epl_slopeX    = xNew<IssmDouble>(numnodes);
		IssmDouble* epl_slopeY    = xNew<IssmDouble>(numnodes);
		IssmDouble* old_thickness = xNew<IssmDouble>(numnodes);
		IssmDouble* ice_thickness = xNew<IssmDouble>(numnodes);
		IssmDouble* bed           = xNew<IssmDouble>(numnodes);

		Input* 	active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);		
		active_element_input->GetInputValue(&active_element);
		element->FindParam(&dt,TimesteppingTimeStepEnum);
	
		/*For now, assuming just one way to compute EPL thickness*/
		IssmDouble gravity          = element->GetMaterialParameter(ConstantsGEnum);
		IssmDouble rho_water        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
		IssmDouble rho_ice          = element->GetMaterialParameter(MaterialsRhoIceEnum);
		IssmDouble n                =	element->GetMaterialParameter(MaterialsRheologyNEnum);
		IssmDouble latentheat       = element->GetMaterialParameter(MaterialsLatentheatEnum);
		IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
		IssmDouble init_thick       =	element->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
		
		A=pow(B,-n);
		
		element->GetInputListOnVertices(&eplhead[0],EplHeadEnum);
		element->GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); 
		element->GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
		element->GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
		element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
		element->GetInputListOnVertices(&bed[0],BaseEnum);
			
		if(!active_element){
			
			/*Keeping thickness to initial value if EPL is not active*/
			for(int i=0;i<numnodes;i++){
				thickness[i]=init_thick;
			}
		}
		else{
			for(int i=0;i<numnodes;i++){
				
				/*Compute first the effective pressure in the EPL*/
				EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
				if(EPL_N<0.0)EPL_N=0.0;
				/*Get then the square of the gradient of EPL heads*/
				EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i]+epl_slopeY[i]*epl_slopeY[i]);
				
				/*And proceed to the real thing*/
				thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
				
				/*Take care of otherthikening*/
				if(thickness[i]>10.0*init_thick){
					thickness[i] = 10.0*init_thick;
				}
			}
		}
		element->AddInput(HydrologydcEplThicknessEnum,thickness,P1Enum);
		xDelete<IssmDouble>(thickness);
		xDelete<IssmDouble>(eplhead);
		xDelete<IssmDouble>(epl_slopeX);
		xDelete<IssmDouble>(epl_slopeY);
		xDelete<IssmDouble>(old_thickness);
		xDelete<IssmDouble>(ice_thickness);
		xDelete<IssmDouble>(bed);
	}
}
/*}}}*/
void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dN/dx ]
	 *          [ dN/dy ]
	 * where N is the finiteelement function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

	/*Fetch number of nodes for this finite element*/
	int numnodes = element->GetNumberOfNodes();

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B: */
	for(int i=0;i<numnodes;i++){
		B[numnodes*0+i] = dbasis[0*numnodes+i];
		B[numnodes*1+i] = dbasis[1*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}/*}}}*/
void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element){

	bool        active_element;
	int         i,j;
	int         meshtype;
	IssmDouble  h_max;
	IssmDouble  sedheadmin;
	Element*   basalelement=NULL;

	/*Get basal element*/
	element->FindParam(&meshtype,MeshTypeEnum);
	switch(meshtype){
		case Mesh2DhorizontalEnum:
			basalelement = element;
			break;
		case Mesh3DEnum:
			if(!element->IsOnBed()) return;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
	}

	/*Intermediaries*/

	int         numnodes      =basalelement->GetNumberOfNodes();
	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
	IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
	IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
	IssmDouble* eplhead       =xNew<IssmDouble>(numnodes);
	IssmDouble* residual      =xNew<IssmDouble>(numnodes);

	IssmDouble init_thick =	basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);

	Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
	active_element_input->GetInputValue(&active_element);

	basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);	
	basalelement-> GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);	
	basalelement-> GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
	basalelement-> GetInputListOnVertices(&eplhead[0],EplHeadEnum);
	basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);

	/*Get minimum sediment head of the element*/
	sedheadmin=sedhead[0];
	for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];

	for(i=0;i<numnodes;i++){
		/*Activate EPL if residual is >0 */
		if(residual[i]>0.){
			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
		}

		/*If mask was already one, keep one*/
		else if(old_active[i]>0.){
			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
			/*If epl thickness gets under , close the layer*/
			if(epl_thickness[i]<0.001*init_thick){
				vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
				epl_thickness[i]=init_thick;
			}
		}
		/*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
		GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
		if(eplhead[i]>=h_max && active_element){
			for(j=0;j<numnodes;j++){
				if(old_active[j]>0.){
					vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
				}
				/*Increase of the domain is on the downstream node in term of sediment head*/
				if(sedhead[j] == sedheadmin){
					vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
				}
			}
		}
	}
	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	xDelete<IssmDouble>(epl_thickness);
	xDelete<IssmDouble>(old_active);
	xDelete<IssmDouble>(sedhead);
	xDelete<IssmDouble>(eplhead);
	xDelete<IssmDouble>(residual);
}
/*}}}*/
void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
	/*Constants*/

	int      meshtype;
	Element*   basalelement=NULL;

	/*Get basal element*/
	element->FindParam(&meshtype,MeshTypeEnum);
	switch(meshtype){
		case Mesh2DhorizontalEnum:
			basalelement = element;
			break;
		case Mesh3DEnum:
			if(!element->IsOnBed()) return;
			basalelement = element->SpawnBasalElement();
			break;
		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
	}
	
	const int   numnodes = basalelement->GetNumberOfNodes();
	IssmDouble  flag     = 0.;
	IssmDouble* active   = xNew<IssmDouble>(numnodes);
		
	basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
	
	for(int i=0;i<numnodes;i++) flag+=active[i];

	if(flag>0.){
		for(int i=0;i<numnodes;i++){
			active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
		}
	}
	else{
		/*Do not do anything: at least one node is active for this element but this element is not solved for*/
	}
	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
	xDelete<IssmDouble>(active);
}

void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
	
	int        hmax_flag;
	IssmDouble h_max;
	IssmDouble rho_ice,rho_water;
	IssmDouble thickness,bed;
	/*Get the flag to the limitation method*/
	element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
	
	/*Switch between the different cases*/
	switch(hmax_flag){
	case 0:
		h_max=1.0e+10;
		break;
	case 1:
		element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
		break;
	case 2:
		/*Compute max*/
		rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
		rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
		element-> GetInputValue(&thickness,innode,ThicknessEnum);
		element-> GetInputValue(&bed,innode,BaseEnum);
		h_max=((rho_ice*thickness)/rho_water)+bed;
		break;
	case 3:
		_error_("Using normal stress  not supported yet");
		break;
	default:
		_error_("no case higher than 3 for SedimentlimitFlag");
	}
	/*Assign output pointer*/
	*ph_max=h_max;
}
/*}}}*/
