Index: issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 24344)
+++ issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 24385)
@@ -383,4 +383,31 @@
 	return pe;
 }/*}}}*/
+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::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum);
@@ -437,4 +464,84 @@
 
 /*Intermediaries*/
+IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
+	IssmDouble epl_storing;
+	IssmDouble water_sheet,storing;
+	IssmDouble epl_thickness,prestep_head,base_elev;
+	IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
+	IssmDouble g                     = element->FindParam(ConstantsGEnum);
+	IssmDouble epl_porosity					 = element->FindParam(HydrologydcEplPorosityEnum);
+	IssmDouble epl_compressibility	 = element->FindParam(HydrologydcEplCompressibilityEnum);
+	IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
+
+	epl_thick_input->GetInputValue(&epl_thickness,gauss);
+	epl_head_input->GetInputValue(&prestep_head,gauss);
+	base_input->GetInputValue(&base_elev,gauss);
+	water_sheet=max(0.0,(prestep_head-base_elev));
+	storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
+
+	/* //porosity for unconfined region */
+	/* if (water_sheet<=0.9*epl_thickness){ */
+	/* 	epl_storing=epl_porosity; */
+	/* } */
+	/* //continuity ramp */
+	/* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */
+	/* 	epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */
+	/* } */
+	/* //storing coefficient for confined */
+	/* else{ */
+	/* 	epl_storing=storing; */
+	/* } */
+ 	/* return epl_storing; */
+	return storing;
+}/*}}}*/
+IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
+	IssmDouble epl_transmitivity;
+	IssmDouble water_sheet;
+	IssmDouble epl_thickness,base_elev,prestep_head;
+	IssmDouble epl_conductivity      = element->FindParam(HydrologydcEplConductivityEnum);
+	epl_thick_input->GetInputValue(&epl_thickness,gauss);
+	epl_head_input->GetInputValue(&prestep_head,gauss);
+	base_input->GetInputValue(&base_elev,gauss);
+
+	water_sheet=max(0.0,(prestep_head-base_elev));
+	epl_transmitivity=epl_conductivity*epl_thickness;
+	//epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
+	return epl_transmitivity;
+}/*}}}*/
+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->FindParam(MaterialsRhoFreshwaterEnum);
+		rho_ice   = element->FindParam(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;
+}
+/*}}}*/
 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/
 
@@ -582,31 +689,4 @@
 	}
 }/*}}}*/
-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, Vector<IssmDouble>* recurence, Element* element){/*{{{*/
 
@@ -703,49 +783,4 @@
 }
 /*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
-	IssmDouble epl_storing;
-	IssmDouble water_sheet,storing;
-	IssmDouble epl_thickness,prestep_head,base_elev;
-	IssmDouble rho_freshwater        = element->FindParam(MaterialsRhoFreshwaterEnum);
-	IssmDouble g                     = element->FindParam(ConstantsGEnum);
-	IssmDouble epl_porosity					 = element->FindParam(HydrologydcEplPorosityEnum);
-	IssmDouble epl_compressibility	 = element->FindParam(HydrologydcEplCompressibilityEnum);
-	IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
-
-	epl_thick_input->GetInputValue(&epl_thickness,gauss);
-	epl_head_input->GetInputValue(&prestep_head,gauss);
-	base_input->GetInputValue(&base_elev,gauss);
-	water_sheet=max(0.0,(prestep_head-base_elev));
-	storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
-
-	/* //porosity for unconfined region */
-	/* if (water_sheet<=0.9*epl_thickness){ */
-	/* 	epl_storing=epl_porosity; */
-	/* } */
-	/* //continuity ramp */
-	/* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */
-	/* 	epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */
-	/* } */
-	/* //storing coefficient for confined */
-	/* else{ */
-	/* 	epl_storing=storing; */
-	/* } */
- 	/* return epl_storing; */
-	return storing;
-}/*}}}*/
-IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/
-	IssmDouble epl_transmitivity;
-	IssmDouble water_sheet;
-	IssmDouble epl_thickness,base_elev,prestep_head;
-	IssmDouble epl_conductivity      = element->FindParam(HydrologydcEplConductivityEnum);
-	epl_thick_input->GetInputValue(&epl_thickness,gauss);
-	epl_head_input->GetInputValue(&prestep_head,gauss);
-	base_input->GetInputValue(&base_elev,gauss);
-
-	water_sheet=max(0.0,(prestep_head-base_elev));
-	epl_transmitivity=epl_conductivity*epl_thickness;
-	//epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));
-	return epl_transmitivity;
-}/*}}}*/
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/
 	/*Constants*/
@@ -796,37 +831,2 @@
 }
 /*}}}*/
-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->FindParam(MaterialsRhoFreshwaterEnum);
-		rho_ice   = element->FindParam(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;
-}
-/*}}}*/
