Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 24485)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 24486)
@@ -502,9 +502,10 @@
 
 	/*feed updated variables back into model*/
+	int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
 	if(dt!=0.){
-		element->AddInput2(enthalpy_enum,enthalpies,P1DGEnum);
-		element->AddInput2(WatercolumnEnum,watercolumns,P1DGEnum);
-	}
-	element->AddInput2(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1DGEnum);
+		element->AddInput2(enthalpy_enum,enthalpies,finite_element);
+		element->AddInput2(WatercolumnEnum,watercolumns,finite_element);
+	}
+	element->AddInput2(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,finite_element);
 
 	/*Clean up and return*/
@@ -1088,5 +1089,6 @@
 			drainage[k]=DrainageFunctionWaterfraction(waterfractions[k], dt);
 		}
-		element->AddInput2(WaterfractionDrainageEnum,drainage,P1DGEnum);
+		int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
+		element->AddInput2(WaterfractionDrainageEnum,drainage,finite_element);
 
 		xDelete<IssmDouble>(waterfractions);
@@ -1121,5 +1123,6 @@
 			drainage_int[k]*=thicknesses[k];
 		}
-		element->AddInput2(WaterfractionDrainageIntegratedEnum, drainage_int, P1DGEnum);
+		int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
+		element->AddInput2(WaterfractionDrainageIntegratedEnum, drainage_int,finite_element);
 
 		xDelete<IssmDouble>(drainage_int);
@@ -1144,5 +1147,6 @@
 			watercolumn[basalnodeindices[k]]+=dt*drainage_int[basalnodeindices[k]];
 		}
-		element->AddInput2(WatercolumnEnum, watercolumn,P1DGEnum);
+		int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
+		element->AddInput2(WatercolumnEnum, watercolumn,finite_element);
 
 		xDelete<IssmDouble>(watercolumn);
@@ -1180,6 +1184,7 @@
 			element->ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
 		}
-		element->AddInput2(WaterfractionEnum,waterfractions,P1DGEnum);
-		element->AddInput2(EnthalpyEnum,enthalpies,P1DGEnum);
+		int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
+		element->AddInput2(WaterfractionEnum,waterfractions,finite_element);
+		element->AddInput2(EnthalpyEnum,enthalpies,finite_element);
 
 		xDelete<IssmDouble>(enthalpies);
@@ -1602,4 +1607,5 @@
 	element->GetInputValue(&converged,ConvergedEnum);
 	element->GetInputListOnNodes(&pressure[0],PressureEnum);
+	int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
 	if(converged){
 		for(i=0;i<numnodes;i++){
@@ -1608,8 +1614,7 @@
 			//if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
 		}
-		_assert_(element->GetElementType()==P1Enum);
-		element->AddInput2(EnthalpyEnum,values,P1DGEnum);
-		element->AddInput2(WaterfractionEnum,waterfraction,P1DGEnum);
-		element->AddInput2(TemperatureEnum,temperature,P1DGEnum);
+		element->AddInput2(EnthalpyEnum,values,finite_element);
+		element->AddInput2(WaterfractionEnum,waterfraction,finite_element);
+		element->AddInput2(TemperatureEnum,temperature,finite_element);
 
 		IssmDouble* n = xNew<IssmDouble>(numnodes);
@@ -1631,35 +1636,35 @@
 			case BuddJackaEnum:
 				for(i=0;i<numnodes;i++) B[i]=BuddJacka(temperature[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case CuffeyEnum:
 				for(i=0;i<numnodes;i++) B[i]=Cuffey(temperature[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case CuffeyTemperateEnum:
 				for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i],n[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case PatersonEnum:
 				for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case NyeH2OEnum:
 				for(i=0;i<numnodes;i++) B[i]=NyeH2O(values[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case NyeCO2Enum:
 				for(i=0;i<numnodes;i++) B[i]=NyeCO2(values[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 			case ArrheniusEnum:{
 				element->GetVerticesCoordinates(&xyz_list);
 				for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n[i]);
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum);
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element);
 				break;
 				}
 			case LliboutryDuvalEnum:{
 				for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],n[i],element->FindParam(MaterialsBetaEnum),element->FindParam(ConstantsReferencetemperatureEnum),element->FindParam(MaterialsHeatcapacityEnum),element->FindParam(MaterialsLatentheatEnum)); 
-				element->AddInput2(MaterialsRheologyBEnum,&B[0],P1DGEnum); 
+				element->AddInput2(MaterialsRheologyBEnum,&B[0],finite_element); 
 				break; 
 				}
@@ -1669,5 +1674,5 @@
 	}
 	else{
-		element->AddInput2(EnthalpyPicardEnum,values,P1DGEnum);
+		element->AddInput2(EnthalpyPicardEnum,values,finite_element);
 	}