Changeset 24486

Timestamp:

12/19/19 15:26:13 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: making enthalpy work with HO elements

File:

• issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (modified) (9 diffs)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -502 +502 @@
 
         /*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 +1089 @@
                         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 +1123 @@
                         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 +1147 @@
                         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 +1184 @@
                         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 +1607 @@
         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 +1614 @@
                         //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 +1636 @@
                         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 +1674 @@
         }
         else{
-                element->AddInput2(EnthalpyPicardEnum,values,P1DGEnum);
+                element->AddInput2(EnthalpyPicardEnum,values,finite_element);
         }