Changeset 11425


Ignore:
Timestamp:
02/14/12 09:27:32 (13 years ago)
Author:
cborstad
Message:

reverted to revision 11411: working to resolve some missed changes during initial branch

Location:
issm/trunk-jpl-damage
Files:
28 deleted
83 edited
5 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl-damage/INSTALL

    r11414 r11425  
    1 For the installation process please go to the ISSM website:
    2 
     1For the installation process please go to the ISSM website
    32http://issm.jpl.nasa.gov/

  • issm/trunk-jpl-damage/cron/cronfiles/linux_cronfile

    r11416 r11425  
    66
    77#cronjob
    8 00 09 * * 1-5 cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_daily
    9 30 12 * * 1-5 cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_daily
    10 00 15 * * 1-5 cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_daily
    11 00 18 * * 1-5 cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_daily
    12 00 23 * * 1-5 cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_nightly
    13 00 20 * * 6   cd /u/astrid-r1b/schlegel/ExecutionNightlyRun && rm -r test*
    14 00 21 * * 6   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/trunk-jpl/test/NightlyRun && make clean
    15 00 22 * * 6   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/trunk-jpl/test/NightlyRun && rm -r qmu*
    16 00 23 * * 6   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_validation
    17 00 21 * * 7   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/trunk/test/NightlyRun && make clean
    18 00 22 * * 7   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/trunk/test/NightlyRun && rm -r qmu*
    19 00 23 * * 7   cd /u/astrid-r1b/schlegel/issmuci/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_schlegel_ucitrunk
     800 09 * * 1-5 cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_daily
     930 12 * * 1-5 cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_daily
     1000 15 * * 1-5 cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_daily
     1100 18 * * 1-5 cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_daily
     1200 23 * * 1-5 cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_nightly
     1300 23 * * 6   cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_validation
     1400 23 * * 7   cd /u/astrid-r1b/seroussi/issm/trunk-jpl/cron/ && ./nightlyrun.sh configs/linux64_seroussi_ucitrunk

  • issm/trunk-jpl-damage/externalpackages/matlab/install.sh

    r11413 r11425  
    55
    66#Select or create a new simlink
    7 #ln -s /usr/local/pkgs/matlab-7.6/ install
     7ln -s /usr/local/pkgs/matlab-7.6/ install
    88#ln -s /usr/local/matlab704/ install
    99#ln -s /usr/local/matlab711/ install
    10 ln -s /usr/local/matlab712/ install
     10#ln -s /usr/local/matlab712/ install
    1111#ln -s /usr/local/pkgs/matlab-7.6/ install
    1212#ln -s /Applications/MATLAB_R2008a/ install

  • issm/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h

    r11417 r11425  
    162162        ThermalSpctemperatureEnum,
    163163        ThermalStabilizationEnum,
    164         ThermalIsenthalpyEnum,
    165164        ThicknessEnum,
    166165        TimesteppingCflCoefficientEnum,

  • issm/trunk-jpl-damage/src/c/EnumDefinitions/Synchronize.sh

    r11417 r11425  
    99rm $ISSM_TIER/src/c/modules/EnumToStringx/EnumToStringx.cpp
    1010rm $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp
    11 
    12 #Get number of enums
    13 NUMENUMS=$(wc -l temp | awk '{printf("%s",$1);}');
    1411
    1512#Take care of EnumToModelField.m first (easy)
    … …  
    9895int  StringToEnumx(const char* name){
    9996
    100    int  stage=1;
     97END
     98#core
     99cat temp |  awk '{print "\t" ((NR==1)?"if":"else if") " (strcmp(name,\"" substr($2,1,length($2)-4) "\")==0) return " $2 ";"}' >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp
     100#Footer
     101cat <<END >> $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp
     102        else _error_("Enum %s not found",name);
    101103
    102 END
    103 
    104 #core
    105 i1=1;
    106 i2=120;
    107 for (( i=1 ; i<=100 ; i++ )); do
    108         echo "   if(stage==$i){" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp
    109         awk -v i1=$i1 -v i2=$i2 '{if(NR>=i1 && NR<=i2) print $0 }' temp |
    110         awk '{print "\t" ((NR==1)?"      if":"      else if") " (strcmp(name,\"" substr($2,1,length($2)-4) "\")==0) return " $2 ";"}' >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp
    111         echo "         else stage=$(($i+1));" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp
    112         echo "   }" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp
    113        
    114         if [ $i2 -ge $NUMENUMS ]; then break; fi
    115         let i1=$i1+120
    116         let i2=$i2+120
    117 done
    118 
    119 #footer
    120 cat <<END >> $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp
    121         /*If we reach this point, the string provided has not been found*/
    122    _error_("Enum %s not found",name);
    123104}
    124105END
    125106#}}}
    126107
     108#get number of lines in temp
     109NUMBEROFLINES=$(wc -l temp | awk '{printf("%s",$1);}');
     110
    127111# go through the lines of temp
    128 for (( i=1 ; i<=$NUMENUMS ; i++ )); do
     112for (( i=1 ; i<=$NUMBEROFLINES ; i++ )); do
    129113
    130114        #Get name and enum of the line i
    … …  
    139123        then
    140124                printf "\r                                                                      "
    141                 printf "\r  $i/$NUMENUMS Adding "$NAME"..."
     125                printf "\r  $i/$NUMBEROFLINES Adding "$NAME"..."
    142126        else
    143127                if [ $i -lt 100 ]
    144128                then
    145129                        printf "\r                                                                      "
    146                         printf "\r $i/$NUMENUMS Adding "$NAME"..."
     130                        printf "\r $i/$NUMBEROFLINES Adding "$NAME"..."
    147131                else
    148132                        printf "\r                                                                      "
    149                         printf "\r$i/$NUMENUMS Adding "$NAME"..."
     133                        printf "\r$i/$NUMBEROFLINES Adding "$NAME"..."
    150134                fi
    151135        fi
    … …  
    169153done
    170154
     155
    171156#clean up{{{
    172157rm temp

  • issm/trunk-jpl-damage/src/c/modules/BamgTriangulatex/BamgTriangulatex.cpp

    r11417 r11425  
    1717        Th.WriteIndex(pindex,pnels);
    1818        //delete &Th;
    19         return 0;
    2019
    2120}

  • issm/trunk-jpl-damage/src/c/modules/Chacox/Chacox.cpp

    r11417 r11425  
    191191
    192192
    193         #else //ifdef _HAVE_CHACO_
    194         return (0);
    195         #endif
     193        #endif //ifdef _HAVE_CHACO_
    196194}

  • issm/trunk-jpl-damage/src/c/modules/Chacox/input_parse.cpp

    r11417 r11425  
    262262        return(0);
    263263
    264         #else //#ifdef _HAVE_CHACO_
    265         return(0);
    266         #endif
     264        #endif //#ifdef _HAVE_CHACO_
    267265}

  • issm/trunk-jpl-damage/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp

    r11417 r11425  
    6969        }
    7070
    71         return NULL;
    72 
    7371}

  • issm/trunk-jpl-damage/src/c/modules/ContourToNodesx/ContourToNodesx.cpp

    r11417 r11425  
    3636        /*Assign output pointers: */
    3737        *pflags=flags;
    38        
    39         return 1;
     38
    4039}

  • issm/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r11417 r11425  
    166166                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
    167167                case ThermalStabilizationEnum : return "ThermalStabilization";
    168                 case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
    169168                case ThicknessEnum : return "Thickness";
    170169                case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient";

  • issm/trunk-jpl-damage/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp

    r11417 r11425  
    126126        }
    127127        if (debug && my_thread==0) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
    128        
    129         return NULL;
    130128
    131129}

  • issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r11417 r11425  
    7878        parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
    7979        parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
    80         parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
    8180        parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    8281        parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));

  • issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/ModelProcessorx.cpp

    r11417 r11425  
    2020
    2121        int   i,analysis_type,dim,verbose;
    22         bool  isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy;
     22        bool  isthermal,isprognostic,isdiagnostic,isgroundingline;
    2323       
    2424        /*output: */
    … …  
    3939        iomodel->Constant(&verbose,VerboseEnum);
    4040        iomodel->Constant(&isthermal,TransientIsthermalEnum);
    41         iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
    4241        iomodel->Constant(&isprognostic,TransientIsprognosticEnum);
    4342        iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum);
    … …  
    5352                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && dim==2) continue;
    5453                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && dim==2) continue;
    55                 if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && dim==2) continue;
    5654                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isthermal==false) continue;
    5755                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isthermal==false) continue;
    58                 if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isthermal==false) continue;
    59                 if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isenthalpy==true) continue;
    60                 if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isenthalpy==true) continue;
    61                 if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isenthalpy==false) continue;
    6256                if(solution_type==TransientSolutionEnum && analysis_type==PrognosticAnalysisEnum && isprognostic==false && isgroundingline==false) continue;
    6357                if(solution_type==TransientSolutionEnum && analysis_type==DiagnosticHorizAnalysisEnum && isdiagnostic==false) continue;

  • issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp

    r11417 r11425  
    3636        /*Assign output pointers: */
    3737        *pflags=flags;
    38 
    39         return 1;
    4038}

  • issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp

    r11417 r11425  
    7474        /*Free ressources:*/
    7575        xfree((void**)&already);
    76        
    77         return NULL;
    7876}

  • issm/trunk-jpl-damage/src/c/modules/Scotchx/Scotchx.cpp

    r11417 r11425  
    339339  return (0);
    340340
    341 #else //#ifdef _HAVE_SCOTCH_
    342   return(0);
    343 #endif
     341#endif //#ifdef _HAVE_SCOTCH_
    344342}

  • issm/trunk-jpl-damage/src/c/modules/Shp2Kmlx/Shp2Kmlx.cpp

    r11417 r11425  
    2424                                        sgn,cm,sp));
    2525
    26         #else //ifdef _HAVE_SHAPELIB_
    27         return 0;
    28         #endif
     26        #endif //ifdef _HAVE_SHAPELIB_
    2927}
    3028
    … …  
    609607        return(iret);
    610608
    611         #else //ifdef _HAVE_SHAPELIB_
    612         return 0;
    613         #endif
     609        #endif //ifdef _HAVE_SHAPELIB_
    614610}
    615611

  • issm/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r11418 r11425  
    1414int  StringToEnumx(const char* name){
    1515
    16    int  stage=1;
     16        if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum;
     17        else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum;
     18        else if (strcmp(name,"AutodiffReverse")==0) return AutodiffReverseEnum;
     19        else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
     20        else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum;
     21        else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum;
     22        else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
     23        else if (strcmp(name,"BasalforcingsMeltingRateCorrection")==0) return BasalforcingsMeltingRateCorrectionEnum;
     24        else if (strcmp(name,"BasalforcingsMeltingRate")==0) return BasalforcingsMeltingRateEnum;
     25        else if (strcmp(name,"Bathymetry")==0) return BathymetryEnum;
     26        else if (strcmp(name,"Bed")==0) return BedEnum;
     27        else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
     28        else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
     29        else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
     30        else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum;
     31        else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum;
     32        else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum;
     33        else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum;
     34        else if (strcmp(name,"DiagnosticPenaltyFactor")==0) return DiagnosticPenaltyFactorEnum;
     35        else if (strcmp(name,"DiagnosticReferential")==0) return DiagnosticReferentialEnum;
     36        else if (strcmp(name,"DiagnosticReltol")==0) return DiagnosticReltolEnum;
     37        else if (strcmp(name,"DiagnosticRequestedOutputs")==0) return DiagnosticRequestedOutputsEnum;
     38        else if (strcmp(name,"DiagnosticRestol")==0) return DiagnosticRestolEnum;
     39        else if (strcmp(name,"DiagnosticRiftPenaltyLock")==0) return DiagnosticRiftPenaltyLockEnum;
     40        else if (strcmp(name,"DiagnosticRiftPenaltyThreshold")==0) return DiagnosticRiftPenaltyThresholdEnum;
     41        else if (strcmp(name,"DiagnosticShelfDampening")==0) return DiagnosticShelfDampeningEnum;
     42        else if (strcmp(name,"DiagnosticSpcvx")==0) return DiagnosticSpcvxEnum;
     43        else if (strcmp(name,"DiagnosticSpcvy")==0) return DiagnosticSpcvyEnum;
     44        else if (strcmp(name,"DiagnosticSpcvz")==0) return DiagnosticSpcvzEnum;
     45        else if (strcmp(name,"DiagnosticStokesreconditioning")==0) return DiagnosticStokesreconditioningEnum;
     46        else if (strcmp(name,"DiagnosticVertexPairing")==0) return DiagnosticVertexPairingEnum;
     47        else if (strcmp(name,"DiagnosticViscosityOvershoot")==0) return DiagnosticViscosityOvershootEnum;
     48        else if (strcmp(name,"FlowequationBordermacayeal")==0) return FlowequationBordermacayealEnum;
     49        else if (strcmp(name,"FlowequationBorderpattyn")==0) return FlowequationBorderpattynEnum;
     50        else if (strcmp(name,"FlowequationBorderstokes")==0) return FlowequationBorderstokesEnum;
     51        else if (strcmp(name,"FlowequationElementEquation")==0) return FlowequationElementEquationEnum;
     52        else if (strcmp(name,"FlowequationIshutter")==0) return FlowequationIshutterEnum;
     53        else if (strcmp(name,"FlowequationIsmacayealpattyn")==0) return FlowequationIsmacayealpattynEnum;
     54        else if (strcmp(name,"FlowequationIsstokes")==0) return FlowequationIsstokesEnum;
     55        else if (strcmp(name,"FlowequationVertexEquation")==0) return FlowequationVertexEquationEnum;
     56        else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
     57        else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
     58        else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
     59        else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
     60        else if (strcmp(name,"HydrologyCR")==0) return HydrologyCREnum;
     61        else if (strcmp(name,"HydrologyKn")==0) return HydrologyKnEnum;
     62        else if (strcmp(name,"HydrologyN")==0) return HydrologyNEnum;
     63        else if (strcmp(name,"HydrologyP")==0) return HydrologyPEnum;
     64        else if (strcmp(name,"HydrologyQ")==0) return HydrologyQEnum;
     65        else if (strcmp(name,"HydrologySpcwatercolumn")==0) return HydrologySpcwatercolumnEnum;
     66        else if (strcmp(name,"HydrologyStabilization")==0) return HydrologyStabilizationEnum;
     67        else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
     68        else if (strcmp(name,"InversionCostFunction")==0) return InversionCostFunctionEnum;
     69        else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
     70        else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
     71        else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
     72        else if (strcmp(name,"InversionGradientOnly")==0) return InversionGradientOnlyEnum;
     73        else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
     74        else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
     75        else if (strcmp(name,"InversionTao")==0) return InversionTaoEnum;
     76        else if (strcmp(name,"InversionMaxParameters")==0) return InversionMaxParametersEnum;
     77        else if (strcmp(name,"InversionMaxiterPerStep")==0) return InversionMaxiterPerStepEnum;
     78        else if (strcmp(name,"InversionMinParameters")==0) return InversionMinParametersEnum;
     79        else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
     80        else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
     81        else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
     82        else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
     83        else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
     84        else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
     85        else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
     86        else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
     87        else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
     88        else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
     89        else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
     90        else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum;
     91        else if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum;
     92        else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
     93        else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
     94        else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
     95        else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
     96        else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
     97        else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
     98        else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
     99        else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
     100        else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
     101        else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
     102        else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
     103        else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum;
     104        else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum;
     105        else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
     106        else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
     107        else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
     108        else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
     109        else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
     110        else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
     111        else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum;
     112        else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum;
     113        else if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum;
     114        else if (strcmp(name,"MeshElementonbed")==0) return MeshElementonbedEnum;
     115        else if (strcmp(name,"MeshElementonsurface")==0) return MeshElementonsurfaceEnum;
     116        else if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum;
     117        else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
     118        else if (strcmp(name,"MeshLowerelements")==0) return MeshLowerelementsEnum;
     119        else if (strcmp(name,"MeshNumberofedges")==0) return MeshNumberofedgesEnum;
     120        else if (strcmp(name,"MeshNumberofelements2d")==0) return MeshNumberofelements2dEnum;
     121        else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum;
     122        else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
     123        else if (strcmp(name,"MeshNumberofvertices2d")==0) return MeshNumberofvertices2dEnum;
     124        else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
     125        else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum;
     126        else if (strcmp(name,"MeshVertexonbed")==0) return MeshVertexonbedEnum;
     127        else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
     128        else if (strcmp(name,"MeshX")==0) return MeshXEnum;
     129        else if (strcmp(name,"MeshY")==0) return MeshYEnum;
     130        else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
     131        else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
     132        else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum;
     133        else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
     134        else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
     135        else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
     136        else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
     137        else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
     138        else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
     139        else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum;
     140        else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
     141        else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
     142        else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
     143        else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
     144        else if (strcmp(name,"QmuVariabledescriptors")==0) return QmuVariabledescriptorsEnum;
     145        else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum;
     146        else if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum;
     147        else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
     148        else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
     149        else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
     150        else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum;
     151        else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
     152        else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
     153        else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
     154        else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
     155        else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
     156        else if (strcmp(name,"Surface")==0) return SurfaceEnum;
     157        else if (strcmp(name,"SurfaceforcingsAblationRate")==0) return SurfaceforcingsAblationRateEnum;
     158        else if (strcmp(name,"SurfaceforcingsAccumulationRate")==0) return SurfaceforcingsAccumulationRateEnum;
     159        else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
     160        else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
     161        else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
     162        else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     163        else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
     164        else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
     165        else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
     166        else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
     167        else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
     168        else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
     169        else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
     170        else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
     171        else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
     172        else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
     173        else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
     174        else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
     175        else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
     176        else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
     177        else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
     178        else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
     179        else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
     180        else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     181        else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
     182        else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
     183        else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
     184        else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
     185        else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
     186        else if (strcmp(name,"BedSlopeAnalysis")==0) return BedSlopeAnalysisEnum;
     187        else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     188        else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum;
     189        else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum;
     190        else if (strcmp(name,"DiagnosticHorizAnalysis")==0) return DiagnosticHorizAnalysisEnum;
     191        else if (strcmp(name,"DiagnosticHutterAnalysis")==0) return DiagnosticHutterAnalysisEnum;
     192        else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum;
     193        else if (strcmp(name,"DiagnosticVertAnalysis")==0) return DiagnosticVertAnalysisEnum;
     194        else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
     195        else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum;
     196        else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
     197        else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
     198        else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum;
     199        else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
     200        else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
     201        else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum;
     202        else if (strcmp(name,"PrognosticAnalysis")==0) return PrognosticAnalysisEnum;
     203        else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum;
     204        else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
     205        else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
     206        else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
     207        else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
     208        else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
     209        else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
     210        else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
     211        else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
     212        else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
     213        else if (strcmp(name,"HutterApproximation")==0) return HutterApproximationEnum;
     214        else if (strcmp(name,"MacAyealApproximation")==0) return MacAyealApproximationEnum;
     215        else if (strcmp(name,"MacAyealPattynApproximation")==0) return MacAyealPattynApproximationEnum;
     216        else if (strcmp(name,"MacAyealStokesApproximation")==0) return MacAyealStokesApproximationEnum;
     217        else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
     218        else if (strcmp(name,"PattynApproximation")==0) return PattynApproximationEnum;
     219        else if (strcmp(name,"PattynStokesApproximation")==0) return PattynStokesApproximationEnum;
     220        else if (strcmp(name,"StokesApproximation")==0) return StokesApproximationEnum;
     221        else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     222        else if (strcmp(name,"Loads")==0) return LoadsEnum;
     223        else if (strcmp(name,"Materials")==0) return MaterialsEnum;
     224        else if (strcmp(name,"Nodes")==0) return NodesEnum;
     225        else if (strcmp(name,"Parameters")==0) return ParametersEnum;
     226        else if (strcmp(name,"Vertices")==0) return VerticesEnum;
     227        else if (strcmp(name,"Results")==0) return ResultsEnum;
     228        else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     229        else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     230        else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
     231        else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
     232        else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
     233        else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
     234        else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
     235        else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
     236        else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
     237        else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
     238        else if (strcmp(name,"Element")==0) return ElementEnum;
     239        else if (strcmp(name,"ElementResult")==0) return ElementResultEnum;
     240        else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
     241        else if (strcmp(name,"FileParam")==0) return FileParamEnum;
     242        else if (strcmp(name,"Hook")==0) return HookEnum;
     243        else if (strcmp(name,"Icefront")==0) return IcefrontEnum;
     244        else if (strcmp(name,"Input")==0) return InputEnum;
     245        else if (strcmp(name,"IntInput")==0) return IntInputEnum;
     246        else if (strcmp(name,"IntParam")==0) return IntParamEnum;
     247        else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
     248        else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
     249        else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
     250        else if (strcmp(name,"Matice")==0) return MaticeEnum;
     251        else if (strcmp(name,"Matpar")==0) return MatparEnum;
     252        else if (strcmp(name,"Node")==0) return NodeEnum;
     253        else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
     254        else if (strcmp(name,"Param")==0) return ParamEnum;
     255        else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
     256        else if (strcmp(name,"Pengrid")==0) return PengridEnum;
     257        else if (strcmp(name,"Penpair")==0) return PenpairEnum;
     258        else if (strcmp(name,"Penta")==0) return PentaEnum;
     259        else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum;
     260        else if (strcmp(name,"PetscMatParam")==0) return PetscMatParamEnum;
     261        else if (strcmp(name,"PetscVecParam")==0) return PetscVecParamEnum;
     262        else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     263        else if (strcmp(name,"Segment")==0) return SegmentEnum;
     264        else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     265        else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
     266        else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
     267        else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
     268        else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
     269        else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
     270        else if (strcmp(name,"StringParam")==0) return StringParamEnum;
     271        else if (strcmp(name,"Tria")==0) return TriaEnum;
     272        else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
     273        else if (strcmp(name,"Vertex")==0) return VertexEnum;
     274        else if (strcmp(name,"Air")==0) return AirEnum;
     275        else if (strcmp(name,"Ice")==0) return IceEnum;
     276        else if (strcmp(name,"Melange")==0) return MelangeEnum;
     277        else if (strcmp(name,"Water")==0) return WaterEnum;
     278        else if (strcmp(name,"Closed")==0) return ClosedEnum;
     279        else if (strcmp(name,"Free")==0) return FreeEnum;
     280        else if (strcmp(name,"Open")==0) return OpenEnum;
     281        else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
     282        else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
     283        else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
     284        else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
     285        else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
     286        else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
     287        else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     288        else if (strcmp(name,"Constant")==0) return ConstantEnum;
     289        else if (strcmp(name,"Converged")==0) return ConvergedEnum;
     290        else if (strcmp(name,"ExtToIu")==0) return ExtToIuEnum;
     291        else if (strcmp(name,"Fill")==0) return FillEnum;
     292        else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
     293        else if (strcmp(name,"Friction")==0) return FrictionEnum;
     294        else if (strcmp(name,"GroundinglineMeltingRate")==0) return GroundinglineMeltingRateEnum;
     295        else if (strcmp(name,"Internal")==0) return InternalEnum;
     296        else if (strcmp(name,"IuToExt")==0) return IuToExtEnum;
     297        else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
     298        else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum;
     299        else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
     300        else if (strcmp(name,"Misfit")==0) return MisfitEnum;
     301        else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum;
     302        else if (strcmp(name,"Pressure")==0) return PressureEnum;
     303        else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
     304        else if (strcmp(name,"QmuPressure")==0) return QmuPressureEnum;
     305        else if (strcmp(name,"QmuVx")==0) return QmuVxEnum;
     306        else if (strcmp(name,"QmuVy")==0) return QmuVyEnum;
     307        else if (strcmp(name,"QmuVz")==0) return QmuVzEnum;
     308        else if (strcmp(name,"QmuThickness")==0) return QmuThicknessEnum;
     309        else if (strcmp(name,"QmuBed")==0) return QmuBedEnum;
     310        else if (strcmp(name,"QmuSurface")==0) return QmuSurfaceEnum;
     311        else if (strcmp(name,"QmuMelting")==0) return QmuMeltingEnum;
     312        else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
     313        else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
     314        else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
     315        else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
     316        else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
     317        else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
     318        else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     319        else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
     320        else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
     321        else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
     322        else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
     323        else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
     324        else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
     325        else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
     326        else if (strcmp(name,"Type")==0) return TypeEnum;
     327        else if (strcmp(name,"Vel")==0) return VelEnum;
     328        else if (strcmp(name,"Velocity")==0) return VelocityEnum;
     329        else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     330        else if (strcmp(name,"Vx")==0) return VxEnum;
     331        else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
     332        else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     333        else if (strcmp(name,"Vy")==0) return VyEnum;
     334        else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
     335        else if (strcmp(name,"Vz")==0) return VzEnum;
     336        else if (strcmp(name,"VzMacAyeal")==0) return VzMacAyealEnum;
     337        else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum;
     338        else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
     339        else if (strcmp(name,"VzStokes")==0) return VzStokesEnum;
     340        else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
     341        else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
     342        else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
     343        else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
     344        else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
     345        else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
     346        else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
     347        else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
     348        else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
     349        else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
     350        else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
     351        else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum;
     352        else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
     353        else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
     354        else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
     355        else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
     356        else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
     357        else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
     358        else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
     359        else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
     360        else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
     361        else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
     362        else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
     363        else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     364        else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
     365        else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
     366        else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     367        else if (strcmp(name,"P0")==0) return P0Enum;
     368        else if (strcmp(name,"P1")==0) return P1Enum;
     369        else if (strcmp(name,"P1DG")==0) return P1DGEnum;
     370        else if (strcmp(name,"BoolElementResult")==0) return BoolElementResultEnum;
     371        else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     372        else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
     373        else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     374        else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     375        else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum;
     376        else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     377        else if (strcmp(name,"J")==0) return JEnum;
     378        else if (strcmp(name,"Patch")==0) return PatchEnum;
     379        else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum;
     380        else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum;
     381        else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
     382        else if (strcmp(name,"PetscVecExternalResult")==0) return PetscVecExternalResultEnum;
     383        else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
     384        else if (strcmp(name,"Time")==0) return TimeEnum;
     385        else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
     386        else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
     387        else if (strcmp(name,"MinVel")==0) return MinVelEnum;
     388        else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
     389        else if (strcmp(name,"MinVx")==0) return MinVxEnum;
     390        else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
     391        else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
     392        else if (strcmp(name,"MinVy")==0) return MinVyEnum;
     393        else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
     394        else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
     395        else if (strcmp(name,"MinVz")==0) return MinVzEnum;
     396        else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
     397        else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
     398        else if (strcmp(name,"Relative")==0) return RelativeEnum;
     399        else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
     400        else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     401        else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
     402        else if (strcmp(name,"None")==0) return NoneEnum;
     403        else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
     404        else if (strcmp(name,"StokesSolver")==0) return StokesSolverEnum;
     405        else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
     406        else if (strcmp(name,"Colinear")==0) return ColinearEnum;
     407        else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
     408        else if (strcmp(name,"Fset")==0) return FsetEnum;
     409        else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
     410        else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
     411        else if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
     412        else if (strcmp(name,"Gradient")==0) return GradientEnum;
     413        else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
     414        else if (strcmp(name,"Gset")==0) return GsetEnum;
     415        else if (strcmp(name,"Index")==0) return IndexEnum;
     416        else if (strcmp(name,"Indexed")==0) return IndexedEnum;
     417        else if (strcmp(name,"Intersect")==0) return IntersectEnum;
     418        else if (strcmp(name,"Nodal")==0) return NodalEnum;
     419        else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
     420        else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
     421        else if (strcmp(name,"PetscOptionsAnalyses")==0) return PetscOptionsAnalysesEnum;
     422        else if (strcmp(name,"PetscOptionsStrings")==0) return PetscOptionsStringsEnum;
     423        else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
     424        else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
     425        else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
     426        else if (strcmp(name,"Regular")==0) return RegularEnum;
     427        else if (strcmp(name,"Scaled")==0) return ScaledEnum;
     428        else if (strcmp(name,"Separate")==0) return SeparateEnum;
     429        else if (strcmp(name,"Sset")==0) return SsetEnum;
     430        else if (strcmp(name,"Verbose")==0) return VerboseEnum;
     431        else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
     432        else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     433        else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
     434        else if (strcmp(name,"XY")==0) return XYEnum;
     435        else if (strcmp(name,"XYZP")==0) return XYZPEnum;
     436        else if (strcmp(name,"Option")==0) return OptionEnum;
     437        else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
     438        else if (strcmp(name,"OptionChar")==0) return OptionCharEnum;
     439        else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
     440        else if (strcmp(name,"OptionDouble")==0) return OptionDoubleEnum;
     441        else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum;
     442        else if (strcmp(name,"Paterson")==0) return PatersonEnum;
     443        else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
     444        else _error_("Enum %s not found",name);
    17445
    18    if(stage==1){
    19               if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum;
    20               else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum;
    21               else if (strcmp(name,"AutodiffReverse")==0) return AutodiffReverseEnum;
    22               else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
    23               else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum;
    24               else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum;
    25               else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
    26               else if (strcmp(name,"BasalforcingsMeltingRateCorrection")==0) return BasalforcingsMeltingRateCorrectionEnum;
    27               else if (strcmp(name,"BasalforcingsMeltingRate")==0) return BasalforcingsMeltingRateEnum;
    28               else if (strcmp(name,"Bathymetry")==0) return BathymetryEnum;
    29               else if (strcmp(name,"Bed")==0) return BedEnum;
    30               else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
    31               else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
    32               else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
    33               else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum;
    34               else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum;
    35               else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum;
    36               else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum;
    37               else if (strcmp(name,"DiagnosticPenaltyFactor")==0) return DiagnosticPenaltyFactorEnum;
    38               else if (strcmp(name,"DiagnosticReferential")==0) return DiagnosticReferentialEnum;
    39               else if (strcmp(name,"DiagnosticReltol")==0) return DiagnosticReltolEnum;
    40               else if (strcmp(name,"DiagnosticRequestedOutputs")==0) return DiagnosticRequestedOutputsEnum;
    41               else if (strcmp(name,"DiagnosticRestol")==0) return DiagnosticRestolEnum;
    42               else if (strcmp(name,"DiagnosticRiftPenaltyLock")==0) return DiagnosticRiftPenaltyLockEnum;
    43               else if (strcmp(name,"DiagnosticRiftPenaltyThreshold")==0) return DiagnosticRiftPenaltyThresholdEnum;
    44               else if (strcmp(name,"DiagnosticShelfDampening")==0) return DiagnosticShelfDampeningEnum;
    45               else if (strcmp(name,"DiagnosticSpcvx")==0) return DiagnosticSpcvxEnum;
    46               else if (strcmp(name,"DiagnosticSpcvy")==0) return DiagnosticSpcvyEnum;
    47               else if (strcmp(name,"DiagnosticSpcvz")==0) return DiagnosticSpcvzEnum;
    48               else if (strcmp(name,"DiagnosticStokesreconditioning")==0) return DiagnosticStokesreconditioningEnum;
    49               else if (strcmp(name,"DiagnosticVertexPairing")==0) return DiagnosticVertexPairingEnum;
    50               else if (strcmp(name,"DiagnosticViscosityOvershoot")==0) return DiagnosticViscosityOvershootEnum;
    51               else if (strcmp(name,"FlowequationBordermacayeal")==0) return FlowequationBordermacayealEnum;
    52               else if (strcmp(name,"FlowequationBorderpattyn")==0) return FlowequationBorderpattynEnum;
    53               else if (strcmp(name,"FlowequationBorderstokes")==0) return FlowequationBorderstokesEnum;
    54               else if (strcmp(name,"FlowequationElementEquation")==0) return FlowequationElementEquationEnum;
    55               else if (strcmp(name,"FlowequationIshutter")==0) return FlowequationIshutterEnum;
    56               else if (strcmp(name,"FlowequationIsmacayealpattyn")==0) return FlowequationIsmacayealpattynEnum;
    57               else if (strcmp(name,"FlowequationIsstokes")==0) return FlowequationIsstokesEnum;
    58               else if (strcmp(name,"FlowequationVertexEquation")==0) return FlowequationVertexEquationEnum;
    59               else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
    60               else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
    61               else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    62               else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
    63               else if (strcmp(name,"HydrologyCR")==0) return HydrologyCREnum;
    64               else if (strcmp(name,"HydrologyKn")==0) return HydrologyKnEnum;
    65               else if (strcmp(name,"HydrologyN")==0) return HydrologyNEnum;
    66               else if (strcmp(name,"HydrologyP")==0) return HydrologyPEnum;
    67               else if (strcmp(name,"HydrologyQ")==0) return HydrologyQEnum;
    68               else if (strcmp(name,"HydrologySpcwatercolumn")==0) return HydrologySpcwatercolumnEnum;
    69               else if (strcmp(name,"HydrologyStabilization")==0) return HydrologyStabilizationEnum;
    70               else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    71               else if (strcmp(name,"InversionCostFunction")==0) return InversionCostFunctionEnum;
    72               else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
    73               else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
    74               else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
    75               else if (strcmp(name,"InversionGradientOnly")==0) return InversionGradientOnlyEnum;
    76               else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
    77               else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
    78               else if (strcmp(name,"InversionTao")==0) return InversionTaoEnum;
    79               else if (strcmp(name,"InversionMaxParameters")==0) return InversionMaxParametersEnum;
    80               else if (strcmp(name,"InversionMaxiterPerStep")==0) return InversionMaxiterPerStepEnum;
    81               else if (strcmp(name,"InversionMinParameters")==0) return InversionMinParametersEnum;
    82               else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
    83               else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
    84               else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
    85               else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    86               else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
    87               else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    88               else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    89               else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    90               else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    91               else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
    92               else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
    93               else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum;
    94               else if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum;
    95               else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
    96               else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
    97               else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    98               else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
    99               else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
    100               else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
    101               else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
    102               else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
    103               else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
    104               else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
    105               else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
    106               else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum;
    107               else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum;
    108               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    109               else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
    110               else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
    111               else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
    112               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
    113               else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
    114               else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum;
    115               else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum;
    116               else if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum;
    117               else if (strcmp(name,"MeshElementonbed")==0) return MeshElementonbedEnum;
    118               else if (strcmp(name,"MeshElementonsurface")==0) return MeshElementonsurfaceEnum;
    119               else if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum;
    120               else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
    121               else if (strcmp(name,"MeshLowerelements")==0) return MeshLowerelementsEnum;
    122               else if (strcmp(name,"MeshNumberofedges")==0) return MeshNumberofedgesEnum;
    123               else if (strcmp(name,"MeshNumberofelements2d")==0) return MeshNumberofelements2dEnum;
    124               else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum;
    125               else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
    126               else if (strcmp(name,"MeshNumberofvertices2d")==0) return MeshNumberofvertices2dEnum;
    127               else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
    128               else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum;
    129               else if (strcmp(name,"MeshVertexonbed")==0) return MeshVertexonbedEnum;
    130               else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
    131               else if (strcmp(name,"MeshX")==0) return MeshXEnum;
    132               else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    133               else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    134               else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
    135               else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum;
    136               else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
    137               else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
    138               else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
    139          else stage=2;
    140    }
    141    if(stage==2){
    142               if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
    143               else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
    144               else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
    145               else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum;
    146               else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
    147               else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
    148               else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
    149               else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
    150               else if (strcmp(name,"QmuVariabledescriptors")==0) return QmuVariabledescriptorsEnum;
    151               else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum;
    152               else if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum;
    153               else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
    154               else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
    155               else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
    156               else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum;
    157               else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    158               else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
    159               else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    160               else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    161               else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    162               else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    163               else if (strcmp(name,"SurfaceforcingsAblationRate")==0) return SurfaceforcingsAblationRateEnum;
    164               else if (strcmp(name,"SurfaceforcingsAccumulationRate")==0) return SurfaceforcingsAccumulationRateEnum;
    165               else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
    166               else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
    167               else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    168               else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    169               else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    170               else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
    171               else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
    172               else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    173               else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    174               else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    175               else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    176               else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
    177               else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
    178               else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
    179               else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    180               else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
    181               else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
    182               else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
    183               else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    184               else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    185               else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
    186               else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
    187               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    188               else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    189               else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
    190               else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
    191               else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    192               else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    193               else if (strcmp(name,"BedSlopeAnalysis")==0) return BedSlopeAnalysisEnum;
    194               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    195               else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum;
    196               else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum;
    197               else if (strcmp(name,"DiagnosticHorizAnalysis")==0) return DiagnosticHorizAnalysisEnum;
    198               else if (strcmp(name,"DiagnosticHutterAnalysis")==0) return DiagnosticHutterAnalysisEnum;
    199               else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum;
    200               else if (strcmp(name,"DiagnosticVertAnalysis")==0) return DiagnosticVertAnalysisEnum;
    201               else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    202               else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum;
    203               else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
    204               else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
    205               else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum;
    206               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    207               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    208               else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum;
    209               else if (strcmp(name,"PrognosticAnalysis")==0) return PrognosticAnalysisEnum;
    210               else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum;
    211               else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    212               else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
    213               else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    214               else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
    215               else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
    216               else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    217               else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    218               else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    219               else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
    220               else if (strcmp(name,"HutterApproximation")==0) return HutterApproximationEnum;
    221               else if (strcmp(name,"MacAyealApproximation")==0) return MacAyealApproximationEnum;
    222               else if (strcmp(name,"MacAyealPattynApproximation")==0) return MacAyealPattynApproximationEnum;
    223               else if (strcmp(name,"MacAyealStokesApproximation")==0) return MacAyealStokesApproximationEnum;
    224               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    225               else if (strcmp(name,"PattynApproximation")==0) return PattynApproximationEnum;
    226               else if (strcmp(name,"PattynStokesApproximation")==0) return PattynStokesApproximationEnum;
    227               else if (strcmp(name,"StokesApproximation")==0) return StokesApproximationEnum;
    228               else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    229               else if (strcmp(name,"Loads")==0) return LoadsEnum;
    230               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    231               else if (strcmp(name,"Nodes")==0) return NodesEnum;
    232               else if (strcmp(name,"Parameters")==0) return ParametersEnum;
    233               else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    234               else if (strcmp(name,"Results")==0) return ResultsEnum;
    235               else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    236               else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    237               else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    238               else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    239               else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
    240               else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    241               else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    242               else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    243               else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    244               else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    245               else if (strcmp(name,"Element")==0) return ElementEnum;
    246               else if (strcmp(name,"ElementResult")==0) return ElementResultEnum;
    247               else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
    248               else if (strcmp(name,"FileParam")==0) return FileParamEnum;
    249               else if (strcmp(name,"Hook")==0) return HookEnum;
    250               else if (strcmp(name,"Icefront")==0) return IcefrontEnum;
    251               else if (strcmp(name,"Input")==0) return InputEnum;
    252               else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    253               else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    254               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    255               else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
    256               else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
    257               else if (strcmp(name,"Matice")==0) return MaticeEnum;
    258               else if (strcmp(name,"Matpar")==0) return MatparEnum;
    259               else if (strcmp(name,"Node")==0) return NodeEnum;
    260               else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    261               else if (strcmp(name,"Param")==0) return ParamEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
    266               else if (strcmp(name,"Pengrid")==0) return PengridEnum;
    267               else if (strcmp(name,"Penpair")==0) return PenpairEnum;
    268               else if (strcmp(name,"Penta")==0) return PentaEnum;
    269               else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum;
    270               else if (strcmp(name,"PetscMatParam")==0) return PetscMatParamEnum;
    271               else if (strcmp(name,"PetscVecParam")==0) return PetscVecParamEnum;
    272               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    273               else if (strcmp(name,"Segment")==0) return SegmentEnum;
    274               else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    275               else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
    276               else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    277               else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    278               else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
    279               else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    280               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    281               else if (strcmp(name,"Tria")==0) return TriaEnum;
    282               else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
    283               else if (strcmp(name,"Vertex")==0) return VertexEnum;
    284               else if (strcmp(name,"Air")==0) return AirEnum;
    285               else if (strcmp(name,"Ice")==0) return IceEnum;
    286               else if (strcmp(name,"Melange")==0) return MelangeEnum;
    287               else if (strcmp(name,"Water")==0) return WaterEnum;
    288               else if (strcmp(name,"Closed")==0) return ClosedEnum;
    289               else if (strcmp(name,"Free")==0) return FreeEnum;
    290               else if (strcmp(name,"Open")==0) return OpenEnum;
    291               else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    292               else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    293               else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    294               else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
    295               else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    296               else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    297               else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    298               else if (strcmp(name,"Constant")==0) return ConstantEnum;
    299               else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    300               else if (strcmp(name,"ExtToIu")==0) return ExtToIuEnum;
    301               else if (strcmp(name,"Fill")==0) return FillEnum;
    302               else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    303               else if (strcmp(name,"Friction")==0) return FrictionEnum;
    304               else if (strcmp(name,"GroundinglineMeltingRate")==0) return GroundinglineMeltingRateEnum;
    305               else if (strcmp(name,"Internal")==0) return InternalEnum;
    306               else if (strcmp(name,"IuToExt")==0) return IuToExtEnum;
    307               else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    308               else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum;
    309               else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    310               else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    311               else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum;
    312               else if (strcmp(name,"Pressure")==0) return PressureEnum;
    313               else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
    314               else if (strcmp(name,"QmuPressure")==0) return QmuPressureEnum;
    315               else if (strcmp(name,"QmuVx")==0) return QmuVxEnum;
    316               else if (strcmp(name,"QmuVy")==0) return QmuVyEnum;
    317               else if (strcmp(name,"QmuVz")==0) return QmuVzEnum;
    318               else if (strcmp(name,"QmuThickness")==0) return QmuThicknessEnum;
    319               else if (strcmp(name,"QmuBed")==0) return QmuBedEnum;
    320               else if (strcmp(name,"QmuSurface")==0) return QmuSurfaceEnum;
    321               else if (strcmp(name,"QmuMelting")==0) return QmuMeltingEnum;
    322               else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
    323               else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
    324               else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    325               else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    326               else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    327               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    328               else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    329               else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    330               else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
    331               else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
    332               else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
    333               else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
    334               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    335               else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    336               else if (strcmp(name,"Type")==0) return TypeEnum;
    337               else if (strcmp(name,"Vel")==0) return VelEnum;
    338               else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    339               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    340               else if (strcmp(name,"Vx")==0) return VxEnum;
    341               else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
    342               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    343               else if (strcmp(name,"Vy")==0) return VyEnum;
    344               else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
    345               else if (strcmp(name,"Vz")==0) return VzEnum;
    346               else if (strcmp(name,"VzMacAyeal")==0) return VzMacAyealEnum;
    347               else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum;
    348               else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
    349               else if (strcmp(name,"VzStokes")==0) return VzStokesEnum;
    350               else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    351               else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    352               else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    353               else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    354               else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    355               else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    356               else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
    357               else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    358               else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    359               else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
    360               else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    361               else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum;
    362               else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    363               else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
    364               else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
    365               else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    366               else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
    367               else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    368               else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    369               else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
    370               else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
    371               else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    372               else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    373               else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    374               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    375               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    376               else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    377               else if (strcmp(name,"P0")==0) return P0Enum;
    378               else if (strcmp(name,"P1")==0) return P1Enum;
    379               else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    380               else if (strcmp(name,"BoolElementResult")==0) return BoolElementResultEnum;
    381               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    382               else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
    383               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    384               else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum;
    389               else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    390               else if (strcmp(name,"J")==0) return JEnum;
    391               else if (strcmp(name,"Patch")==0) return PatchEnum;
    392               else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum;
    393               else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum;
    394               else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
    395               else if (strcmp(name,"PetscVecExternalResult")==0) return PetscVecExternalResultEnum;
    396               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    397               else if (strcmp(name,"Time")==0) return TimeEnum;
    398               else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
    399               else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    400               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    401               else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    402               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    403               else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    404               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    405               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    406               else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    407               else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    408               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    409               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
    410               else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
    411               else if (strcmp(name,"Relative")==0) return RelativeEnum;
    412               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    413               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    414               else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
    415               else if (strcmp(name,"None")==0) return NoneEnum;
    416               else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
    417               else if (strcmp(name,"StokesSolver")==0) return StokesSolverEnum;
    418               else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
    419               else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    420               else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
    421               else if (strcmp(name,"Fset")==0) return FsetEnum;
    422               else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
    423               else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
    424               else if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
    425               else if (strcmp(name,"Gradient")==0) return GradientEnum;
    426               else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
    427               else if (strcmp(name,"Gset")==0) return GsetEnum;
    428               else if (strcmp(name,"Index")==0) return IndexEnum;
    429               else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    430               else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    431               else if (strcmp(name,"Nodal")==0) return NodalEnum;
    432               else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
    433               else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
    434               else if (strcmp(name,"PetscOptionsAnalyses")==0) return PetscOptionsAnalysesEnum;
    435               else if (strcmp(name,"PetscOptionsStrings")==0) return PetscOptionsStringsEnum;
    436               else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
    437               else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
    438               else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
    439               else if (strcmp(name,"Regular")==0) return RegularEnum;
    440               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
    441               else if (strcmp(name,"Separate")==0) return SeparateEnum;
    442               else if (strcmp(name,"Sset")==0) return SsetEnum;
    443               else if (strcmp(name,"Verbose")==0) return VerboseEnum;
    444               else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    445               else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    446               else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    447               else if (strcmp(name,"XY")==0) return XYEnum;
    448               else if (strcmp(name,"XYZP")==0) return XYZPEnum;
    449               else if (strcmp(name,"Option")==0) return OptionEnum;
    450               else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
    451               else if (strcmp(name,"OptionChar")==0) return OptionCharEnum;
    452               else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
    453               else if (strcmp(name,"OptionDouble")==0) return OptionDoubleEnum;
    454               else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum;
    455               else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    456               else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    457          else stage=5;
    458    }
    459         /*If we reach this point, the string provided has not been found*/
    460    _error_("Enum %s not found",name);
    461446}

  • issm/trunk-jpl-damage/src/c/objects/Elements/Penta.cpp

    r11419 r11425  
    769769}
    770770/*}}}*/
    771 /*FUNCTION Penta::CreateJacobianMatrix{{{1*/
    772 void  Penta::CreateJacobianMatrix(Mat Jff){
    773 
    774         /*retrieve parameters: */
    775         ElementMatrix* Ke=NULL;
    776         int analysis_type;
    777         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    778 
    779         /*Checks in debugging {{{2*/
    780         _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    781         /*}}}*/
    782 
    783         /*Skip if water element*/
    784         if(IsOnWater()) return;
    785 
    786         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    787         switch(analysis_type){
    788 #ifdef _HAVE_DIAGNOSTIC_
    789                 case DiagnosticHorizAnalysisEnum:
    790                         Ke=CreateJacobianDiagnosticHoriz();
    791                         break;
    792 #endif
    793                 default:
    794                         _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    795         }
    796 
    797         /*Add to global matrix*/
    798         if(Ke){
    799                 Ke->AddToGlobal(Jff);
    800                 delete Ke;
    801         }
    802 }
    803 /*}}}*/
    804771/*FUNCTION Penta::DeepEcho{{{1*/
    805772void Penta::DeepEcho(void){
    … …  
    11321099/*}}}*/
    11331100/*FUNCTION Penta::GetStabilizationParameter {{{1*/
    1134 double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){
     1101double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
    11351102        /*Compute stabilization parameter*/
    1136         /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
    1137         /*kappa=enthalpydiffusionparameter for enthalpy model*/
    11381103
    11391104        double normu;
    … …  
    11411106
    11421107        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
    1143         if(normu*diameter/(3*2*kappa)<1){
    1144                 tau_parameter=pow(diameter,2)/(3*2*2*kappa);
     1108        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
     1109                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
    11451110        }
    11461111        else tau_parameter=diameter/(2*normu);
    … …  
    33183283                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    33193284
    3320                 vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    3321                 vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    3322                 vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
     3285                vx_input->GetInputValue(&u, gauss);
     3286                vy_input->GetInputValue(&v, gauss);
     3287                vz_input->GetInputValue(&w, gauss);
     3288                vxm_input->GetInputValue(&um,gauss);
     3289                vym_input->GetInputValue(&vm,gauss);
     3290                vzm_input->GetInputValue(&wm,gauss);
     3291                vx=u-um; vy=v-vm; vz=w-wm;
    33233292
    33243293                D_scalar_advec=gauss->weight*Jdet;
    … …  
    33523321                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    33533322                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
    3354                         K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
    3355                         K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
    3356                         K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
     3323                        K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
     3324                        K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
     3325                        K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
    33573326                        D_scalar_stab=gauss->weight*Jdet;
    33583327                        if(dt) D_scalar_stab=D_scalar_stab*dt;
    … …  
    33683337                else if(stabilization==2){
    33693338                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3370                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
     3339                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    33713340
    33723341                        for(i=0;i<numdof;i++){
    … …  
    34823451        double     Jdet,u,v,w,um,vm,wm,vel;
    34833452        double     h,hx,hy,hz,vx,vy,vz;
    3484         double     gravity,rho_ice,rho_water,kappa;
     3453        double     gravity,rho_ice,rho_water;
    34853454        double     heatcapacity,thermalconductivity,dt;
    34863455        double     tau_parameter,diameter;
    … …  
    35083477        heatcapacity=matpar->GetHeatCapacity();
    35093478        thermalconductivity=matpar->GetThermalConductivity();
    3510         kappa=thermalconductivity/(rho_ice*heatcapacity);
    35113479        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    35123480        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    … …  
    35313499                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    35323500
    3533                 D_scalar_conduct=gauss->weight*Jdet*kappa;
     3501                D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
    35343502                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
    35353503
    … …  
    36003568                else if(stabilization==2){
    36013569                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3602                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
     3570                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    36033571
    36043572                        for(i=0;i<numdof;i++){
    … …  
    37053673        double Jdet,phi,dt;
    37063674        double rho_ice,heatcapacity;
    3707         double thermalconductivity,kappa;
    3708         double viscosity,pressure;
    3709         double enthalpy,enthalpypicard;
     3675        double thermalconductivity;
     3676        double viscosity,enthalpy;
    37103677        double tau_parameter,diameter;
    37113678        double u,v,w;
    … …  
    37283695        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    37293696        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3730         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    3731         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    3732         Input* vz_input=inputs->GetInput(VzEnum);                                  _assert_(vz_input);
    3733         Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
    3734         Input* enthalpy_input=NULL;
    3735         Input* enthalpypicard_input=NULL;
    3736         if(dt){
    3737                 enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
    3738         }
    3739         if (stabilization==2){
    3740                 diameter=MinEdgeLength(xyz_list);
    3741                 enthalpypicard_input=inputs->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
    3742         }
     3697        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3698        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3699        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     3700        Input* enthalpy_input=NULL;
     3701        if (dt) enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(inputs);
     3702        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
    37433703
    37443704        /* Start  looping on the number of gaussian points: */
    … …  
    37553715                GetPhi(&phi, &epsilon[0], viscosity);
    37563716
    3757                 scalar_def=phi/rho_ice*Jdet*gauss->weight;
     3717                scalar_def=phi/(rho_ice)*Jdet*gauss->weight;
    37583718                if(dt) scalar_def=scalar_def*dt;
    37593719
    … …  
    37733733                        vy_input->GetInputValue(&v, gauss);
    37743734                        vz_input->GetInputValue(&w, gauss);
    3775                         pressure_input->GetInputValue(&pressure, gauss);
    3776                         enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
    3777                         kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    3778                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
     3735
     3736                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
    37793737
    37803738                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    … …  
    38613819        double     rho_ice,heatcapacity,geothermalflux_value;
    38623820        double     basalfriction,alpha2,vx,vy;
    3863         double     scalar,enthalpy,enthalpyup;
    3864         double     pressure,pressureup;
     3821        double     scalar;
    38653822        double     basis[NUMVERTICES];
    3866         Friction*  friction=NULL;
    3867         GaussPenta* gauss=NULL;
    3868         GaussPenta* gaussup=NULL;
    3869 
    3870         /* Geothermal flux on ice sheet base and basal friction */
    3871         if (!IsOnBed() || IsFloating()) return NULL;
    3872 
    3873         /*Initialize Element vector*/
    3874         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    3875 
    3876         /*Retrieve all inputs and parameters*/
    3877         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3878         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    3879         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3880         rho_ice=matpar->GetRhoIce();
    3881         heatcapacity=matpar->GetHeatCapacity();
    3882         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3883         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    3884         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    3885         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
    3886         Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
    3887         Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    3888         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    3889 
    3890         /*Build frictoin element, needed later: */
    3891         friction=new Friction("3d",inputs,matpar,analysis_type);
    3892 
    3893         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    3894         gauss=new GaussPenta(0,1,2,2);
    3895         gaussup=new GaussPenta(3,4,5,2);
    3896         for(ig=gauss->begin();ig<gauss->end();ig++){
    3897 
    3898                 gauss->GaussPoint(ig);
    3899                 gaussup->GaussPoint(ig);
    3900 
    3901                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3902                 GetNodalFunctionsP1(&basis[0], gauss);
    3903 
    3904                 enthalpy_input->GetInputValue(&enthalpy,gauss);
    3905                 pressure_input->GetInputValue(&pressure,gauss);
    3906 //              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
    3907 //                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    3908 //                      pressure_input->GetInputValue(&pressureup,gaussup);
    3909 //                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
    3910 //                              //do nothing, don't add heatflux
    3911 //                      }
    3912 //                      else{
    3913 //                              //need to change spcenthalpy according to Aschwanden
    3914 //                              //NEED TO UPDATE
    3915 //                      }
    3916 //              }
    3917 //              else{
    3918                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    3919                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    3920                         vx_input->GetInputValue(&vx,gauss);
    3921                         vy_input->GetInputValue(&vy,gauss);
    3922                         basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    3923 
    3924                         scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
    3925                         if(dt) scalar=dt*scalar;
    3926 
    3927                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    3928 //              }
    3929         }
    3930 
    3931         /*Clean up and return*/
    3932         delete gauss;
    3933         delete gaussup;
    3934         delete friction;
    3935         return pe;
    3936 }
    3937 /*}}}*/
    3938 /*FUNCTION Penta::CreatePVectorMelting {{{1*/
    3939 ElementVector* Penta::CreatePVectorMelting(void){
    3940         return NULL;
    3941 }
    3942 /*}}}*/
    3943 /*FUNCTION Penta::CreatePVectorThermal {{{1*/
    3944 ElementVector* Penta::CreatePVectorThermal(void){
    3945 
    3946         /*compute all load vectors for this element*/
    3947         ElementVector* pe1=CreatePVectorThermalVolume();
    3948         ElementVector* pe2=CreatePVectorThermalSheet();
    3949         ElementVector* pe3=CreatePVectorThermalShelf();
    3950         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    3951 
    3952         /*clean-up and return*/
    3953         delete pe1;
    3954         delete pe2;
    3955         delete pe3;
    3956         return pe;
    3957 }
    3958 /*}}}*/
    3959 /*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
    3960 ElementVector* Penta::CreatePVectorThermalVolume(void){
    3961 
    3962         /*Constants*/
    3963         const int  numdof=NUMVERTICES*NDOF1;
    3964 
    3965         /*Intermediaries*/
    3966         int    i,j,ig,found=0;
    3967         int    friction_type,stabilization;
    3968         double Jdet,phi,dt;
    3969         double rho_ice,heatcapacity;
    3970         double thermalconductivity,kappa;
    3971         double viscosity,temperature;
    3972         double tau_parameter,diameter;
    3973         double u,v,w;
    3974         double scalar_def,scalar_transient;
    3975         double temperature_list[NUMVERTICES];
    3976         double xyz_list[NUMVERTICES][3];
    3977         double L[numdof];
    3978         double dbasis[3][6];
    3979         double epsilon[6];
    3980         GaussPenta *gauss=NULL;
    3981 
    3982         /*Initialize Element vector*/
    3983         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    3984 
    3985         /*Retrieve all inputs and parameters*/
    3986         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3987         rho_ice=matpar->GetRhoIce();
    3988         heatcapacity=matpar->GetHeatCapacity();
    3989         thermalconductivity=matpar->GetThermalConductivity();
    3990         kappa=thermalconductivity/(rho_ice*heatcapacity);
    3991         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3992         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3993         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    3994         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    3995         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    3996         Input* temperature_input=NULL;
    3997         if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
    3998         if (stabilization==2) diameter=MinEdgeLength(xyz_list);
    3999 
    4000         /* Start  looping on the number of gaussian points: */
    4001         gauss=new GaussPenta(2,3);
    4002         for (ig=gauss->begin();ig<gauss->end();ig++){
    4003 
    4004                 gauss->GaussPoint(ig);
    4005 
    4006                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4007                 GetNodalFunctionsP1(&L[0], gauss);
    4008 
    4009                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    4010                 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    4011                 GetPhi(&phi, &epsilon[0], viscosity);
    4012 
    4013                 scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
    4014                 if(dt) scalar_def=scalar_def*dt;
    4015 
    4016                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    4017 
    4018                 /* Build transient now */
    4019                 if(dt){
    4020                         temperature_input->GetInputValue(&temperature, gauss);
    4021                         scalar_transient=temperature*Jdet*gauss->weight;
    4022                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
    4023                 }
    4024 
    4025                 if(stabilization==2){
    4026                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    4027 
    4028                         vx_input->GetInputValue(&u, gauss);
    4029                         vy_input->GetInputValue(&v, gauss);
    4030                         vz_input->GetInputValue(&w, gauss);
    4031 
    4032                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
    4033 
    4034                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    4035                         if(dt){
    4036                                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    4037                         }
    4038                 }
    4039         }
    4040 
    4041         /*Clean up and return*/
    4042         delete gauss;
    4043         return pe;
    4044 }
    4045 /*}}}*/
    4046 /*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
    4047 ElementVector* Penta::CreatePVectorThermalShelf(void){
    4048 
    4049         /*Constants*/
    4050         const int  numdof=NUMVERTICES*NDOF1;
    4051 
    4052         /*Intermediaries */
    4053         int        i,j,ig;
    4054         double     Jdet2d;
    4055         double     mixed_layer_capacity,thermal_exchange_velocity;
    4056         double     rho_ice,rho_water,pressure,dt,scalar_ocean;
    4057         double     heatcapacity,t_pmp;
    4058         double     xyz_list[NUMVERTICES][3];
    4059         double     xyz_list_tria[NUMVERTICES2D][3];
    4060         double     basis[NUMVERTICES];
    4061         GaussPenta* gauss=NULL;
    4062 
    4063         /* Ice/ocean heat exchange flux on ice shelf base */
    4064         if (!IsOnBed() || !IsFloating()) return NULL;
    4065 
    4066         /*Initialize Element vector*/
    4067         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4068 
    4069         /*Retrieve all inputs and parameters*/
    4070         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4071         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4072         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    4073         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    4074         rho_water=matpar->GetRhoWater();
    4075         rho_ice=matpar->GetRhoIce();
    4076         heatcapacity=matpar->GetHeatCapacity();
    4077         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4078         Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    4079 
    4080         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4081         gauss=new GaussPenta(0,1,2,2);
    4082         for(ig=gauss->begin();ig<gauss->end();ig++){
    4083 
    4084                 gauss->GaussPoint(ig);
    4085 
    4086                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4087                 GetNodalFunctionsP1(&basis[0], gauss);
    4088 
    4089                 pressure_input->GetInputValue(&pressure,gauss);
    4090                 t_pmp=matpar->TMeltingPoint(pressure);
    4091 
    4092                 scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    4093                 if(dt) scalar_ocean=dt*scalar_ocean;
    4094 
    4095                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    4096         }
    4097 
    4098         /*Clean up and return*/
    4099         delete gauss;
    4100         return pe;
    4101 }
    4102 /*}}}*/
    4103 /*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
    4104 ElementVector* Penta::CreatePVectorThermalSheet(void){
    4105 
    4106         /*Constants*/
    4107         const int  numdof=NUMVERTICES*NDOF1;
    4108 
    4109         /*Intermediaries */
    4110         int        i,j,ig;
    4111         int        analysis_type;
    4112         double     xyz_list[NUMVERTICES][3];
    4113         double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
    4114         double     Jdet2d,dt;
    4115         double     rho_ice,heatcapacity,geothermalflux_value;
    4116         double     basalfriction,alpha2,vx,vy;
    4117         double     basis[NUMVERTICES];
    4118         double     scalar;
    41193823        Friction*  friction=NULL;
    41203824        GaussPenta* gauss=NULL;
    … …  
    41503854                GetNodalFunctionsP1(&basis[0], gauss);
    41513855
    4152                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4153                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4154                         vx_input->GetInputValue(&vx,gauss);
    4155                         vy_input->GetInputValue(&vy,gauss);
    4156                         basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4157 
    4158                         scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    4159                         if(dt) scalar=dt*scalar;
    4160 
    4161                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     3856                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     3857                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3858                vx_input->GetInputValue(&vx,gauss);
     3859                vy_input->GetInputValue(&vy,gauss);
     3860                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     3861               
     3862                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     3863                if(dt) scalar=dt*scalar;
     3864
     3865                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     3866        }
     3867
     3868        /*Clean up and return*/
     3869        delete gauss;
     3870        delete friction;
     3871        return pe;
     3872}
     3873/*}}}*/
     3874/*FUNCTION Penta::CreatePVectorMelting {{{1*/
     3875ElementVector* Penta::CreatePVectorMelting(void){
     3876        return NULL;
     3877}
     3878/*}}}*/
     3879/*FUNCTION Penta::CreatePVectorThermal {{{1*/
     3880ElementVector* Penta::CreatePVectorThermal(void){
     3881
     3882        /*compute all load vectors for this element*/
     3883        ElementVector* pe1=CreatePVectorThermalVolume();
     3884        ElementVector* pe2=CreatePVectorThermalSheet();
     3885        ElementVector* pe3=CreatePVectorThermalShelf();
     3886        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3887
     3888        /*clean-up and return*/
     3889        delete pe1;
     3890        delete pe2;
     3891        delete pe3;
     3892        return pe;
     3893}
     3894/*}}}*/
     3895/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
     3896ElementVector* Penta::CreatePVectorThermalVolume(void){
     3897
     3898        /*Constants*/
     3899        const int  numdof=NUMVERTICES*NDOF1;
     3900
     3901        /*Intermediaries*/
     3902        int    i,j,ig,found=0;
     3903        int    friction_type,stabilization;
     3904        double Jdet,phi,dt;
     3905        double rho_ice,heatcapacity;
     3906        double thermalconductivity;
     3907        double viscosity,temperature;
     3908        double tau_parameter,diameter;
     3909        double u,v,w;
     3910        double scalar_def,scalar_transient;
     3911        double temperature_list[NUMVERTICES];
     3912        double xyz_list[NUMVERTICES][3];
     3913        double L[numdof];
     3914        double dbasis[3][6];
     3915        double epsilon[6];
     3916        GaussPenta *gauss=NULL;
     3917
     3918        /*Initialize Element vector*/
     3919        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3920
     3921        /*Retrieve all inputs and parameters*/
     3922        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3923        rho_ice=matpar->GetRhoIce();
     3924        heatcapacity=matpar->GetHeatCapacity();
     3925        thermalconductivity=matpar->GetThermalConductivity();
     3926        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     3927        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
     3928        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3929        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3930        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     3931        Input* temperature_input=NULL;
     3932        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
     3933        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
     3934
     3935        /* Start  looping on the number of gaussian points: */
     3936        gauss=new GaussPenta(2,3);
     3937        for (ig=gauss->begin();ig<gauss->end();ig++){
     3938
     3939                gauss->GaussPoint(ig);
     3940
     3941                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3942                GetNodalFunctionsP1(&L[0], gauss);
     3943
     3944                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     3945                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3946                GetPhi(&phi, &epsilon[0], viscosity);
     3947
     3948                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
     3949                if(dt) scalar_def=scalar_def*dt;
     3950
     3951                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
     3952
     3953                /* Build transient now */
     3954                if(dt){
     3955                        temperature_input->GetInputValue(&temperature, gauss);
     3956                        scalar_transient=temperature*Jdet*gauss->weight;
     3957                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     3958                }
     3959
     3960                if(stabilization==2){
     3961                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
     3962
     3963                        vx_input->GetInputValue(&u, gauss);
     3964                        vy_input->GetInputValue(&v, gauss);
     3965                        vz_input->GetInputValue(&w, gauss);
     3966
     3967                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
     3968
     3969                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
     3970                        if(dt){
     3971                                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
     3972                        }
     3973                }
     3974        }
     3975
     3976        /*Clean up and return*/
     3977        delete gauss;
     3978        return pe;
     3979}
     3980/*}}}*/
     3981/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
     3982ElementVector* Penta::CreatePVectorThermalShelf(void){
     3983
     3984        /*Constants*/
     3985        const int  numdof=NUMVERTICES*NDOF1;
     3986
     3987        /*Intermediaries */
     3988        int        i,j,ig;
     3989        double     Jdet2d;
     3990        double     mixed_layer_capacity,thermal_exchange_velocity;
     3991        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     3992        double     heatcapacity,t_pmp;
     3993        double     xyz_list[NUMVERTICES][3];
     3994        double     xyz_list_tria[NUMVERTICES2D][3];
     3995        double     basis[NUMVERTICES];
     3996        GaussPenta* gauss=NULL;
     3997
     3998        /* Ice/ocean heat exchange flux on ice shelf base */
     3999        if (!IsOnBed() || !IsFloating()) return NULL;
     4000
     4001        /*Initialize Element vector*/
     4002        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4003
     4004        /*Retrieve all inputs and parameters*/
     4005        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4006        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4007        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     4008        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     4009        rho_water=matpar->GetRhoWater();
     4010        rho_ice=matpar->GetRhoIce();
     4011        heatcapacity=matpar->GetHeatCapacity();
     4012        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4013        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
     4014
     4015        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     4016        gauss=new GaussPenta(0,1,2,2);
     4017        for(ig=gauss->begin();ig<gauss->end();ig++){
     4018
     4019                gauss->GaussPoint(ig);
     4020
     4021                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     4022                GetNodalFunctionsP1(&basis[0], gauss);
     4023
     4024                pressure_input->GetInputValue(&pressure,gauss);
     4025                t_pmp=matpar->TMeltingPoint(pressure);
     4026
     4027                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
     4028                if(dt) scalar_ocean=dt*scalar_ocean;
     4029
     4030                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
     4031        }
     4032
     4033        /*Clean up and return*/
     4034        delete gauss;
     4035        return pe;
     4036}
     4037/*}}}*/
     4038/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
     4039ElementVector* Penta::CreatePVectorThermalSheet(void){
     4040
     4041        /*Constants*/
     4042        const int  numdof=NUMVERTICES*NDOF1;
     4043
     4044        /*Intermediaries */
     4045        int        i,j,ig;
     4046        int        analysis_type;
     4047        double     xyz_list[NUMVERTICES][3];
     4048        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4049        double     Jdet2d,dt;
     4050        double     rho_ice,heatcapacity,geothermalflux_value;
     4051        double     basalfriction,alpha2,vx,vy;
     4052        double     scalar;
     4053        double     basis[NUMVERTICES];
     4054        Friction*  friction=NULL;
     4055        GaussPenta* gauss=NULL;
     4056
     4057        /* Geothermal flux on ice sheet base and basal friction */
     4058        if (!IsOnBed() || IsFloating()) return NULL;
     4059
     4060        /*Initialize Element vector*/
     4061        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4062
     4063        /*Retrieve all inputs and parameters*/
     4064        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4065        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4066        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4067        rho_ice=matpar->GetRhoIce();
     4068        heatcapacity=matpar->GetHeatCapacity();
     4069        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4070        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
     4071        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
     4072        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     4073        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
     4074
     4075        /*Build frictoin element, needed later: */
     4076        friction=new Friction("3d",inputs,matpar,analysis_type);
     4077
     4078        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     4079        gauss=new GaussPenta(0,1,2,2);
     4080        for(ig=gauss->begin();ig<gauss->end();ig++){
     4081
     4082                gauss->GaussPoint(ig);
     4083
     4084                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     4085                GetNodalFunctionsP1(&basis[0], gauss);
     4086
     4087                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4088                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4089                vx_input->GetInputValue(&vx,gauss);
     4090                vy_input->GetInputValue(&vy,gauss);
     4091                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4092               
     4093                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     4094                if(dt) scalar=dt*scalar;
     4095
     4096                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    41624097        }
    41634098
    … …  
    43374272                /*Convert enthalpy into temperature and water fraction*/
    43384273                for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
    4339                 for(i=0;i<numdof;i++){
    4340                         matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
    4341                         if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
    4342                         //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
    4343                 }
    43444274                       
    43454275                this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values));
    … …  
    73897319}
    73907320/*}}}*/
    7391 /*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/
    7392 ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){
    7393 
    7394         int approximation;
    7395         inputs->GetInputValue(&approximation,ApproximationEnum);
    7396 
    7397         switch(approximation){
    7398                 case MacAyealApproximationEnum:
    7399                         return CreateJacobianDiagnosticMacayeal2d();
    7400                 case PattynApproximationEnum:
    7401                         return CreateJacobianDiagnosticPattyn();
    7402                 case StokesApproximationEnum:
    7403                         return CreateJacobianDiagnosticStokes();
    7404                 case NoneApproximationEnum:
    7405                         return NULL;
    7406                 default:
    7407                         _error_("Approximation %s not supported yet",EnumToStringx(approximation));
    7408         }
    7409 }
    7410 /*}}}*/
    7411 /*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{1*/
    7412 ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){
    7413 
    7414         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    7415           bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
    7416           the stiffness matrix. */
    7417         if (!IsOnBed()) return NULL;
    7418 
    7419         /*Depth Averaging B*/
    7420         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    7421 
    7422         /*Call Tria function*/
    7423         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    7424         ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal();
    7425         delete tria->matice; delete tria;
    7426 
    7427         /*Delete B averaged*/
    7428         this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    7429 
    7430         /*clean up and return*/
    7431         return Ke;
    7432 }
    7433 /*}}}*/
    7434 /*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
    7435 ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
    7436 
    7437         /*Constants*/
    7438         const int    numdof=NDOF2*NUMVERTICES;
    7439 
    7440         /*Intermediaries */
    7441         int        i,j,ig;
    7442         double     xyz_list[NUMVERTICES][3];
    7443         double     Jdet;
    7444         double     eps1dotdphii,eps1dotdphij;
    7445         double     eps2dotdphii,eps2dotdphij;
    7446         double     mu_prime;
    7447         double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    7448         double     eps1[3],eps2[3];
    7449         double     phi[NUMVERTICES];
    7450         double     dphi[3][NUMVERTICES];
    7451         GaussPenta *gauss=NULL;
    7452 
    7453         /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
    7454         ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
    7455 
    7456         /*Retrieve all inputs and parameters*/
    7457         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    7458         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    7459         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    7460 
    7461         /* Start  looping on the number of gaussian points: */
    7462         gauss=new GaussPenta(5,5);
    7463         for (ig=gauss->begin();ig<gauss->end();ig++){
    7464 
    7465                 gauss->GaussPoint(ig);
    7466 
    7467                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7468                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    7469 
    7470                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    7471                 matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    7472                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    7473                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    7474                 eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
    7475 
    7476                 for(i=0;i<6;i++){
    7477                         for(j=0;j<6;j++){
    7478                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    7479                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    7480                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    7481                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    7482 
    7483                                 Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    7484                                 Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    7485                                 Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    7486                                 Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    7487                         }
    7488                 }
    7489         }
    7490 
    7491         /*Transform Coordinate System*/
    7492         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
    7493 
    7494         /*Clean up and return*/
    7495         delete gauss;
    7496         return Ke;
    7497 }
    7498 /*}}}*/
    7499 /*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{1*/
    7500 ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
    7501 
    7502         /*Constants*/
    7503         const int    numdof=NDOF4*NUMVERTICES;
    7504 
    7505         /*Intermediaries */
    7506         int        i,j,ig;
    7507         double     xyz_list[NUMVERTICES][3];
    7508         double     Jdet;
    7509         double     eps1dotdphii,eps1dotdphij;
    7510         double     eps2dotdphii,eps2dotdphij;
    7511         double     eps3dotdphii,eps3dotdphij;
    7512         double     mu_prime;
    7513         double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    7514         double     eps1[3],eps2[3],eps3[3];
    7515         double     phi[NUMVERTICES];
    7516         double     dphi[3][NUMVERTICES];
    7517         GaussPenta *gauss=NULL;
    7518 
    7519         /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
    7520         ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
    7521 
    7522         /*Retrieve all inputs and parameters*/
    7523         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    7524         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    7525         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    7526         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    7527 
    7528         /* Start  looping on the number of gaussian points: */
    7529         gauss=new GaussPenta(5,5);
    7530         for (ig=gauss->begin();ig<gauss->end();ig++){
    7531 
    7532                 gauss->GaussPoint(ig);
    7533 
    7534                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7535                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    7536 
    7537                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    7538                 matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    7539                 eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    7540                 eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
    7541                 eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
    7542 
    7543                 for(i=0;i<6;i++){
    7544                         for(j=0;j<6;j++){
    7545                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    7546                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    7547                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    7548                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    7549                                 eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
    7550                                 eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
    7551 
    7552                                 Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    7553                                 Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    7554                                 Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
    7555 
    7556                                 Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    7557                                 Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    7558                                 Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
    7559 
    7560                                 Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
    7561                                 Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
    7562                                 Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
    7563                         }
    7564                 }
    7565         }
    7566 
    7567         /*Transform Coordinate System*/
    7568         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
    7569 
    7570         /*Clean up and return*/
    7571         delete gauss;
    7572         return Ke;
    7573 }
    7574 /*}}}*/
    75757321/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    75767322void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){

  • issm/trunk-jpl-damage/src/c/objects/Elements/Penta.h

    r11422 r11425  
    8888                void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
    8989                void   CreatePVector(Vec pf);
    90                 void   CreateJacobianMatrix(Mat Jff);
    9190                void   DeleteResults(void);
    9291                int    GetNodeIndex(Node* node);
    … …  
    180179                void      GetPhi(double* phi, double*  epsilon, double viscosity);
    181180                void      GetSolutionFromInputsEnthalpy(Vec solutiong);
    182                 double  GetStabilizationParameter(double u, double v, double w, double diameter, double kappa);
     181                double  GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
    183182                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    184183                void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    … …  
    232231                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
    233232                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
    234                 ElementMatrix* CreateJacobianDiagnosticHoriz(void);
    235                 ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void);
    236                 ElementMatrix* CreateJacobianDiagnosticPattyn(void);
    237                 ElementMatrix* CreateJacobianDiagnosticStokes(void);
    238233                void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
    239234                void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);

  • issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp

    r11422 r11425  
    886886        delete gauss;
    887887        return pe;
    888 }
    889 /*}}}*/
    890 /*FUNCTION Tria::CreateJacobianMatrix{{{1*/
    891 void  Tria::CreateJacobianMatrix(Mat Jff){
    892 
    893         /*retrieve parameters: */
    894         ElementMatrix* Ke=NULL;
    895         int analysis_type;
    896         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    897 
    898         /*Checks in debugging {{{2*/
    899         _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
    900         /*}}}*/
    901 
    902         /*Skip if water element*/
    903         if(IsOnWater()) return;
    904 
    905         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    906         switch(analysis_type){
    907 #ifdef _HAVE_DIAGNOSTIC_
    908                 case DiagnosticHorizAnalysisEnum:
    909                         Ke=CreateJacobianDiagnosticMacayeal();
    910                         break;
    911 #endif
    912                 default:
    913                         _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    914         }
    915 
    916         /*Add to global matrix*/
    917         if(Ke){
    918                 Ke->AddToGlobal(Jff);
    919                 delete Ke;
    920         }
    921888}
    922889/*}}}*/
    … …  
    30643031}
    30653032/*}}}*/
    3066 /*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/
    3067 ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
    3068 
    3069         /*Constants*/
    3070         const int    numdof=NDOF2*NUMVERTICES;
    3071 
    3072         /*Intermediaries */
    3073         int        i,j,ig;
    3074         double     xyz_list[NUMVERTICES][3];
    3075         double     Jdet,thickness;
    3076         double     eps1dotdphii,eps1dotdphij;
    3077         double     eps2dotdphii,eps2dotdphij;
    3078         double     mu_prime;
    3079         double     epsilon[3];/* epsilon=[exx,eyy,exy];*/
    3080         double     eps1[2],eps2[2];
    3081         double     phi[NUMVERTICES];
    3082         double     dphi[2][NUMVERTICES];
    3083         GaussTria *gauss=NULL;
    3084 
    3085         /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
    3086         ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
    3087 
    3088         /*Retrieve all inputs and parameters*/
    3089         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3090         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    3091         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    3092         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3093 
    3094         /* Start  looping on the number of gaussian points: */
    3095         gauss=new GaussTria(2);
    3096         for (ig=gauss->begin();ig<gauss->end();ig++){
    3097 
    3098                 gauss->GaussPoint(ig);
    3099 
    3100                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3101                 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    3102 
    3103                 thickness_input->GetInputValue(&thickness, gauss);
    3104                 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    3105                 matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    3106                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    3107                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    3108 
    3109                 for(i=0;i<3;i++){
    3110                         for(j=0;j<3;j++){
    3111                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
    3112                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
    3113                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
    3114                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
    3115 
    3116                                 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
    3117                                 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
    3118                                 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
    3119                                 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
    3120                         }
    3121                 }
    3122         }
    3123 
    3124         /*Transform Coordinate System*/
    3125         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
    3126 
    3127         /*Clean up and return*/
    3128         delete gauss;
    3129         return Ke;
    3130 }
    3131 /*}}}*/
    31323033/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    31333034void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
    … …  
    34313332        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    34323333
    3433         /*If on water, grad = 0: (WHY???) */
    3434         if(IsOnWater()) return; 
     3334        /*If on water, grad = 0: */
     3335        if(IsOnWater()) return;
    34353336
    34363337        /*First deal with ∂/∂alpha(KU-F)*/

  • issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h

    r11422 r11425  
    8484                void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
    8585                void   CreatePVector(Vec pf);
    86                 void   CreateJacobianMatrix(Mat Jff);
    8786                int    GetNodeIndex(Node* node);
    8887                int    Sid();
    … …  
    141140                #ifdef _HAVE_CONTROL_
    142141                double DragCoefficientAbsGradient(bool process_units,int weight_index);
    143                 void   GradientIndexing(int* indexing,int control_index);
    144142                void   Gradj(Vec gradient,int control_type);
    145143                void   GradjBGradient(Vec gradient,int weight_index);
    … …  
    208206                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    209207                ElementVector* CreatePVectorDiagnosticHutter(void);
    210                 ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
    211208                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    212209                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);

  • issm/trunk-jpl-damage/src/c/objects/KML/KML_File.cpp

    r11422 r11425  
    2424KML_File::KML_File(){
    2525
    26         ;
     26        kmlobj    =new DataSet;
    2727
    2828}
    … …  
    3131KML_File::~KML_File(){
    3232
    33         ;
     33        if (kmlobj) {
     34                delete kmlobj;
     35                kmlobj    =NULL;
     36        }
    3437
    3538}
    … …  
    4548        KML_Object::Echo();
    4649
     50        _printf_(flag,"        kmlobj: (size=%d)\n" ,kmlobj->Size());
     51
    4752        return;
    4853}
    … …  
    6166void  KML_File::DeepEcho(const char* indent){
    6267
     68        int   i;
     69        char  indent2[81];
    6370        bool  flag=true;
    6471
    … …  
    6673        KML_Object::DeepEcho(indent);
    6774
     75/*  loop over the kml objects for the file  */
     76
     77        memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
     78
     79        strcat(indent2,"  ");
     80
     81        if (kmlobj->Size())
     82                for (i=0; i<kmlobj->Size(); i++) {
     83                        _printf_(flag,"%s        kmlobj: -------- begin [%d] --------\n" ,indent,i);
     84                        ((KML_Object *)kmlobj->GetObjectByOffset(i))->DeepEcho(indent2);
     85                        _printf_(flag,"%s        kmlobj: --------  end  [%d] --------\n" ,indent,i);
     86                }
     87        else
     88                _printf_(flag,"%s        kmlobj: [empty]\n"    ,indent);
     89
    6890        return;
    6991}
    … …  
    7193/*FUNCTION KML_File::Write {{{1*/
    7294void  KML_File::Write(FILE* filout,const char* indent){
     95
     96        int   i;
     97        char  indent2[81];
    7398
    7499        fprintf(filout,"%s<kml",indent);
    … …  
    78103
    79104        KML_Object::Write(filout,indent);
     105
     106/*  loop over the kml objects for the file  */
     107
     108        memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
     109
     110        strcat(indent2,"  ");
     111
     112        for (i=0; i<kmlobj->Size(); i++)
     113                ((KML_Object *)kmlobj->GetObjectByOffset(i))->Write(filout,indent2);
    80114
    81115        fprintf(filout,"%s</kml>\n",indent);
    … …  
    111145                        _error_("KML_File::Read -- Unexpected field \"%s\".\n",kstri);
    112146
     147                else if (!strncmp(kstri,"<Placemark",10)) {
     148                        kobj=(KML_Object*)new KML_Placemark();
     149                        kobj->Read(fid,kstri);
     150                        kmlobj    ->AddObject((Object*)kobj);
     151                }
     152
     153                else if (!strncmp(kstri,"<Folder", 7)) {
     154                        kobj=(KML_Object*)new KML_Folder();
     155                        kobj->Read(fid,kstri);
     156                        kmlobj    ->AddObject((Object*)kobj);
     157                }
     158
     159                else if (!strncmp(kstri,"<Document", 9)) {
     160                        kobj=(KML_Object*)new KML_Document();
     161                        kobj->Read(fid,kstri);
     162                        kmlobj    ->AddObject((Object*)kobj);
     163                }
     164
     165                else if (!strncmp(kstri,"<GroundOverlay",14)) {
     166                        kobj=(KML_Object*)new KML_GroundOverlay();
     167                        kobj->Read(fid,kstri);
     168                        kmlobj    ->AddObject((Object*)kobj);
     169                }
     170
     171                else if (!strncmp(kstri,"<LatLonBox",10)) {
     172                        kobj=(KML_Object*)new KML_LatLonBox();
     173                        kobj->Read(fid,kstri);
     174                        kmlobj    ->AddObject((Object*)kobj);
     175                }
     176
     177                else if (!strncmp(kstri,"<Icon", 5)) {
     178                        kobj=(KML_Object*)new KML_Icon();
     179                        kobj->Read(fid,kstri);
     180                        kmlobj    ->AddObject((Object*)kobj);
     181                }
     182
     183                else if (!strncmp(kstri,"<Point", 6)) {
     184                        kobj=(KML_Object*)new KML_Point();
     185                        kobj->Read(fid,kstri);
     186                        kmlobj    ->AddObject((Object*)kobj);
     187                }
     188
     189                else if (!strncmp(kstri,"<LineString",11)) {
     190                        kobj=(KML_Object*)new KML_LineString();
     191                        kobj->Read(fid,kstri);
     192                        kmlobj    ->AddObject((Object*)kobj);
     193                }
     194
     195                else if (!strncmp(kstri,"<LinearRing",11)) {
     196                        kobj=(KML_Object*)new KML_LinearRing();
     197                        kobj->Read(fid,kstri);
     198                        kmlobj    ->AddObject((Object*)kobj);
     199                }
     200
     201                else if (!strncmp(kstri,"<Polygon", 8)) {
     202                        kobj=(KML_Object*)new KML_Polygon();
     203                        kobj->Read(fid,kstri);
     204                        kmlobj    ->AddObject((Object*)kobj);
     205                }
     206
     207                else if (!strncmp(kstri,"<MultiGeometry",14)) {
     208                        kobj=(KML_Object*)new KML_MultiGeometry();
     209                        kobj->Read(fid,kstri);
     210                        kmlobj    ->AddObject((Object*)kobj);
     211                }
     212
     213//              else if (!strncmp(kstri,"<IconStyle",10)) {
     214//                      kobj=(KML_Object*)new KML_IconStyle();
     215//                      kobj->Read(fid,kstri);
     216//                      kmlobj    ->AddObject((Object*)kobj);
     217//              }
     218
     219//              else if (!strncmp(kstri,"<LabelStyle",11)) {
     220//                      kobj=(KML_Object*)new KML_LabelStyle();
     221//                      kobj->Read(fid,kstri);
     222//                      kmlobj    ->AddObject((Object*)kobj);
     223//              }
     224
     225                else if (!strncmp(kstri,"<LineStyle",10)) {
     226                        kobj=(KML_Object*)new KML_LineStyle();
     227                        kobj->Read(fid,kstri);
     228                        kmlobj    ->AddObject((Object*)kobj);
     229                }
     230
     231                else if (!strncmp(kstri,"<PolyStyle",10)) {
     232                        kobj=(KML_Object*)new KML_PolyStyle();
     233                        kobj->Read(fid,kstri);
     234                        kmlobj    ->AddObject((Object*)kobj);
     235                }
     236
     237//              else if (!strncmp(kstri,"<BalloonStyle",13)) {
     238//                      kobj=(KML_Object*)new KML_BalloonStyle();
     239//                      kobj->Read(fid,kstri);
     240//                      kmlobj    ->AddObject((Object*)kobj);
     241//              }
     242
     243//              else if (!strncmp(kstri,"<ListStyle",10)) {
     244//                      kobj=(KML_Object*)new KML_ListStyle();
     245//                      kobj->Read(fid,kstri);
     246//                      kmlobj    ->AddObject((Object*)kobj);
     247//              }
     248
    113249                else if (!strncmp(kstri,"<",1))
    114250                        KML_Object::Read(fid,kstri);

  • issm/trunk-jpl-damage/src/c/objects/KML/KML_File.h

    r11417 r11425  
    1919
    2020        public:
     21
     22                DataSet* kmlobj;
    2123
    2224                /*KML_File constructors, destructors {{{1*/

  • issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.cpp

    r11422 r11425  
    2626        attrib    =new DataSet;
    2727        commnt    =new DataSet;
    28         kmlobj    =new DataSet;
    2928
    3029}
    … …  
    4140                commnt    =NULL;
    4241        }
    43         if (kmlobj) {
    44                 delete kmlobj;
    45                 kmlobj    =NULL;
    46         }
    4742
    4843}
    … …  
    5752        _printf_(flag,"        attrib: (size=%d)\n" ,attrib->Size());
    5853        _printf_(flag,"        commnt: (size=%d)\n" ,commnt->Size());
    59         _printf_(flag,"        kmlobj: (size=%d)\n" ,kmlobj->Size());
    6054
    6155        return;
    … …  
    7670
    7771        int   i;
    78         char  indent2[81];
    7972        bool  flag=true;
    8073
    … …  
    9790                _printf_(flag,"%s        commnt: [empty]\n"    ,indent);
    9891
    99 /*  loop over the unknown objects for the object  */
    100 
    101         memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
    102         strcat(indent2,"  ");
    103 
    104         if (kmlobj->Size())
    105                 for (i=0; i<kmlobj->Size(); i++) {
    106             _printf_(flag,"%s        kmlobj: -------- begin [%d] --------\n" ,indent,i);
    107                         ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->DeepEcho(indent2);
    108             _printf_(flag,"%s        kmlobj: --------  end  [%d] --------\n" ,indent,i);
    109                 }
    110         else
    111                 _printf_(flag,"%s        kmlobj: [empty]\n"    ,indent);
    112 
    11392        return;
    11493}
    … …  
    11796void  KML_Object::Write(FILE* filout,const char* indent){
    11897
    119         int   i;
    120         char  indent2[81];
    121 
    12298//  attributes always written in keyword line of derived classes
    12399//  comments always written after keyword line of derived classes
    124100
    125 /*  loop over the unknown objects for the object  */
    126 
    127         memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
    128         strcat(indent2,"  ");
    129 
    130         if (kmlobj->Size())
    131                 for (i=0; i<kmlobj->Size(); i++) {
    132                         ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->Write(filout,indent2);
    133                 }
     101        ;
    134102
    135103        return;
    … …  
    138106/*FUNCTION KML_Object::Read {{{1*/
    139107void  KML_Object::Read(FILE* fid,char* kstr){
    140 
    141         KML_Object*  kobj;
    142108
    143109/*  process field within opening and closing tags  */
    … …  
    150116                _error_("KML_Object::Read -- Unexpected field \"%s\".\n",kstr);
    151117
    152         else if (!strncmp(kstr,"<Placemark",10)) {
    153                 kobj=(KML_Object*)new KML_Placemark();
    154                 kobj->Read(fid,kstr);
    155                 kmlobj    ->AddObject((Object*)kobj);
    156         }
    157 
    158         else if (!strncmp(kstr,"<Folder", 7)) {
    159                 kobj=(KML_Object*)new KML_Folder();
    160                 kobj->Read(fid,kstr);
    161                 kmlobj    ->AddObject((Object*)kobj);
    162         }
    163 
    164         else if (!strncmp(kstr,"<Document", 9)) {
    165                 kobj=(KML_Object*)new KML_Document();
    166                 kobj->Read(fid,kstr);
    167                 kmlobj    ->AddObject((Object*)kobj);
    168         }
    169 
    170         else if (!strncmp(kstr,"<GroundOverlay",14)) {
    171                 kobj=(KML_Object*)new KML_GroundOverlay();
    172                 kobj->Read(fid,kstr);
    173                 kmlobj    ->AddObject((Object*)kobj);
    174         }
    175 
    176         else if (!strncmp(kstr,"<LatLonBox",10)) {
    177                 kobj=(KML_Object*)new KML_LatLonBox();
    178                 kobj->Read(fid,kstr);
    179                 kmlobj    ->AddObject((Object*)kobj);
    180         }
    181 
    182         else if (!strncmp(kstr,"<Icon", 5)) {
    183                 kobj=(KML_Object*)new KML_Icon();
    184                 kobj->Read(fid,kstr);
    185                 kmlobj    ->AddObject((Object*)kobj);
    186         }
    187 
    188         else if (!strncmp(kstr,"<Point", 6)) {
    189                 kobj=(KML_Object*)new KML_Point();
    190                 kobj->Read(fid,kstr);
    191                 kmlobj    ->AddObject((Object*)kobj);
    192         }
    193 
    194         else if (!strncmp(kstr,"<LineString",11)) {
    195                 kobj=(KML_Object*)new KML_LineString();
    196                 kobj->Read(fid,kstr);
    197                 kmlobj    ->AddObject((Object*)kobj);
    198         }
    199 
    200         else if (!strncmp(kstr,"<LinearRing",11)) {
    201                 kobj=(KML_Object*)new KML_LinearRing();
    202                 kobj->Read(fid,kstr);
    203                 kmlobj    ->AddObject((Object*)kobj);
    204         }
    205 
    206         else if (!strncmp(kstr,"<Polygon", 8)) {
    207                 kobj=(KML_Object*)new KML_Polygon();
    208                 kobj->Read(fid,kstr);
    209                 kmlobj    ->AddObject((Object*)kobj);
    210         }
    211 
    212         else if (!strncmp(kstr,"<MultiGeometry",14)) {
    213                 kobj=(KML_Object*)new KML_MultiGeometry();
    214                 kobj->Read(fid,kstr);
    215                 kmlobj    ->AddObject((Object*)kobj);
    216         }
    217 
    218 //      else if (!strncmp(kstr,"<IconStyle",10)) {
    219 //              kobj=(KML_Object*)new KML_IconStyle();
    220 //              kobj->Read(fid,kstr);
    221 //              kmlobj    ->AddObject((Object*)kobj);
    222 //      }
    223 
    224 //      else if (!strncmp(kstr,"<LabelStyle",11)) {
    225 //              kobj=(KML_Object*)new KML_LabelStyle();
    226 //              kobj->Read(fid,kstr);
    227 //              kmlobj    ->AddObject((Object*)kobj);
    228 //      }
    229 
    230         else if (!strncmp(kstr,"<LineStyle",10)) {
    231                 kobj=(KML_Object*)new KML_LineStyle();
    232                 kobj->Read(fid,kstr);
    233                 kmlobj    ->AddObject((Object*)kobj);
    234         }
    235 
    236         else if (!strncmp(kstr,"<PolyStyle",10)) {
    237                 kobj=(KML_Object*)new KML_PolyStyle();
    238                 kobj->Read(fid,kstr);
    239                 kmlobj    ->AddObject((Object*)kobj);
    240         }
    241 
    242 //      else if (!strncmp(kstr,"<BalloonStyle",13)) {
    243 //              kobj=(KML_Object*)new KML_BalloonStyle();
    244 //              kobj->Read(fid,kstr);
    245 //              kmlobj    ->AddObject((Object*)kobj);
    246 //      }
    247 
    248 //      else if (!strncmp(kstr,"<ListStyle",10)) {
    249 //              kobj=(KML_Object*)new KML_ListStyle();
    250 //              kobj->Read(fid,kstr);
    251 //              kmlobj    ->AddObject((Object*)kobj);
    252 //      }
    253 
    254118        else if (!strncmp(kstr,"<",1)) {
    255119                _printf_(true,"KML_Object::Read -- Unrecognized opening tag %s.\n",kstr);
    256 //              KMLFileTagSkip(kstr,
    257 //                                         fid);
    258                 kobj=(KML_Object*)new KML_Unknown();
    259                 kobj->Read(fid,kstr);
    260                 kmlobj    ->AddObject((Object*)kobj);
     120                KMLFileTagSkip(kstr,
     121                                           fid);
    261122        }
    262123

  • issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.h

    r11419 r11425  
    2121                DataSet* attrib;
    2222                DataSet* commnt;
    23                 DataSet* kmlobj;
    2423
    2524                /*KML_Object constructors, destructors {{{1*/

  • issm/trunk-jpl-damage/src/c/objects/KML/KML_Unknown.cpp

    r11422 r11425  
    4444        bool  flag=true;
    4545
    46         _printf_(flag,"KML_Unknown %s:\n",name);
     46        _printf_(flag,"KML_Unknown_%s:\n",name);
    4747        KML_Object::Echo();
    4848
    4949        if (value     )
    5050                _printf_(flag,"         value: \"%s\"\n"     ,value);
    51     else
    52         _printf_(flag,"         value: [none]\n"     );
    5351
    5452        return;
    … …  
    6866void  KML_Unknown::DeepEcho(const char* indent){
    6967
    70         char*        valuei;
    71         char*        vtoken;
    72         char         nl[]={'\n','\0'};
    73         bool         flag=true;
     68        bool  flag=true;
    7469
    75         _printf_(flag,"%sKML_Unknown %s:\n",indent,name);
     70        _printf_(flag,"%sKML_Unknown_%s:\n",indent,name);
    7671        KML_Object::DeepEcho(indent);
    7772
    78         if (value     ) {
    79                 valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char));
    80                 memcpy(valuei,value,(strlen(value)+1)*sizeof(char));
    81        
    82                 vtoken=strtok(valuei,nl);
    83                 _printf_(flag,"%s         value: \"%s"     ,indent,vtoken);
    84    
    85                 while (vtoken=strtok(NULL,nl))
    86                         _printf_(flag,"\n%s                 %s"     ,indent,vtoken);
    87                 _printf_(flag,"\"\n");
    88 
    89                 xfree((void**)&valuei);
    90         }
    91     else
    92         _printf_(flag,"%s         value: [none]\n"     ,indent);
     73        if (value     )
     74                _printf_(flag,"%s         value: \"%s\"\n"     ,indent,value);
    9375
    9476        return;
    … …  
    9880void  KML_Unknown::Write(FILE* filout,const char* indent){
    9981
    100         char*        valuei;
    101         char*        vtoken;
    102         char         nl[]={'\n','\0'};
    103 
    10482        fprintf(filout,"%s<%s",indent,name);
    10583        WriteAttrib(filout," ");
    … …  
    10785        WriteCommnt(filout,indent);
    10886
    109         if (value     ) {
    110                 valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char));
    111                 memcpy(valuei,value,(strlen(value)+1)*sizeof(char));
    112        
    113                 vtoken=strtok(valuei,nl);
    114                 fprintf(filout,"%s  %s\n",indent,vtoken);
    115    
    116                 while (vtoken=strtok(NULL,nl))
    117                         fprintf(filout,"%s  %s\n",indent,vtoken);
     87        KML_Object::Write(filout,indent);
    11888
    119                 xfree((void**)&valuei);
    120         }
    121 
    122         KML_Object::Write(filout,indent);
     89        if (value     )
     90                fprintf(filout,"%s  %s\n",indent,value);
    12391
    12492        fprintf(filout,"%s</%s>\n",indent,name);
    … …  
    133101        int          ncom=0;
    134102        char**       pcom=NULL;
    135         char         nl[]={'\n','\0'};
    136 
    137 /*  get object name  */
    138 
    139         /*name=KMLFileTagName(NULL,
    140                                                 kstr);*/
    141 //      _printf_(true,"KML_Unknown::Read -- opening name=%s.\n",name);
    142103
    143104/*  get object attributes and check for solo tag  */
    … …  
    151112        while (kstri=KMLFileToken(fid,
    152113                                                          &ncom,&pcom)) {
    153 //              _printf_(true,"KML_Unknown::Read -- kstri=%s.\n",kstri);
    154114                if      (!strncmp(&kstri[0],"</", 2) &&
    155115                                 !strncmp(&kstri[2],name,strlen(name))) {
    156 //                      _printf_(true,"KML_Unknown::Read -- closing name=%s.\n",name);
    157116                        xfree((void**)&kstri);
    158117                        break;
    … …  
    161120                        _error_("KML_Unknown::Read -- Unexpected closing tag %s.\n",kstri);
    162121
    163                 else if (strncmp(kstri,"<",1)) {
    164                         if (value) {
    165                                 value=(char *) xrealloc(value,(strlen(value)+1+strlen(kstri)+1)*sizeof(char));
    166                                 strcat(value,nl);
    167                                 strcat(value,kstri);
    168                         }
    169                         else {
    170                                 value=(char *) xmalloc((strlen(kstri)+1)*sizeof(char));
    171                                 memcpy(value,kstri,(strlen(kstri)+1)*sizeof(char));
    172                         }
    173                 }
     122                else if (strncmp(kstri,"<",1))
     123                        KMLFileTokenParse( value     ,NULL,0,
     124                                                          NULL,
     125                                                          fid);
    174126
    175127                else if (!strncmp(kstri,"<",1))

  • issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.cpp

    r11417 r11425  
    361361}
    362362/*}}}*/
    363 /*FUNCTION Icefront::CreateJacobianMatrix{{{1*/
    364 void  Icefront::CreateJacobianMatrix(Mat Jff){
    365         this->CreateKMatrix(Jff,NULL);
    366 }
    367 /*}}}1*/
    368363/*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/
    369364void  Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){
    … …  
    378373}
    379374/*}}}*/
    380 /*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/
    381 void  Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
    382         this->PenaltyCreateKMatrix(Jff,NULL,kmax);
    383 }
    384 /*}}}1*/
    385375/*FUNCTION Icefront::InAnalysis{{{1*/
    386376bool Icefront::InAnalysis(int in_analysis_type){

  • issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.h

    r11419 r11425  
    7676                void  CreateKMatrix(Mat Kff, Mat Kfs);
    7777                void  CreatePVector(Vec pf);
    78                 void  CreateJacobianMatrix(Mat Jff);
    7978                void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
    8079                void  PenaltyCreatePVector(Vec pf, double kmax);
    81                 void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
    8280                bool  InAnalysis(int analysis_type);
    8381                /*}}}*/

  • issm/trunk-jpl-damage/src/c/objects/Loads/Penpair.cpp

    r11419 r11425  
    212212        return;
    213213
    214 }
    215 /*}}}1*/
    216 /*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
    217 void  Penpair::CreateJacobianMatrix(Mat Jff){
    218         this->CreateKMatrix(Jff,NULL);
    219214}
    220215/*}}}1*/

  • issm/trunk-jpl-damage/src/c/objects/Loads/Penpair.h

    r11422 r11425  
    6464                void  CreateKMatrix(Mat Kff, Mat Kfs);
    6565                void  CreatePVector(Vec pf);
    66                 void  CreateJacobianMatrix(Mat Jff);
    67                 void  PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax);
     66                void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
    6867                void  PenaltyCreatePVector(Vec pf, double kmax);
    69                 void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
    7068                bool  InAnalysis(int analysis_type);
    7169                /*}}}*/

  • issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp

    r11422 r11425  
    606606        /*Return: */
    607607        *pviscosity_complement=viscosity_complement;
    608 }
    609 /*}}}*/
    610 /*FUNCTION Matice::GetViscosityDerivativeEpsSquare{{{1*/
    611 void  Matice::GetViscosityDerivativeEpsSquare(double* pmu_prime, double* epsilon){
    612 
    613         /*output: */
    614         double mu_prime;
    615         double mu,n,eff2;
    616 
    617         /*input strain rate: */
    618         double exx,eyy,exy,exz;
    619 
    620         /*Get visocisty and n*/
    621         GetViscosity2d(&mu,epsilon);
    622         n=GetN();
    623 
    624         if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    625                 mu_prime=0.5*pow((double)10,(double)14);
    626         }
    627         else{
    628                 /*Retrive strain rate components: */
    629                 exx=epsilon[0];
    630                 eyy=epsilon[1];
    631                 exy=epsilon[2];
    632                 eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ;
    633 
    634                 mu_prime=(1-n)/(2*n) * mu/eff2;
    635         }
    636 
    637         /*Assign output pointers:*/
    638         *pmu_prime=mu_prime;
    639608}
    640609/*}}}*/

  • issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h

    r11422 r11425  
    6969                void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
    7070                void   GetViscosityZComplement(double* pviscosity_complement, double* pepsilon);
    71                 void   GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
    7271                double GetB();
    7372                double GetBbar();

  • issm/trunk-jpl-damage/src/c/objects/Numerics/ElementMatrix.cpp

    r11422 r11425  
    251251        double* localvalues=NULL;
    252252
    253         /*If Kfs is not provided, call the other function*/
    254         if(!Kfs){
    255                 this->AddToGlobal(Kff);
    256                 return;
    257         }
    258 
    259253        /*In debugging mode, check consistency (no NaN, and values not too big)*/
    260254        this->CheckConsistency();
    … …  
    300294}
    301295/*}}}*/
    302 /*FUNCTION ElementMatrix::AddToGlobal(Mat Jff){{{1*/
    303 void ElementMatrix::AddToGlobal(Mat Jff){
    304          
    305         int i,j;
    306         double* localvalues=NULL;
    307        
    308         /*Check that Jff is not NULL*/
    309         _assert_(Jff);
    310        
    311         /*In debugging mode, check consistency (no NaN, and values not too big)*/
    312         this->CheckConsistency();
    313        
    314         if(this->dofsymmetrical){
    315                 /*only use row dofs to add values into global matrices: */
    316                
    317                 if(this->row_fsize){
    318                         /*first, retrieve values that are in the f-set from the g-set values matrix: */
    319                         localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
    320                         for(i=0;i<this->row_fsize;i++){
    321                                 for(j=0;j<this->row_fsize;j++){
    322                                         *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
    323                                 }
    324                         }
    325                         /*add local values into global  matrix, using the fglobaldoflist: */
    326                         MatSetValues(Jff,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES);
    327                        
    328                         /*Free ressources:*/
    329                         xfree((void**)&localvalues);
    330                 }
    331                
    332         }
    333         else{
    334                 _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
    335         }
    336        
    337 }
    338 /*}}}*/
    339296/*FUNCTION ElementMatrix::CheckConsistency{{{1*/
    340297void ElementMatrix::CheckConsistency(void){

  • issm/trunk-jpl-damage/src/c/objects/Numerics/ElementMatrix.h

    r11422 r11425  
    5959                /*ElementMatrix specific routines {{{1*/
    6060                void AddToGlobal(Mat Kff, Mat Kfs);
    61                 void AddToGlobal(Mat Jff);
    6261                void Echo(void);
    6362                void CheckConsistency(void);

  • issm/trunk-jpl-damage/src/c/objects/Params/BoolParam.h

    r11417 r11425  
    5757                void  GetParameterValue(double* pdouble){_error_("Bool param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5858                void  GetParameterValue(char** pstring){_error_("Bool param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    59                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("Bool param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     59                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Bool param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6060                void  GetParameterValue(double** pdoublearray,int* pM){_error_("Bool param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("Bool param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatArrayParam.h

    r11417 r11425  
    6060                void  GetParameterValue(double* pdouble){_error_("DoubleMatArray param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(char** pstring){_error_("DoubleMatArray param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    62                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     62                void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6363                void  GetParameterValue(double** pdoublearray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6464                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("DoubleMatArray param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatParam.h

    r11417 r11425  
    5959                void  GetParameterValue(double* pdouble){_error_("DoubleMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    6060                void  GetParameterValue(char** pstring){_error_("DoubleMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    61                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     61                void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6363                void  GetParameterValue(double** pdoublearray,int* pM,int* pN);

  • issm/trunk-jpl-damage/src/c/objects/Params/DoubleParam.h

    r11417 r11425  
    5858                void  GetParameterValue(double* pdouble){*pdouble=value;}
    5959                void  GetParameterValue(char** pstring){_error_("Double param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("Double param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Double param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM);
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN);

  • issm/trunk-jpl-damage/src/c/objects/Params/DoubleVecParam.h

    r11417 r11425  
    5555                void  GetParameterValue(int* pinteger){_error_("DoubleVec param of enum %i (%s) cannot return an integer",enum_type,EnumToStringx(enum_type));}
    5656                void  GetParameterValue(int** pintarray,int* pM);
    57                 void  GetParameterValue(int** pintarray,int* pM,int* pN){_error_("DoubleVec param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));};
     57                void  GetParameterValue(int** pintarray,int* pM,int* pN){_error_("DoubleVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));};
    5858                void  GetParameterValue(double* pdouble){_error_("DoubleVec param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5959                void  GetParameterValue(char** pstring){_error_("DoubleVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleVec param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleVec param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM);
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN);

  • issm/trunk-jpl-damage/src/c/objects/Params/FileParam.h

    r11417 r11425  
    5757                void  GetParameterValue(double* pdouble){_error_("FileParam of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5858                void  GetParameterValue(char** pstring){_error_("FileParam of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    59                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("FileParam of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     59                void  GetParameterValue(char*** pstringarray,int* pM){_error_("FileParam of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6060                void  GetParameterValue(double** pdoublearray,int* pM){_error_("FileParam of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("FileParam of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/IntMatParam.h

    r11417 r11425  
    5959                void  GetParameterValue(double* pdouble){_error_("IntMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    6060                void  GetParameterValue(char** pstring){_error_("IntMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    61                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     61                void  GetParameterValue(char*** pstringarray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6363                void  GetParameterValue(double** pdoublearray,int* pM,int* pN){_error_("IntMat param of enum %i (%s) cannot return a matrix array",enum_type,EnumToStringx(enum_type));};

  • issm/trunk-jpl-damage/src/c/objects/Params/IntParam.h

    r11417 r11425  
    5858                void  GetParameterValue(double* pdouble){_error_("Int param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5959                void  GetParameterValue(char** pstring){_error_("Int param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("Int param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Int param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM){_error_("Int param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("Int param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/PetscMatParam.h

    r11417 r11425  
    5858                void  GetParameterValue(double* pdouble){_error_("PetscMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5959                void  GetParameterValue(char** pstring){_error_("PetscMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("PetscMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/PetscVecParam.h

    r11417 r11425  
    5858                void  GetParameterValue(double* pdouble){_error_("PetscVec param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5959                void  GetParameterValue(char** pstring){_error_("PetscVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("PetscVec param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Params/StringParam.h

    r11417 r11425  
    5858                void  GetParameterValue(double* pdouble){_error_("String param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5959                void  GetParameterValue(char** pstring);
    60                 void  GetParameterValue(char*** pstringarray,int* pM){_error_("String param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("String param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));}
    6161                void  GetParameterValue(double** pdoublearray,int* pM){_error_("String param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}
    6262                void  GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("String param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));}

  • issm/trunk-jpl-damage/src/c/objects/Patch.cpp

    r11417 r11425  
    1212#include <stdio.h>
    1313#include <string.h>
    14 #include <math.h>
    1514#include "./objects.h"
    1615#include "../Container/Container.h"

  • issm/trunk-jpl-damage/src/c/shared/Numerics/Synchronize.sh

    r11417 r11425  
    165165
    166166        verbositylevel = (int)level;
    167         return verbositylevel;
    168167
    169168#else

  • issm/trunk-jpl-damage/src/c/shared/Numerics/Verbosity.cpp

    r11417 r11425  
    6464
    6565        verbositylevel = (int)level;
    66         return verbositylevel;
    6766
    6867#else

  • issm/trunk-jpl-damage/src/c/shared/TriMesh/TriMeshUtils.cpp

    r11417 r11425  
    124124******************************************************************************************************************************/
    125125
    126 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
     126int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
    127127       
    128128        int i,counter;
    … …  
    379379                                   IsInPoly
    380380******************************************************************************************************************************/
    381 //void IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){
    382 //
    383 //      int i;
    384 //      double x0,y0;
    385 //
    386 //      /*Go through all nodes of the mesh:*/
    387 //      for (i=0;i<nods;i++){
    388 //              if (in[i]){
    389 //                      /*this node already is inside one of the contours, continue*/
    390 //                      continue;
    391 //              }
    392 //              /*pick up node: */
    393 //              x0=x[i];
    394 //              y0=y[i];
    395 //              if (pnpoly(numnodes,xc,yc,x0,y0)){
    396 //                      in[i]=1;
    397 //              }
    398 //      }
    399 //}
     381int IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){
     382
     383        int i;
     384        double x0,y0;
     385
     386        /*Go through all nodes of the mesh:*/
     387        for (i=0;i<nods;i++){
     388                if (in[i]){
     389                        /*this node already is inside one of the contours, continue*/
     390                        continue;
     391                }
     392                /*pick up node: */
     393                x0=x[i];
     394                y0=y[i];
     395                if (pnpoly(numnodes,xc,yc,x0,y0)){
     396                        in[i]=1;
     397                }
     398        }
     399}
    400400
    401401/******************************************************************************************************************************

  • issm/trunk-jpl-damage/src/c/shared/TriMesh/trimesh.h

    r11417 r11425  
    2626int IsNeighbor(int el1,int el2,double* index);
    2727int IsOnRift(int el,int nriftsegs,int* riftsegments);
    28 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments);
     28int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments);
    2929int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel);
    3030int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs);

  • issm/trunk-jpl-damage/src/c/solutions/AnalysisConfiguration.cpp

    r11417 r11425  
    9595
    9696                case TransientSolutionEnum:
    97                         numanalyses=9;
     97                        numanalyses=8;
    9898                        analyses=(int*)xmalloc(numanalyses*sizeof(int));
    9999                        analyses[0]=DiagnosticHorizAnalysisEnum;
    … …  
    104104                        analyses[5]=ThermalAnalysisEnum;
    105105                        analyses[6]=MeltingAnalysisEnum;
    106                         analyses[7]=EnthalpyAnalysisEnum;
    107                         analyses[8]=PrognosticAnalysisEnum;
     106                        analyses[7]=PrognosticAnalysisEnum;
    108107                        break;
    109108               

  • issm/trunk-jpl-damage/src/c/solutions/controltao_core.cpp

    r11419 r11425  
    3131        Vec       XL = NULL;
    3232        Vec       XU = NULL;
    33         int        ierr;
    34         int        num_controls,solution_type;
    35         int        nsteps,maxiter;
    36         AppCtx     user;
    37         TaoSolver  tao;
    38         double    *dummy          = NULL;
    39         int       *control_list   = NULL;
    40         Vec        X              = NULL;
    41         Vec        XL             = NULL;
    42         Vec        XU             = NULL;
    4333
    4434        /*Initialize TAO*/
    … …  
    4737        ierr = TaoInitialize(&argc,&args,(char*)0,"");
    4838        if(ierr) _error_("Could not initialize Tao");
    49 
    50         /*Recover some parameters*/
    51         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    52         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    53         femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum);
    54         femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
    55         femmodel->parameters->FindParam(&dummy,NULL,NULL,InversionMaxiterPerStepEnum);
    56         maxiter=nsteps*(int)dummy[0]; xfree((void**)&dummy);
    5739
    5840        /*Initialize TAO*/
    … …  
    6345
    6446        /*Prepare all TAO parameters*/
    65         TaoSetMaximumFunctionEvaluations(tao,maxiter);
    66         TaoSetMaximumIterations(tao,nsteps);
     47        TaoSetMaximumFunctionEvaluations(tao,50);
     48        TaoSetMaximumIterations(tao,10);
    6749        TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28);
    6850
    … …  
    8466        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
    8567        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
    86 
    87         /*Finalize*/
    88         _printf_(VerboseControl(),"%s\n","   preparing final solution");
    89         femmodel->parameters->SetParam(false,InversionIscontrolEnum); //needed to turn control result output in solutioncore
    90         void (*solutioncore)(FemModel*)=NULL;
    91         CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
    92         solutioncore(femmodel);
    9368
    9469        /*Clean up and return*/
    … …  
    11590        VecFree(&gradient);
    11691        CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    117 
    118         /*Compute gradient*/
    119         Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    120         VecCopy(gradient,G); VecFree(&gradient);
    121         VecScale(G,-1.);
    122 
    123         /*Clean-up and return*/
    124         xfree((void**)&cost_functions);
    125         xfree((void**)&cost_functionsd);
    12692        return 0;
    12793}

  • issm/trunk-jpl-damage/src/c/solutions/transient_core.cpp

    r11417 r11425  
    2424        /*parameters: */
    2525        double finaltime,dt,yts;
    26         bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy;
     26        bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline;
    2727        bool   dakota_analysis=false;
    2828        bool   time_adapt=false;
    … …  
    5151        femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
    5252        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    53         femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    5453        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5554        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    … …  
    9291                        _printf_(VerboseSolution(),"   computing temperatures:\n");
    9392                        #ifdef _HAVE_THERMAL_
    94                         if(isenthalpy==0){
    95                                 thermal_core_step(femmodel,step,time);
    96                         }
    97                         else{
    98                                 enthalpy_core_step(femmodel,step,time);
    99                         }
     93                        thermal_core_step(femmodel,step,time);
    10094                        #else
    10195                        _error_("ISSM was not compiled with thermal capabilities. Exiting");
    … …  
    141135                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
    142136                        if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
    143                         if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum,step,time);
    144                         if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,step,time);
    145                         if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time);
     137                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time);
    146138                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time);
    147139                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum,step,time);

  • issm/trunk-jpl-damage/src/c/solvers/solver_nonlinear.cpp

    r11417 r11425  
    8585                if(count>=max_nonlinear_iterations){
    8686                        _printf_(true,"   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
    87                         converged=true;
    88                 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
    89                 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    9087                        break;
    9188                }

  • issm/trunk-jpl-damage/src/m/classes/clusters/castor.m

    r11417 r11425  
    1010                 % {{{1
    1111                 name='castor'
    12                  login='username';
    13                  np   =128;
     12                 login='larour';
     13                 np   =128; %number of processors
    1414                 port=0;
    1515                 queue='shortc';
    1616                 time=180;
    17                  codepath='/workp/edw/issm-2.0/bin'
    18                  executionpath='/workp/edw/Testing/Execution'
     17                 codepath='/workp/edw/larour/issm-2.0/bin'
     18                 executionpath='/workp/edw/larour/Testing/Execution'
    1919                 %}}}
    2020         end

  • issm/trunk-jpl-damage/src/m/classes/clusters/cosmos.m

    r11417 r11425  
    1010                 % {{{1
    1111                 name='cosmos'
    12                  login='username';
     12                 login='larour';
    1313                 np=128;
    1414                 port=0;
    1515                 queue='shortq';
    1616                 time=3*60;
    17                  codepath='/work00/edw/issm-2.0/bin';
    18                  executionpath='/work00/edw/Execution';
     17                 codepath='/work00/edw/larour/issm-2.0/bin';
     18                 executionpath='/work00/edw/larour/Execution';
    1919                 %}}}
    2020         end

  • issm/trunk-jpl-damage/src/m/classes/clusters/gemini.m

    r11417 r11425  
    1010        % {{{1
    1111                name='gemini'
    12                 login='username';
     12                login='larour';
    1313                np=50;
    1414                port=0;
    1515                queue='debug';
    1616                time=60;
    17                 codepath='/workg/edw/issm-2.0/bin'
    18                 executionpath='/workg/edw/Testing/Execution'
     17                codepath='/workg/edw/larour/issm-2.0/bin'
     18                executionpath='/workg/edw/larour/Testing/Execution'
    1919        %}}}
    2020    end

  • issm/trunk-jpl-damage/src/m/classes/clusters/pollux.m

    r11417 r11425  
    1010                 % {{{1
    1111                 name='pollux'
    12                  login='username';
     12                 login='larour';
    1313                 np=128;
    1414                 port=0;
    1515                 queue='shortp';
    1616                 time=180;
    17                  codepath='/workc/edw/issm-2.0/bin'
    18                  executionpath='/workc/edw/Testing/Execution'
     17                 codepath='/workc/edw/larour/issm-2.0/bin'
     18                 executionpath='/workc/edw/larour/Testing/Execution'
    1919                 %}}}
    2020         end

  • issm/trunk-jpl-damage/src/m/classes/initialization.m

    r11417 r11425  
    6262                                checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
    6363                        end
    64                         if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy),
     64                        if ismember(EnthalpyAnalysisEnum,analyses),
    6565                                checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
    6666                        end

  • issm/trunk-jpl-damage/src/m/classes/solver.m

    r11417 r11425  
    66classdef solver
    77    properties (SetAccess=public)
    8                  options=cell(0,0);
     8                 options={NoneAnalysisEnum,mumpsoptions};
    99         end
    1010         methods
    … …  
    2626                 function obj = setdefaultparameters(obj) % {{{
    2727
    28                          %MUMPS is the default solver
    29                          obj.options={'NoneAnalysis',mumpsoptions};
    30 
    3128                 end % }}}
    32                  function obj = addoptions(obj,analysis,solveroptions) % {{{1
    33 
    34                          %Convert analysis from enum to string
    35                          analysis=EnumToString(analysis);
    36 
     29                 function obj=addoptions(obj,analysis,solveroptions) % {{{1
    3730                         %first, find out if analysis has already been supplied
    3831                         found=false;
    3932                         for i=1:size(obj.options,1),
    4033                                 inanalysis=obj.options{i,1};
    41                                  if strcmp(inanalysis,analysis),
     34                                 if inanalysis==analysis,
    4235                                         found=true;
    43                                          obj.options{i,1} = analysis;
    44                                          obj.options{i,2} = solveroptions;
     36                                         obj.options{i,1}=analysis;
     37                                         obj.options{i,2}=solveroptions;
    4538                                         break;
    4639                                 end
    4740                         end
    48 
    4941                         if ~found,
    50                                  obj.options{end+1,1}= analysis;
    51                                  obj.options{end,2}  = solveroptions;
     42                                 obj.options{end+1,1}=analysis;
     43                                 obj.options{end,2}=solveroptions;
    5244                         end
    5345                 end
    5446                 %}}}
    5547                 function checkconsistency(obj,md,solution,analyses) % {{{
    56                          for i=1:size(obj.options,1),
    57                                  if ~ischar(obj.options{i,1}),
    58                                          checkmessage('solver is not well formatted: Analyses are not strings');
    59                                  end
    60                          end
     48
    6149                 end % }}}
    6250                 function PetscFile(solver,filename) % {{{
    … …  
    8270
    8371                                 %first write analysis:
    84                                  fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
     72                                 fprintf(fid,'\n+%s\n',EnumToString(analysis)); %append a + to recognize it's an analysis enum
    8573
    8674                                 %now, write options
    … …  
    138126                                end
    139127
    140                                 disp(sprintf('   %s -> ''%s''',analysis,string));
     128                                disp(sprintf('   %s -> ''%s''',EnumToString(analysis),string));
    141129                        end
    142130                 end

  • issm/trunk-jpl-damage/src/m/classes/thermal.m

    r11417 r11425  
    1212                penalty_lock      = 0;
    1313                penalty_factor    = 0;
    14                 isenthalpy        = 0;
    1514        end
    1615        methods
    … …  
    4342                        %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
    4443                        obj.penalty_factor=3;
    45 
    46                         %Should we use cold ice (default) or enthalpy formulation
    47                         obj.isenthalpy=0;
    4844                end % }}}
    4945                function checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5450                        checkfield(md,'thermal.stabilization','numel',1,'values',[0 1 2]);
    5551                        checkfield(md,'thermal.spctemperature','forcing',1);
    56                         if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy),
     52                        if ismember(EnthalpyAnalysisEnum,analyses),
    5753                        checkfield(md,'thermal.spctemperature','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*md.geometry.thickness,'message','spctemperature should be below the adjusted melting point');
    58                         checkfield(md,'thermal.isenthalpy','numel',1,'values',[0 1]);
    5954                        end
    6055                end % }}}
    … …  
    6762                        fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
    6863                        fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
    69                         fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    7064
    7165                end % }}}
    … …  
    7771                        WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
    7872                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
    79                         WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
    8073                end % }}}
    8174        end

  • issm/trunk-jpl-damage/src/m/model/collapse.m

    r11417 r11425  
    2929if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
    3030if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
    31 if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
    32 if ~isnan(md.inversion.min_parameters), md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
    33 if ~isnan(md.inversion.max_parameters), md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
    3431if ~isnan(md.surfaceforcings.mass_balance),
    3532        md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers);
    … …  
    5249if ~isnan(md.flowequation.element_equation)
    5350        md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
    54         md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
    55         md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1);
    56         md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1);
    57         md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1);
    5851end     
    5952
    … …  
    6255md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers);
    6356md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers);
    64 md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers);
    6557md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers);
    6658md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
    … …  
    6860%Extrusion of Neumann BC
    6961if ~isnan(md.diagnostic.icefront),
    70         numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1);
     62        numberofneumann2d=size(md.diagnostic.icefront,1)/md.mesh.numberoflayers;
    7163        md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
    7264end
    … …  
    9789md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1);
    9890md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1);
    99 md.mask.elementonwater=project2d(md,md.mask.elementonwater,1);
    100 md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1);
    10191
    10292%Initialize with the 2d mesh

  • issm/trunk-jpl-damage/src/m/model/loadresultsfromdisk.m

    r11417 r11425  
    1212
    1313if ~md.qmu.isdakota,
    14 
    15         %Check that file exists
    16         if ~exist(filename,'file'),
    17                 error(['binary file ' filename ' not found.']);
    18         end
    1914
    2015        %initialize md.results if not a structure yet

  • issm/trunk-jpl-damage/src/m/model/radarpower.m

    r11417 r11425  
    1111
    1212%If gdal does not work, uncomment the following line
    13 %setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
     13setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
    1414%Parse inputs
    1515if nargin==1,

  • issm/trunk-jpl-damage/src/m/model/tres.m

    r11417 r11425  
    1010if strcmpi(string,'diagnostic'),
    1111        if md.mesh.dimension==2,
    12                 md.initialization.vx=md.results.DiagnosticSolution.Vx;
    13                 md.initialization.vy=md.results.DiagnosticSolution.Vy;
     12                if isfield(md.results.DiagnosticSolution,'VxAverage'),
     13                        md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.VxAverage);
     14                else
     15                        md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.Vx);
     16                end
     17                if isfield(md.results.DiagnosticSolution,'VyAverage'),
     18                        md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.VyAverage);
     19                else
     20                        md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.Vy);
     21                end
    1422        else
    15                 md.initialization.vx=md.results.DiagnosticSolution.Vx;
    16                 md.initialization.vy=md.results.DiagnosticSolution.Vy;
    17                 md.initialization.vz=md.results.DiagnosticSolution.Vz;
     23                md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.Vx);
     24                md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.Vy);
     25                if isfield(md.results.DiagnosticSolution,'Vz'),
     26                        md.initialization.vz=PatchToVec(md.results.DiagnosticSolution.Vz);
     27                else
     28                        md.initialization.vz=zeros(md.mesh.numberofvertices,1);
     29                end
    1830        end
    19         md.initialization.vel=md.results.DiagnosticSolution.Vel;
     31        md.initialization.vel=PatchToVec(md.results.DiagnosticSolution.Vel);
    2032
    2133        if isfield(md.results.DiagnosticSolution,'Pressure'),
    22                 md.initialization.pressure=md.results.DiagnosticSolution.Pressure;
     34                md.initialization.pressure=PatchToVec(md.results.DiagnosticSolution.Pressure);
    2335        end
    2436        if md.rifts.numrifts,
    … …  
    3042                for control_parameters=md.inversion.control_parameters
    3143                        %Will need to be updated... good luck ;)
    32                         md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters));
     44                        md.(EnumToModelField(control_parameters))=PatchToVec(md.results.DiagnosticSolution.(EnumToString(control_parameters)));
    3345                end
    3446        end
    … …  
    4759        for i=1:length(results),
    4860                if ~isempty(md.results.TransientSolution(i).Vel),
    49                         results2(count).Vel=md.results.TransientSolution(i).Vel;
    50                         results2(count).Surface=md.results.TransientSolution(i).Surface;
    51                         results2(count).Thickness=md.results.TransientSolution(i).Thickness;
    52                         results2(count).Bed=md.results.TransientSolution(i).Bed;
    53                         results2(count).Vx=md.results.TransientSolution(i).Vx;
    54                         results2(count).Vy=md.results.TransientSolution(i).Vy;
     61                        results2(count).Vel=PatchToVec(md.results.TransientSolution(i).Vel);
     62                        results2(count).Surface=PatchToVec(md.results.TransientSolution(i).Surface);
     63                        results2(count).Thickness=PatchToVec(md.results.TransientSolution(i).Thickness);
     64                        results2(count).Bed=PatchToVec(md.results.TransientSolution(i).Bed);
     65                        results2(count).Vx=PatchToVec(md.results.TransientSolution(i).Vx);
     66                        results2(count).Vy=PatchToVec(md.results.TransientSolution(i).Vy);
    5567                        results2(count).time=md.results.TransientSolution(i).time;
    5668                        results2(count).step=md.results.TransientSolution(i).step;
    … …  
    6476        clear results,results2;
    6577elseif strcmpi(string,'steadystate'),
    66         md.initialization.vx=md.results.SteadystateSolution.Vx;
    67         md.initialization.vy=md.results.SteadystateSolution.Vy;
     78        md.initialization.vx=PatchToVec(md.results.SteadystateSolution.Vx);
     79        md.initialization.vy=PatchToVec(md.results.SteadystateSolution.Vy);
    6880        if isfield(md.results.SteadystateSolution,'Vz'),
    69                 md.initialization.vz=md.results.SteadystateSolution.Vz;
     81                md.initialization.vz=PatchToVec(md.results.SteadystateSolution.Vz);
    7082        end
    7183
    72         md.initialization.vel=md.results.SteadystateSolution.Vel;
    73         md.initialization.pressure=md.results.SteadystateSolution.Pressure;
    74         md.initialization.temperature=md.results.SteadystateSolution.Temperature;
    75         md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate;
     84        md.initialization.vel=PatchToVec(md.results.SteadystateSolution.Vel);
     85        md.initialization.pressure=PatchToVec(md.results.SteadystateSolution.Pressure);
     86        md.initialization.temperature=PatchToVec(md.results.SteadystateSolution.Temperature);
     87        md.basalforcings.melting_rate=PatchToVec(md.results.SteadystateSolution.BasalforcingsMeltingRate);
    7688
    7789        if md.inversion.iscontrol==1,
    7890                for control_parameters=md.inversion.control_parameters
    79                         md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters));
     91                        md.(EnumToModelField(control_parameters))=PatchToVec(md.results.SteadystateSolution.(EnumToString(control_parameters)));
    8092                end
    8193        end
    8294
    8395elseif strcmpi(string,'thermal'),
    84         md.initialization.temperature=md.results.ThermalSolution.Temperature;
    85         md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate;
     96        md.initialization.temperature=PatchToVec(md.results.ThermalSolution.Temperature);
     97        md.basalforcings.melting_rate=PatchToVec(md.results.ThermalSolution.BasalMeltingRate);
    8698elseif strcmpi(string,'hydrology'),
    87         md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn;
     99        md.initialization.watercolumn=PatchToVec(md.results.HydrologySolution.Watercolumn);
    88100
    89101else

  • issm/trunk-jpl-damage/src/m/qmu/expandresponses.m

    r11417 r11425  
    11function dresp=expandresponses(md,responses)
    2 %EXPANDRESPONSES - expand responses
    32
    43fnames=fieldnames(responses);

  • issm/trunk-jpl-damage/src/m/solutions/AnalysisConfiguration.m

    r11417 r11425  
    4242
    4343        case TransientSolutionEnum,
    44                 numanalyses=9;
    45                 analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum];
     44                numanalyses=8;
     45                analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum];
    4646
    4747        case FlaimSolutionEnum,

  • issm/trunk-jpl-damage/src/m/solutions/diagnostic_core.m

    r11417 r11425  
    1515        ismacayealpattyn=femmodel.parameters.FlowequationIsmacayealpattyn;
    1616        isstokes=femmodel.parameters.FlowequationIsstokes;
    17         isnewton=femmodel.parameters.DiagnosticIsnewton;
    1817        dakota_analysis=femmodel.parameters.QmuIsdakota;
    1918        control_analysis=femmodel.parameters.InversionIscontrol;
    … …  
    5453                issmprintf(VerboseSolution,'\n%s',['   computing horizontal velocities']);
    5554                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    56                 if isnewton,
    57                         femmodel=solver_newton(femmodel);
    58                 else
    59                         femmodel=solver_nonlinear(femmodel,modify_loads);
    60                 end
     55                femmodel=solver_nonlinear(femmodel,modify_loads);
    6156        end
    6257       

  • issm/trunk-jpl-damage/src/m/solutions/transient_core.m

    r11417 r11425  
    1919        isthermal=femmodel.parameters.TransientIsthermal;
    2020        isgroundingline=femmodel.parameters.TransientIsgroundingline;
    21         isenthalpy=femmodel.parameters.ThermalIsenthalpy;
    2221        groundinglinemigration=femmodel.parameters.GroundinglineMigration;
    2322
    … …  
    5857                if (isthermal & dim==3)
    5958                        issmprintf(VerboseSolution,'\n%s',['   computing temperature']);
    60                         if (isenthalpy==0),
    61                                 femmodel=thermal_core_step(femmodel);
    62                         else
    63                                 femmodel=enthalpy_core_step(femmodel);
    64                         end
     59                        femmodel=thermal_core_step(femmodel);
    6560                end
    6661
    … …  
    9590                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time);
    9691                        if (dim==3 & isthermal), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,step,time);end
    97                         if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,WaterfractionEnum,step,time);end
    98                         if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,EnthalpyEnum,step,time);end
    9992                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time);
    10093                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceforcingsMassBalanceEnum,step,time);

  • issm/trunk-jpl-damage/src/mex/Makefile.am

    r11417 r11425  
    2222                                CostFunction \
    2323                                CreateNodalConstraints\
    24                                 CreateJacobianMatrix\
    2524                                Echo\
    2625                                ElementConnectivity\
    … …  
    328327                          SystemMatrices/SystemMatrices.h
    329328
    330 
    331 CreateJacobianMatrix_SOURCES = CreateJacobianMatrix/CreateJacobianMatrix.cpp\
    332                                                                  CreateJacobianMatrix/CreateJacobianMatrix.h
    333 
    334329SurfaceArea_SOURCES = SurfaceArea/SurfaceArea.cpp\
    335330                                                                 SurfaceArea/SurfaceArea.h

  • issm/trunk-jpl-damage/test/NightlyRun/IdToName.m

    r11412 r11425  
    4848        case 141, name='SquareShelfConstrainedEnthalpyTranSerial';
    4949        case 142, name='SquareShelfConstrainedEnthalpyTranParallel';
    50         case 143, name='SquareShelfConstrainedTransP3dEnthSerial';
    51         case 144, name='SquareShelfConstrainedTransP3dEnthParallel';
    5250        case 201, name='SquareShelfDiagM2dSerial';
    5351        case 202, name='SquareShelfDiagM2dParallel';
    … …  
    170168        case 351, name='SquareSheetConstrainedEnthalpyTranSerial';
    171169        case 352, name='SquareSheetConstrainedEnthalpyTranParallel';
    172         case 353, name='SquareSheetConstrainedTransP3dEnthSerial';
    173         case 354, name='SquareSheetConstrainedTransP3dEnthParallel';
    174170        case 401, name='SquareSheetShelfDiagM2dSerial';
    175171        case 402, name='SquareSheetShelfDiagM2dParallel';

  • issm/trunk-jpl-damage/test/NightlyRun/test314.m

    r11412 r11425  
    99%Fields and tolerances to track changes
    1010field_names     ={'Vx','Vy','Vz','Vel','Pressure'};
    11 field_tolerances={1e-09,1e-09,1e-10,2e-10,1e-10};
     11field_tolerances={1e-09,1e-09,1e-10,1e-10,1e-10};
    1212field_values={...
    1313        (md.results.DiagnosticSolution.Vx),...

  • issm/trunk-jpl-damage/test/NightlyRun/test350.m

    r11412 r11425  
    77md.cluster=generic('name',oshostname(),'np',3);
    88md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
     9%md.debug.valgrind=1;
    910md=solve(md,EnthalpySolutionEnum);
     11%md=solve(md,ThermalSolutionEnum);
    1012
    1113%Fields and tolerances to track changes
Note: See TracChangeset for help on using the changeset viewer.