Changeset 11417


Ignore:
Timestamp:
02/13/12 15:50:19 (13 years ago)
Author:
cborstad
Message:

merged src changes 11330:11410 from trunk-jpl

Location:
issm/trunk-jpl-damage/src
Files:
66 edited
6 copied

Legend:

Unmodified
Added
Removed

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

    r11368 r11417  
    162162        ThermalSpctemperatureEnum,
    163163        ThermalStabilizationEnum,
     164        ThermalIsenthalpyEnum,
    164165        ThicknessEnum,
    165166        TimesteppingCflCoefficientEnum,

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

    r11225 r11417  
    99rm $ISSM_TIER/src/c/modules/EnumToStringx/EnumToStringx.cpp
    1010rm $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp
     11
     12#Get number of enums
     13NUMENUMS=$(wc -l temp | awk '{printf("%s",$1);}');
    1114
    1215#Take care of EnumToModelField.m first (easy)
    … …  
    9598int  StringToEnumx(const char* name){
    9699
     100   int  stage=1;
     101
    97102END
     103
    98104#core
    99 cat 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
     105i1=1;
     106i2=120;
     107for (( 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
     117done
     118
     119#footer
    101120cat <<END >> $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp
    102         else _error_("Enum %s not found",name);
    103 
     121        /*If we reach this point, the string provided has not been found*/
     122   _error_("Enum %s not found",name);
    104123}
    105124END
    106125#}}}
    107126
    108 #get number of lines in temp
    109 NUMBEROFLINES=$(wc -l temp | awk '{printf("%s",$1);}');
    110 
    111127# go through the lines of temp
    112 for (( i=1 ; i<=$NUMBEROFLINES ; i++ )); do
     128for (( i=1 ; i<=$NUMENUMS ; i++ )); do
    113129
    114130        #Get name and enum of the line i
    … …  
    123139        then
    124140                printf "\r                                                                      "
    125                 printf "\r  $i/$NUMBEROFLINES Adding "$NAME"..."
     141                printf "\r  $i/$NUMENUMS Adding "$NAME"..."
    126142        else
    127143                if [ $i -lt 100 ]
    128144                then
    129145                        printf "\r                                                                      "
    130                         printf "\r $i/$NUMBEROFLINES Adding "$NAME"..."
     146                        printf "\r $i/$NUMENUMS Adding "$NAME"..."
    131147                else
    132148                        printf "\r                                                                      "
    133                         printf "\r$i/$NUMBEROFLINES Adding "$NAME"..."
     149                        printf "\r$i/$NUMENUMS Adding "$NAME"..."
    134150                fi
    135151        fi
    … …  
    153169done
    154170
    155 
    156171#clean up{{{
    157172rm temp

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

    r10205 r11417  
    1717        Th.WriteIndex(pindex,pnels);
    1818        //delete &Th;
     19        return 0;
    1920
    2021}

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

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

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

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

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

    r9320 r11417  
    6969        }
    7070
     71        return NULL;
     72
    7173}

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

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

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

    r11368 r11417  
    166166                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
    167167                case ThermalStabilizationEnum : return "ThermalStabilization";
     168                case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
    168169                case ThicknessEnum : return "Thickness";
    169170                case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient";

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

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

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

    r11275 r11417  
    7878        parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
    7979        parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
     80        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
    8081        parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    8182        parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));

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

    r11129 r11417  
    2020
    2121        int   i,analysis_type,dim,verbose;
    22         bool  isthermal,isprognostic,isdiagnostic,isgroundingline;
     22        bool  isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy;
    2323       
    2424        /*output: */
    … …  
    3939        iomodel->Constant(&verbose,VerboseEnum);
    4040        iomodel->Constant(&isthermal,TransientIsthermalEnum);
     41        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
    4142        iomodel->Constant(&isprognostic,TransientIsprognosticEnum);
    4243        iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum);
    … …  
    5253                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && dim==2) continue;
    5354                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && dim==2) continue;
     55                if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && dim==2) continue;
    5456                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isthermal==false) continue;
    5557                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;
    5662                if(solution_type==TransientSolutionEnum && analysis_type==PrognosticAnalysisEnum && isprognostic==false && isgroundingline==false) continue;
    5763                if(solution_type==TransientSolutionEnum && analysis_type==DiagnosticHorizAnalysisEnum && isdiagnostic==false) continue;

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

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

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

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

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

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

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

    r10380 r11417  
    2424                                        sgn,cm,sp));
    2525
    26         #endif //ifdef _HAVE_SHAPELIB_
     26        #else //ifdef _HAVE_SHAPELIB_
     27        return 0;
     28        #endif
    2729}
    2830
    … …  
    607609        return(iret);
    608610
    609         #endif //ifdef _HAVE_SHAPELIB_
     611        #else //ifdef _HAVE_SHAPELIB_
     612        return 0;
     613        #endif
    610614}
    611615

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

    r11368 r11417  
    1414int  StringToEnumx(const char* name){
    1515
     16<<<<<<< .working
    1617        if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum;
    1718        else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum;
    … …  
    443444        else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    444445        else _error_("Enum %s not found",name);
     446=======
     447   int  stage=1;
     448>>>>>>> .merge-right.r11410
    445449
     450   if(stage==1){
     451              if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum;
     452              else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum;
     453              else if (strcmp(name,"AutodiffReverse")==0) return AutodiffReverseEnum;
     454              else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
     455              else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum;
     456              else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum;
     457              else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
     458              else if (strcmp(name,"BasalforcingsMeltingRateCorrection")==0) return BasalforcingsMeltingRateCorrectionEnum;
     459              else if (strcmp(name,"BasalforcingsMeltingRate")==0) return BasalforcingsMeltingRateEnum;
     460              else if (strcmp(name,"Bathymetry")==0) return BathymetryEnum;
     461              else if (strcmp(name,"Bed")==0) return BedEnum;
     462              else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
     463              else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
     464              else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
     465              else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum;
     466              else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum;
     467              else if (strcmp(name,"DiagnosticIsnewton")==0) return DiagnosticIsnewtonEnum;
     468              else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum;
     469              else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum;
     470              else if (strcmp(name,"DiagnosticPenaltyFactor")==0) return DiagnosticPenaltyFactorEnum;
     471              else if (strcmp(name,"DiagnosticReferential")==0) return DiagnosticReferentialEnum;
     472              else if (strcmp(name,"DiagnosticReltol")==0) return DiagnosticReltolEnum;
     473              else if (strcmp(name,"DiagnosticRequestedOutputs")==0) return DiagnosticRequestedOutputsEnum;
     474              else if (strcmp(name,"DiagnosticRestol")==0) return DiagnosticRestolEnum;
     475              else if (strcmp(name,"DiagnosticRiftPenaltyLock")==0) return DiagnosticRiftPenaltyLockEnum;
     476              else if (strcmp(name,"DiagnosticRiftPenaltyThreshold")==0) return DiagnosticRiftPenaltyThresholdEnum;
     477              else if (strcmp(name,"DiagnosticShelfDampening")==0) return DiagnosticShelfDampeningEnum;
     478              else if (strcmp(name,"DiagnosticSpcvx")==0) return DiagnosticSpcvxEnum;
     479              else if (strcmp(name,"DiagnosticSpcvy")==0) return DiagnosticSpcvyEnum;
     480              else if (strcmp(name,"DiagnosticSpcvz")==0) return DiagnosticSpcvzEnum;
     481              else if (strcmp(name,"DiagnosticStokesreconditioning")==0) return DiagnosticStokesreconditioningEnum;
     482              else if (strcmp(name,"DiagnosticVertexPairing")==0) return DiagnosticVertexPairingEnum;
     483              else if (strcmp(name,"DiagnosticViscosityOvershoot")==0) return DiagnosticViscosityOvershootEnum;
     484              else if (strcmp(name,"FlowequationBordermacayeal")==0) return FlowequationBordermacayealEnum;
     485              else if (strcmp(name,"FlowequationBorderpattyn")==0) return FlowequationBorderpattynEnum;
     486              else if (strcmp(name,"FlowequationBorderstokes")==0) return FlowequationBorderstokesEnum;
     487              else if (strcmp(name,"FlowequationElementEquation")==0) return FlowequationElementEquationEnum;
     488              else if (strcmp(name,"FlowequationIshutter")==0) return FlowequationIshutterEnum;
     489              else if (strcmp(name,"FlowequationIsmacayealpattyn")==0) return FlowequationIsmacayealpattynEnum;
     490              else if (strcmp(name,"FlowequationIsstokes")==0) return FlowequationIsstokesEnum;
     491              else if (strcmp(name,"FlowequationVertexEquation")==0) return FlowequationVertexEquationEnum;
     492              else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
     493              else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
     494              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
     495              else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
     496              else if (strcmp(name,"HydrologyCR")==0) return HydrologyCREnum;
     497              else if (strcmp(name,"HydrologyKn")==0) return HydrologyKnEnum;
     498              else if (strcmp(name,"HydrologyN")==0) return HydrologyNEnum;
     499              else if (strcmp(name,"HydrologyP")==0) return HydrologyPEnum;
     500              else if (strcmp(name,"HydrologyQ")==0) return HydrologyQEnum;
     501              else if (strcmp(name,"HydrologySpcwatercolumn")==0) return HydrologySpcwatercolumnEnum;
     502              else if (strcmp(name,"HydrologyStabilization")==0) return HydrologyStabilizationEnum;
     503              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
     504              else if (strcmp(name,"InversionCostFunction")==0) return InversionCostFunctionEnum;
     505              else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
     506              else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
     507              else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
     508              else if (strcmp(name,"InversionGradientOnly")==0) return InversionGradientOnlyEnum;
     509              else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
     510              else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
     511              else if (strcmp(name,"InversionTao")==0) return InversionTaoEnum;
     512              else if (strcmp(name,"InversionMaxParameters")==0) return InversionMaxParametersEnum;
     513              else if (strcmp(name,"InversionMaxiterPerStep")==0) return InversionMaxiterPerStepEnum;
     514              else if (strcmp(name,"InversionMinParameters")==0) return InversionMinParametersEnum;
     515              else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
     516              else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
     517              else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
     518              else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
     519              else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
     520              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
     521              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
     522              else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
     523              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
     524              else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
     525              else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
     526              else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum;
     527              else if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum;
     528              else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
     529              else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
     530              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
     531              else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
     532              else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
     533              else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
     534              else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
     535              else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
     536              else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
     537              else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
     538              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
     539              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
     540              else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
     541              else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
     542              else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
     543              else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
     544              else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
     545              else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum;
     546              else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum;
     547              else if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum;
     548              else if (strcmp(name,"MeshElementonbed")==0) return MeshElementonbedEnum;
     549              else if (strcmp(name,"MeshElementonsurface")==0) return MeshElementonsurfaceEnum;
     550              else if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum;
     551              else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
     552              else if (strcmp(name,"MeshLowerelements")==0) return MeshLowerelementsEnum;
     553              else if (strcmp(name,"MeshNumberofedges")==0) return MeshNumberofedgesEnum;
     554              else if (strcmp(name,"MeshNumberofelements2d")==0) return MeshNumberofelements2dEnum;
     555              else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum;
     556              else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
     557              else if (strcmp(name,"MeshNumberofvertices2d")==0) return MeshNumberofvertices2dEnum;
     558              else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
     559              else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum;
     560              else if (strcmp(name,"MeshVertexonbed")==0) return MeshVertexonbedEnum;
     561              else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
     562              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
     563              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
     564              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
     565              else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
     566              else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum;
     567              else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
     568              else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
     569              else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
     570              else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
     571         else stage=2;
     572   }
     573   if(stage==2){
     574              if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
     575              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
     576              else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum;
     577              else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
     578              else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
     579              else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
     580              else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
     581              else if (strcmp(name,"QmuVariabledescriptors")==0) return QmuVariabledescriptorsEnum;
     582              else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum;
     583              else if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum;
     584              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
     585              else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
     586              else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
     587              else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum;
     588              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
     589              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
     590              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
     591              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
     592              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
     593              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
     594              else if (strcmp(name,"SurfaceforcingsAblationRate")==0) return SurfaceforcingsAblationRateEnum;
     595              else if (strcmp(name,"SurfaceforcingsAccumulationRate")==0) return SurfaceforcingsAccumulationRateEnum;
     596              else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
     597              else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
     598              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
     599              else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     600              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
     601              else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
     602              else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
     603              else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
     604              else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
     605              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
     606              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
     607              else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
     608              else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
     609              else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
     610              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
     611              else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
     612              else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
     613              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
     614              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
     615              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
     616              else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
     617              else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
     618              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     619              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
     620              else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
     621              else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
     622              else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
     623              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
     624              else if (strcmp(name,"BedSlopeAnalysis")==0) return BedSlopeAnalysisEnum;
     625              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     626              else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum;
     627              else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum;
     628              else if (strcmp(name,"DiagnosticHorizAnalysis")==0) return DiagnosticHorizAnalysisEnum;
     629              else if (strcmp(name,"DiagnosticHutterAnalysis")==0) return DiagnosticHutterAnalysisEnum;
     630              else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum;
     631              else if (strcmp(name,"DiagnosticVertAnalysis")==0) return DiagnosticVertAnalysisEnum;
     632              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
     633              else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum;
     634              else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
     635              else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
     636              else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum;
     637              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
     638              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
     639              else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum;
     640              else if (strcmp(name,"PrognosticAnalysis")==0) return PrognosticAnalysisEnum;
     641              else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum;
     642              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
     643              else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
     644              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
     645              else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
     646              else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
     647              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
     648              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
     649              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
     650              else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
     651              else if (strcmp(name,"HutterApproximation")==0) return HutterApproximationEnum;
     652              else if (strcmp(name,"MacAyealApproximation")==0) return MacAyealApproximationEnum;
     653              else if (strcmp(name,"MacAyealPattynApproximation")==0) return MacAyealPattynApproximationEnum;
     654              else if (strcmp(name,"MacAyealStokesApproximation")==0) return MacAyealStokesApproximationEnum;
     655              else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
     656              else if (strcmp(name,"PattynApproximation")==0) return PattynApproximationEnum;
     657              else if (strcmp(name,"PattynStokesApproximation")==0) return PattynStokesApproximationEnum;
     658              else if (strcmp(name,"StokesApproximation")==0) return StokesApproximationEnum;
     659              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     660              else if (strcmp(name,"Loads")==0) return LoadsEnum;
     661              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
     662              else if (strcmp(name,"Nodes")==0) return NodesEnum;
     663              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
     664              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
     665              else if (strcmp(name,"Results")==0) return ResultsEnum;
     666              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     667              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     668              else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
     669              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
     670              else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
     671              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
     672              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
     673              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
     674              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
     675              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
     676              else if (strcmp(name,"Element")==0) return ElementEnum;
     677              else if (strcmp(name,"ElementResult")==0) return ElementResultEnum;
     678              else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
     679              else if (strcmp(name,"FileParam")==0) return FileParamEnum;
     680              else if (strcmp(name,"Hook")==0) return HookEnum;
     681              else if (strcmp(name,"Icefront")==0) return IcefrontEnum;
     682              else if (strcmp(name,"Input")==0) return InputEnum;
     683              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
     684              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
     685              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
     686              else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
     687              else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
     688              else if (strcmp(name,"Matice")==0) return MaticeEnum;
     689              else if (strcmp(name,"Matpar")==0) return MatparEnum;
     690              else if (strcmp(name,"Node")==0) return NodeEnum;
     691              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
     692              else if (strcmp(name,"Param")==0) return ParamEnum;
     693              else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
     694         else stage=3;
     695   }
     696   if(stage==3){
     697              if (strcmp(name,"Pengrid")==0) return PengridEnum;
     698              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
     699              else if (strcmp(name,"Penta")==0) return PentaEnum;
     700              else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum;
     701              else if (strcmp(name,"PetscMatParam")==0) return PetscMatParamEnum;
     702              else if (strcmp(name,"PetscVecParam")==0) return PetscVecParamEnum;
     703              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     704              else if (strcmp(name,"Segment")==0) return SegmentEnum;
     705              else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     706              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
     707              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
     708              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
     709              else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
     710              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
     711              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
     712              else if (strcmp(name,"Tria")==0) return TriaEnum;
     713              else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
     714              else if (strcmp(name,"Vertex")==0) return VertexEnum;
     715              else if (strcmp(name,"Air")==0) return AirEnum;
     716              else if (strcmp(name,"Ice")==0) return IceEnum;
     717              else if (strcmp(name,"Melange")==0) return MelangeEnum;
     718              else if (strcmp(name,"Water")==0) return WaterEnum;
     719              else if (strcmp(name,"Closed")==0) return ClosedEnum;
     720              else if (strcmp(name,"Free")==0) return FreeEnum;
     721              else if (strcmp(name,"Open")==0) return OpenEnum;
     722              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
     723              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
     724              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
     725              else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
     726              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
     727              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
     728              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     729              else if (strcmp(name,"Constant")==0) return ConstantEnum;
     730              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
     731              else if (strcmp(name,"ExtToIu")==0) return ExtToIuEnum;
     732              else if (strcmp(name,"Fill")==0) return FillEnum;
     733              else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
     734              else if (strcmp(name,"Friction")==0) return FrictionEnum;
     735              else if (strcmp(name,"GroundinglineMeltingRate")==0) return GroundinglineMeltingRateEnum;
     736              else if (strcmp(name,"Internal")==0) return InternalEnum;
     737              else if (strcmp(name,"IuToExt")==0) return IuToExtEnum;
     738              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
     739              else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum;
     740              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
     741              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
     742              else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum;
     743              else if (strcmp(name,"Pressure")==0) return PressureEnum;
     744              else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
     745              else if (strcmp(name,"QmuPressure")==0) return QmuPressureEnum;
     746              else if (strcmp(name,"QmuVx")==0) return QmuVxEnum;
     747              else if (strcmp(name,"QmuVy")==0) return QmuVyEnum;
     748              else if (strcmp(name,"QmuVz")==0) return QmuVzEnum;
     749              else if (strcmp(name,"QmuThickness")==0) return QmuThicknessEnum;
     750              else if (strcmp(name,"QmuBed")==0) return QmuBedEnum;
     751              else if (strcmp(name,"QmuSurface")==0) return QmuSurfaceEnum;
     752              else if (strcmp(name,"QmuMelting")==0) return QmuMeltingEnum;
     753              else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
     754              else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
     755              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
     756              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
     757              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
     758              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
     759              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     760              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
     761              else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
     762              else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
     763              else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
     764              else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
     765              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
     766              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
     767              else if (strcmp(name,"Type")==0) return TypeEnum;
     768              else if (strcmp(name,"Vel")==0) return VelEnum;
     769              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
     770              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     771              else if (strcmp(name,"Vx")==0) return VxEnum;
     772              else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
     773              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     774              else if (strcmp(name,"Vy")==0) return VyEnum;
     775              else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
     776              else if (strcmp(name,"Vz")==0) return VzEnum;
     777              else if (strcmp(name,"VzMacAyeal")==0) return VzMacAyealEnum;
     778              else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum;
     779              else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
     780              else if (strcmp(name,"VzStokes")==0) return VzStokesEnum;
     781              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
     782              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
     783              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
     784              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
     785              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
     786              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
     787              else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
     788              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
     789              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
     790              else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
     791              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
     792              else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum;
     793              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
     794              else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
     795              else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
     796              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
     797              else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
     798              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
     799              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
     800              else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
     801              else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
     802              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
     803              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
     804              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     805              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
     806              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
     807              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     808              else if (strcmp(name,"P0")==0) return P0Enum;
     809              else if (strcmp(name,"P1")==0) return P1Enum;
     810              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
     811              else if (strcmp(name,"BoolElementResult")==0) return BoolElementResultEnum;
     812              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     813              else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
     814              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     815              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     816              else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum;
     817         else stage=4;
     818   }
     819   if(stage==4){
     820              if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     821              else if (strcmp(name,"J")==0) return JEnum;
     822              else if (strcmp(name,"Patch")==0) return PatchEnum;
     823              else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum;
     824              else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum;
     825              else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
     826              else if (strcmp(name,"PetscVecExternalResult")==0) return PetscVecExternalResultEnum;
     827              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
     828              else if (strcmp(name,"Time")==0) return TimeEnum;
     829              else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum;
     830              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
     831              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
     832              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
     833              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
     834              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
     835              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
     836              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
     837              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
     838              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
     839              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
     840              else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
     841              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
     842              else if (strcmp(name,"Relative")==0) return RelativeEnum;
     843              else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
     844              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     845              else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
     846              else if (strcmp(name,"None")==0) return NoneEnum;
     847              else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
     848              else if (strcmp(name,"StokesSolver")==0) return StokesSolverEnum;
     849              else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
     850              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
     851              else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
     852              else if (strcmp(name,"Fset")==0) return FsetEnum;
     853              else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
     854              else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
     855              else if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
     856              else if (strcmp(name,"Gradient")==0) return GradientEnum;
     857              else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
     858              else if (strcmp(name,"Gset")==0) return GsetEnum;
     859              else if (strcmp(name,"Index")==0) return IndexEnum;
     860              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
     861              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
     862              else if (strcmp(name,"Nodal")==0) return NodalEnum;
     863              else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
     864              else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
     865              else if (strcmp(name,"PetscOptionsAnalyses")==0) return PetscOptionsAnalysesEnum;
     866              else if (strcmp(name,"PetscOptionsStrings")==0) return PetscOptionsStringsEnum;
     867              else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
     868              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
     869              else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
     870              else if (strcmp(name,"Regular")==0) return RegularEnum;
     871              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
     872              else if (strcmp(name,"Separate")==0) return SeparateEnum;
     873              else if (strcmp(name,"Sset")==0) return SsetEnum;
     874              else if (strcmp(name,"Verbose")==0) return VerboseEnum;
     875              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
     876              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     877              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
     878              else if (strcmp(name,"XY")==0) return XYEnum;
     879              else if (strcmp(name,"XYZP")==0) return XYZPEnum;
     880              else if (strcmp(name,"Option")==0) return OptionEnum;
     881              else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
     882              else if (strcmp(name,"OptionChar")==0) return OptionCharEnum;
     883              else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
     884              else if (strcmp(name,"OptionDouble")==0) return OptionDoubleEnum;
     885              else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum;
     886              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
     887              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
     888         else stage=5;
     889   }
     890        /*If we reach this point, the string provided has not been found*/
     891   _error_("Enum %s not found",name);
    446892}

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

    r11409 r11417  
    769769}
    770770/*}}}*/
     771<<<<<<< .working
     772=======
     773/*FUNCTION Penta::CreateJacobianMatrix{{{1*/
     774void  Penta::CreateJacobianMatrix(Mat Jff){
     775
     776        /*retrieve parameters: */
     777        ElementMatrix* Ke=NULL;
     778        int analysis_type;
     779        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     780
     781        /*Checks in debugging {{{2*/
     782        _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
     783        /*}}}*/
     784
     785        /*Skip if water element*/
     786        if(IsOnWater()) return;
     787
     788        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     789        switch(analysis_type){
     790#ifdef _HAVE_DIAGNOSTIC_
     791                case DiagnosticHorizAnalysisEnum:
     792                        Ke=CreateJacobianDiagnosticHoriz();
     793                        break;
     794#endif
     795                default:
     796                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
     797        }
     798
     799        /*Add to global matrix*/
     800        if(Ke){
     801                Ke->AddToGlobal(Jff);
     802                delete Ke;
     803        }
     804}
     805/*}}}*/
     806>>>>>>> .merge-right.r11410
    771807/*FUNCTION Penta::DeepEcho{{{1*/
    772808void Penta::DeepEcho(void){
    … …  
    10991135/*}}}*/
    11001136/*FUNCTION Penta::GetStabilizationParameter {{{1*/
    1101 double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
     1137double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){
    11021138        /*Compute stabilization parameter*/
     1139        /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
     1140        /*kappa=enthalpydiffusionparameter for enthalpy model*/
    11031141
    11041142        double normu;
    … …  
    11061144
    11071145        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
    1108         if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
    1109                 tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
     1146        if(normu*diameter/(3*2*kappa)<1){
     1147                tau_parameter=pow(diameter,2)/(3*2*2*kappa);
    11101148        }
    11111149        else tau_parameter=diameter/(2*normu);
    … …  
    32833321                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    32843322
    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;
     3323                vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     3324                vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     3325                vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    32923326
    32933327                D_scalar_advec=gauss->weight*Jdet;
    … …  
    33213355                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    33223356                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
    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);
     3357                        K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
     3358                        K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
     3359                        K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
    33263360                        D_scalar_stab=gauss->weight*Jdet;
    33273361                        if(dt) D_scalar_stab=D_scalar_stab*dt;
    … …  
    33373371                else if(stabilization==2){
    33383372                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3339                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
     3373                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    33403374
    33413375                        for(i=0;i<numdof;i++){
    … …  
    34513485        double     Jdet,u,v,w,um,vm,wm,vel;
    34523486        double     h,hx,hy,hz,vx,vy,vz;
    3453         double     gravity,rho_ice,rho_water;
     3487        double     gravity,rho_ice,rho_water,kappa;
    34543488        double     heatcapacity,thermalconductivity,dt;
    34553489        double     tau_parameter,diameter;
    … …  
    34773511        heatcapacity=matpar->GetHeatCapacity();
    34783512        thermalconductivity=matpar->GetThermalConductivity();
     3513        kappa=thermalconductivity/(rho_ice*heatcapacity);
    34793514        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    34803515        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    … …  
    34993534                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    35003535
    3501                 D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
     3536                D_scalar_conduct=gauss->weight*Jdet*kappa;
    35023537                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
    35033538
    … …  
    35683603                else if(stabilization==2){
    35693604                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3570                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
     3605                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    35713606
    35723607                        for(i=0;i<numdof;i++){
    … …  
    36733708        double Jdet,phi,dt;
    36743709        double rho_ice,heatcapacity;
    3675         double thermalconductivity;
    3676         double viscosity,enthalpy;
     3710        double thermalconductivity,kappa;
     3711        double viscosity,pressure;
     3712        double enthalpy,enthalpypicard;
    36773713        double tau_parameter,diameter;
    36783714        double u,v,w;
    … …  
    36953731        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    36963732        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    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);
     3733        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     3734        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     3735        Input* vz_input=inputs->GetInput(VzEnum);                                  _assert_(vz_input);
     3736        Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
     3737        Input* enthalpy_input=NULL;
     3738        Input* enthalpypicard_input=NULL;
     3739        if(dt){
     3740                enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     3741        }
     3742        if (stabilization==2){
     3743                diameter=MinEdgeLength(xyz_list);
     3744                enthalpypicard_input=inputs->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
     3745        }
    37033746
    37043747        /* Start  looping on the number of gaussian points: */
    … …  
    37153758                GetPhi(&phi, &epsilon[0], viscosity);
    37163759
    3717                 scalar_def=phi/(rho_ice)*Jdet*gauss->weight;
     3760                scalar_def=phi/rho_ice*Jdet*gauss->weight;
    37183761                if(dt) scalar_def=scalar_def*dt;
    37193762
    … …  
    37333776                        vy_input->GetInputValue(&v, gauss);
    37343777                        vz_input->GetInputValue(&w, gauss);
    3735 
    3736                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
     3778                        pressure_input->GetInputValue(&pressure, gauss);
     3779                        enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
     3780                        kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
     3781                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
    37373782
    37383783                        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]);
    … …  
    38193864        double     rho_ice,heatcapacity,geothermalflux_value;
    38203865        double     basalfriction,alpha2,vx,vy;
     3866        double     scalar,enthalpy,enthalpyup;
     3867        double     pressure,pressureup;
     3868        double     basis[NUMVERTICES];
     3869        Friction*  friction=NULL;
     3870        GaussPenta* gauss=NULL;
     3871        GaussPenta* gaussup=NULL;
     3872
     3873        /* Geothermal flux on ice sheet base and basal friction */
     3874        if (!IsOnBed() || IsFloating()) return NULL;
     3875
     3876        /*Initialize Element vector*/
     3877        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3878
     3879        /*Retrieve all inputs and parameters*/
     3880        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3881        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3882        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     3883        rho_ice=matpar->GetRhoIce();
     3884        heatcapacity=matpar->GetHeatCapacity();
     3885        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     3886        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
     3887        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
     3888        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     3889        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
     3890        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
     3891        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
     3892
     3893        /*Build frictoin element, needed later: */
     3894        friction=new Friction("3d",inputs,matpar,analysis_type);
     3895
     3896        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3897        gauss=new GaussPenta(0,1,2,2);
     3898        gaussup=new GaussPenta(3,4,5,2);
     3899        for(ig=gauss->begin();ig<gauss->end();ig++){
     3900
     3901                gauss->GaussPoint(ig);
     3902                gaussup->GaussPoint(ig);
     3903
     3904                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3905                GetNodalFunctionsP1(&basis[0], gauss);
     3906
     3907                enthalpy_input->GetInputValue(&enthalpy,gauss);
     3908                pressure_input->GetInputValue(&pressure,gauss);
     3909//              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
     3910//                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     3911//                      pressure_input->GetInputValue(&pressureup,gaussup);
     3912//                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
     3913//                              //do nothing, don't add heatflux
     3914//                      }
     3915//                      else{
     3916//                              //need to change spcenthalpy according to Aschwanden
     3917//                              //NEED TO UPDATE
     3918//                      }
     3919//              }
     3920//              else{
     3921                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     3922                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3923                        vx_input->GetInputValue(&vx,gauss);
     3924                        vy_input->GetInputValue(&vy,gauss);
     3925                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     3926
     3927                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     3928                        if(dt) scalar=dt*scalar;
     3929
     3930                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     3931//              }
     3932        }
     3933
     3934        /*Clean up and return*/
     3935        delete gauss;
     3936        delete gaussup;
     3937        delete friction;
     3938        return pe;
     3939}
     3940/*}}}*/
     3941/*FUNCTION Penta::CreatePVectorMelting {{{1*/
     3942ElementVector* Penta::CreatePVectorMelting(void){
     3943        return NULL;
     3944}
     3945/*}}}*/
     3946/*FUNCTION Penta::CreatePVectorThermal {{{1*/
     3947ElementVector* Penta::CreatePVectorThermal(void){
     3948
     3949        /*compute all load vectors for this element*/
     3950        ElementVector* pe1=CreatePVectorThermalVolume();
     3951        ElementVector* pe2=CreatePVectorThermalSheet();
     3952        ElementVector* pe3=CreatePVectorThermalShelf();
     3953        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3954
     3955        /*clean-up and return*/
     3956        delete pe1;
     3957        delete pe2;
     3958        delete pe3;
     3959        return pe;
     3960}
     3961/*}}}*/
     3962/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
     3963ElementVector* Penta::CreatePVectorThermalVolume(void){
     3964
     3965        /*Constants*/
     3966        const int  numdof=NUMVERTICES*NDOF1;
     3967
     3968        /*Intermediaries*/
     3969        int    i,j,ig,found=0;
     3970        int    friction_type,stabilization;
     3971        double Jdet,phi,dt;
     3972        double rho_ice,heatcapacity;
     3973        double thermalconductivity,kappa;
     3974        double viscosity,temperature;
     3975        double tau_parameter,diameter;
     3976        double u,v,w;
     3977        double scalar_def,scalar_transient;
     3978        double temperature_list[NUMVERTICES];
     3979        double xyz_list[NUMVERTICES][3];
     3980        double L[numdof];
     3981        double dbasis[3][6];
     3982        double epsilon[6];
     3983        GaussPenta *gauss=NULL;
     3984
     3985        /*Initialize Element vector*/
     3986        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3987
     3988        /*Retrieve all inputs and parameters*/
     3989        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3990        rho_ice=matpar->GetRhoIce();
     3991        heatcapacity=matpar->GetHeatCapacity();
     3992        thermalconductivity=matpar->GetThermalConductivity();
     3993        kappa=thermalconductivity/(rho_ice*heatcapacity);
     3994        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     3995        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
     3996        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3997        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3998        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     3999        Input* temperature_input=NULL;
     4000        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
     4001        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
     4002
     4003        /* Start  looping on the number of gaussian points: */
     4004        gauss=new GaussPenta(2,3);
     4005        for (ig=gauss->begin();ig<gauss->end();ig++){
     4006
     4007                gauss->GaussPoint(ig);
     4008
     4009                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     4010                GetNodalFunctionsP1(&L[0], gauss);
     4011
     4012                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     4013                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     4014                GetPhi(&phi, &epsilon[0], viscosity);
     4015
     4016                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
     4017                if(dt) scalar_def=scalar_def*dt;
     4018
     4019                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
     4020
     4021                /* Build transient now */
     4022                if(dt){
     4023                        temperature_input->GetInputValue(&temperature, gauss);
     4024                        scalar_transient=temperature*Jdet*gauss->weight;
     4025                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     4026                }
     4027
     4028                if(stabilization==2){
     4029                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
     4030
     4031                        vx_input->GetInputValue(&u, gauss);
     4032                        vy_input->GetInputValue(&v, gauss);
     4033                        vz_input->GetInputValue(&w, gauss);
     4034
     4035                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
     4036
     4037                        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]);
     4038                        if(dt){
     4039                                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]);
     4040                        }
     4041                }
     4042        }
     4043
     4044        /*Clean up and return*/
     4045        delete gauss;
     4046        return pe;
     4047}
     4048/*}}}*/
     4049/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
     4050ElementVector* Penta::CreatePVectorThermalShelf(void){
     4051
     4052        /*Constants*/
     4053        const int  numdof=NUMVERTICES*NDOF1;
     4054
     4055        /*Intermediaries */
     4056        int        i,j,ig;
     4057        double     Jdet2d;
     4058        double     mixed_layer_capacity,thermal_exchange_velocity;
     4059        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     4060        double     heatcapacity,t_pmp;
     4061        double     xyz_list[NUMVERTICES][3];
     4062        double     xyz_list_tria[NUMVERTICES2D][3];
     4063        double     basis[NUMVERTICES];
     4064        GaussPenta* gauss=NULL;
     4065
     4066        /* Ice/ocean heat exchange flux on ice shelf base */
     4067        if (!IsOnBed() || !IsFloating()) return NULL;
     4068
     4069        /*Initialize Element vector*/
     4070        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4071
     4072        /*Retrieve all inputs and parameters*/
     4073        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4074        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4075        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     4076        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     4077        rho_water=matpar->GetRhoWater();
     4078        rho_ice=matpar->GetRhoIce();
     4079        heatcapacity=matpar->GetHeatCapacity();
     4080        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4081        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
     4082
     4083        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     4084        gauss=new GaussPenta(0,1,2,2);
     4085        for(ig=gauss->begin();ig<gauss->end();ig++){
     4086
     4087                gauss->GaussPoint(ig);
     4088
     4089                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     4090                GetNodalFunctionsP1(&basis[0], gauss);
     4091
     4092                pressure_input->GetInputValue(&pressure,gauss);
     4093                t_pmp=matpar->TMeltingPoint(pressure);
     4094
     4095                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
     4096                if(dt) scalar_ocean=dt*scalar_ocean;
     4097
     4098                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
     4099        }
     4100
     4101        /*Clean up and return*/
     4102        delete gauss;
     4103        return pe;
     4104}
     4105/*}}}*/
     4106/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
     4107ElementVector* Penta::CreatePVectorThermalSheet(void){
     4108
     4109        /*Constants*/
     4110        const int  numdof=NUMVERTICES*NDOF1;
     4111
     4112        /*Intermediaries */
     4113        int        i,j,ig;
     4114        int        analysis_type;
     4115        double     xyz_list[NUMVERTICES][3];
     4116        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4117        double     Jdet2d,dt;
     4118        double     rho_ice,heatcapacity,geothermalflux_value;
     4119        double     basalfriction,alpha2,vx,vy;
     4120        double     basis[NUMVERTICES];
    38214121        double     scalar;
    3822         double     basis[NUMVERTICES];
    38234122        Friction*  friction=NULL;
    38244123        GaussPenta* gauss=NULL;
    … …  
    38544153                GetNodalFunctionsP1(&basis[0], gauss);
    38554154
    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*/
    3875 ElementVector* Penta::CreatePVectorMelting(void){
    3876         return NULL;
    3877 }
    3878 /*}}}*/
    3879 /*FUNCTION Penta::CreatePVectorThermal {{{1*/
    3880 ElementVector* 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*/
    3896 ElementVector* 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*/
    3982 ElementVector* 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*/
    4039 ElementVector* 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];
     4155                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4156                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4157                        vx_input->GetInputValue(&vx,gauss);
     4158                        vy_input->GetInputValue(&vy,gauss);
     4159                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4160
     4161                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     4162                        if(dt) scalar=dt*scalar;
     4163
     4164                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    40974165        }
    40984166
    … …  
    42714339        if(converged){
    42724340                /*Convert enthalpy into temperature and water fraction*/
     4341<<<<<<< .working
    42734342                for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
     4343=======
     4344                for(i=0;i<numdof;i++){
     4345                        matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
     4346                        if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
     4347                        //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
     4348                }
     4349>>>>>>> .merge-right.r11410
    42744350                       
    42754351                this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values));
    … …  
    73197395}
    73207396/*}}}*/
     7397/*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/
     7398ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){
     7399
     7400        int approximation;
     7401        inputs->GetInputValue(&approximation,ApproximationEnum);
     7402
     7403        switch(approximation){
     7404                case MacAyealApproximationEnum:
     7405                        return CreateJacobianDiagnosticMacayeal2d();
     7406                case PattynApproximationEnum:
     7407                        return CreateJacobianDiagnosticPattyn();
     7408                case StokesApproximationEnum:
     7409                        return CreateJacobianDiagnosticStokes();
     7410                case NoneApproximationEnum:
     7411                        return NULL;
     7412                default:
     7413                        _error_("Approximation %s not supported yet",EnumToStringx(approximation));
     7414        }
     7415}
     7416/*}}}*/
     7417/*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{1*/
     7418ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){
     7419
     7420        /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
     7421          bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
     7422          the stiffness matrix. */
     7423        if (!IsOnBed()) return NULL;
     7424
     7425        /*Depth Averaging B*/
     7426        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
     7427
     7428        /*Call Tria function*/
     7429        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     7430        ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal();
     7431        delete tria->matice; delete tria;
     7432
     7433        /*Delete B averaged*/
     7434        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     7435
     7436        /*clean up and return*/
     7437        return Ke;
     7438}
     7439/*}}}*/
     7440/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
     7441ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
     7442
     7443        /*Constants*/
     7444        const int    numdof=NDOF2*NUMVERTICES;
     7445
     7446        /*Intermediaries */
     7447        int        i,j,ig;
     7448        double     xyz_list[NUMVERTICES][3];
     7449        double     Jdet;
     7450        double     eps1dotdphii,eps1dotdphij;
     7451        double     eps2dotdphii,eps2dotdphij;
     7452        double     mu_prime;
     7453        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     7454        double     eps1[3],eps2[3];
     7455        double     phi[NUMVERTICES];
     7456        double     dphi[3][NUMVERTICES];
     7457        GaussPenta *gauss=NULL;
     7458
     7459        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
     7460        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
     7461
     7462        /*Retrieve all inputs and parameters*/
     7463        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     7464        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     7465        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     7466
     7467        /* Start  looping on the number of gaussian points: */
     7468        gauss=new GaussPenta(5,5);
     7469        for (ig=gauss->begin();ig<gauss->end();ig++){
     7470
     7471                gauss->GaussPoint(ig);
     7472
     7473                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7474                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     7475
     7476                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7477                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     7478                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     7479                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     7480                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
     7481
     7482                for(i=0;i<6;i++){
     7483                        for(j=0;j<6;j++){
     7484                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
     7485                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
     7486                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
     7487                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
     7488
     7489                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     7490                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     7491                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     7492                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     7493                        }
     7494                }
     7495        }
     7496
     7497        /*Transform Coordinate System*/
     7498        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
     7499
     7500        /*Clean up and return*/
     7501        delete gauss;
     7502        return Ke;
     7503}
     7504/*}}}*/
     7505/*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{1*/
     7506ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
     7507
     7508        /*Constants*/
     7509        const int    numdof=NDOF4*NUMVERTICES;
     7510
     7511        /*Intermediaries */
     7512        int        i,j,ig;
     7513        double     xyz_list[NUMVERTICES][3];
     7514        double     Jdet;
     7515        double     eps1dotdphii,eps1dotdphij;
     7516        double     eps2dotdphii,eps2dotdphij;
     7517        double     eps3dotdphii,eps3dotdphij;
     7518        double     mu_prime;
     7519        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     7520        double     eps1[3],eps2[3],eps3[3];
     7521        double     phi[NUMVERTICES];
     7522        double     dphi[3][NUMVERTICES];
     7523        GaussPenta *gauss=NULL;
     7524
     7525        /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
     7526        ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
     7527
     7528        /*Retrieve all inputs and parameters*/
     7529        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     7530        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     7531        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     7532        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
     7533
     7534        /* Start  looping on the number of gaussian points: */
     7535        gauss=new GaussPenta(5,5);
     7536        for (ig=gauss->begin();ig<gauss->end();ig++){
     7537
     7538                gauss->GaussPoint(ig);
     7539
     7540                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7541                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     7542
     7543                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7544                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     7545                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
     7546                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
     7547                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
     7548
     7549                for(i=0;i<6;i++){
     7550                        for(j=0;j<6;j++){
     7551                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
     7552                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
     7553                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
     7554                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
     7555                                eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
     7556                                eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
     7557
     7558                                Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     7559                                Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     7560                                Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
     7561
     7562                                Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     7563                                Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     7564                                Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
     7565
     7566                                Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
     7567                                Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
     7568                                Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
     7569                        }
     7570                }
     7571        }
     7572
     7573        /*Transform Coordinate System*/
     7574        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7575
     7576        /*Clean up and return*/
     7577        delete gauss;
     7578        return Ke;
     7579}
     7580/*}}}*/
    73217581/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    73227582void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){

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

    r11260 r11417  
    179179                void      GetPhi(double* phi, double*  epsilon, double viscosity);
    180180                void      GetSolutionFromInputsEnthalpy(Vec solutiong);
    181                 double  GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
     181                double  GetStabilizationParameter(double u, double v, double w, double diameter, double kappa);
    182182                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    183183                void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    … …  
    231231                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
    232232                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
     233<<<<<<< .working
     234=======
     235                ElementMatrix* CreateJacobianDiagnosticHoriz(void);
     236                ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void);
     237                ElementMatrix* CreateJacobianDiagnosticPattyn(void);
     238>>>>>>> .merge-right.r11410
     239                ElementMatrix* CreateJacobianDiagnosticStokes(void);
    233240                void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
    234241                void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);

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

    r11375 r11417  
    886886        delete gauss;
    887887        return pe;
     888}
     889/*}}}*/
     890/*FUNCTION Tria::CreateJacobianMatrix{{{1*/
     891void  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        }
    888921}
    889922/*}}}*/
    … …  
    30313064}
    30323065/*}}}*/
     3066/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/
     3067ElementMatrix* 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->GetViscosity2dDerivativeEpsSquare(&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/*}}}*/
    30333132/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    30343133void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
    … …  
    33323431        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    33333432
    3334         /*If on water, grad = 0: */
    3335         if(IsOnWater()) return;
     3433        /*If on water, grad = 0: (WHY???) */
     3434        if(IsOnWater()) return; 
    33363435
    33373436        /*First deal with ∂/∂alpha(KU-F)*/

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

    r11375 r11417  
    8484                void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
    8585                void   CreatePVector(Vec pf);
     86<<<<<<< .working
     87=======
     88                void   CreateJacobianMatrix(Mat Jff);
     89>>>>>>> .merge-right.r11410
    8690                int    GetNodeIndex(Node* node);
    8791                int    Sid();
    … …  
    206210                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    207211                ElementVector* CreatePVectorDiagnosticHutter(void);
     212                ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
    208213                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    209214                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);

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

    r11202 r11417  
    2424KML_File::KML_File(){
    2525
    26         kmlobj    =new DataSet;
     26        ;
    2727
    2828}
    … …  
    3131KML_File::~KML_File(){
    3232
    33         if (kmlobj) {
    34                 delete kmlobj;
    35                 kmlobj    =NULL;
    36         }
     33        ;
    3734
    3835}
    … …  
    4744        _printf_(flag,"KML_File:\n");
    4845        KML_Object::Echo();
    49 
    50         _printf_(flag,"        kmlobj: (size=%d)\n" ,kmlobj->Size());
    5146
    5247        return;
    … …  
    6661void  KML_File::DeepEcho(const char* indent){
    6762
    68         int   i;
    69         char  indent2[81];
    7063        bool  flag=true;
    7164
    … …  
    7366        KML_Object::DeepEcho(indent);
    7467
     68<<<<<<< .working
    7569/*  loop over the kml objects for the file  */
    7670
    … …  
    8882                _printf_(flag,"%s        kmlobj: [empty]\n"    ,indent);
    8983
     84=======
     85>>>>>>> .merge-right.r11410
    9086        return;
    9187}
    … …  
    9389/*FUNCTION KML_File::Write {{{1*/
    9490void  KML_File::Write(FILE* filout,const char* indent){
    95 
    96         int   i;
    97         char  indent2[81];
    9891
    9992        fprintf(filout,"%s<kml",indent);
    … …  
    10396
    10497        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);
    11498
    11599        fprintf(filout,"%s</kml>\n",indent);
    … …  
    145129                        _error_("KML_File::Read -- Unexpected field \"%s\".\n",kstri);
    146130
    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 
    249131                else if (!strncmp(kstri,"<",1))
    250132                        KML_Object::Read(fid,kstri);

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

    r11202 r11417  
    1919
    2020        public:
    21 
    22                 DataSet* kmlobj;
    2321
    2422                /*KML_File constructors, destructors {{{1*/

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

    r11202 r11417  
    2626        attrib    =new DataSet;
    2727        commnt    =new DataSet;
     28<<<<<<< .working
     29=======
     30        kmlobj    =new DataSet;
     31>>>>>>> .merge-right.r11410
    2832
    2933}
    … …  
    4044                commnt    =NULL;
    4145        }
     46<<<<<<< .working
     47=======
     48        if (kmlobj) {
     49                delete kmlobj;
     50                kmlobj    =NULL;
     51        }
     52>>>>>>> .merge-right.r11410
    4253
    4354}
    … …  
    5263        _printf_(flag,"        attrib: (size=%d)\n" ,attrib->Size());
    5364        _printf_(flag,"        commnt: (size=%d)\n" ,commnt->Size());
     65<<<<<<< .working
     66=======
     67        _printf_(flag,"        kmlobj: (size=%d)\n" ,kmlobj->Size());
     68>>>>>>> .merge-right.r11410
    5469
    5570        return;
    … …  
    7085
    7186        int   i;
     87        char  indent2[81];
    7288        bool  flag=true;
    7389
    … …  
    90106                _printf_(flag,"%s        commnt: [empty]\n"    ,indent);
    91107
     108<<<<<<< .working
     109=======
     110/*  loop over the unknown objects for the object  */
     111
     112        memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
     113        strcat(indent2,"  ");
     114
     115        if (kmlobj->Size())
     116                for (i=0; i<kmlobj->Size(); i++) {
     117            _printf_(flag,"%s        kmlobj: -------- begin [%d] --------\n" ,indent,i);
     118                        ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->DeepEcho(indent2);
     119            _printf_(flag,"%s        kmlobj: --------  end  [%d] --------\n" ,indent,i);
     120                }
     121        else
     122                _printf_(flag,"%s        kmlobj: [empty]\n"    ,indent);
     123
     124>>>>>>> .merge-right.r11410
    92125        return;
    93126}
    … …  
    101134        ;
    102135
     136<<<<<<< .working
     137=======
     138        memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char));
     139        strcat(indent2,"  ");
     140
     141        if (kmlobj->Size())
     142                for (i=0; i<kmlobj->Size(); i++) {
     143                        ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->Write(filout,indent2);
     144                }
     145
     146>>>>>>> .merge-right.r11410
    103147        return;
    104148}
    … …  
    116160                _error_("KML_Object::Read -- Unexpected field \"%s\".\n",kstr);
    117161
     162        else if (!strncmp(kstr,"<Placemark",10)) {
     163                kobj=(KML_Object*)new KML_Placemark();
     164                kobj->Read(fid,kstr);
     165                kmlobj    ->AddObject((Object*)kobj);
     166        }
     167
     168        else if (!strncmp(kstr,"<Folder", 7)) {
     169                kobj=(KML_Object*)new KML_Folder();
     170                kobj->Read(fid,kstr);
     171                kmlobj    ->AddObject((Object*)kobj);
     172        }
     173
     174        else if (!strncmp(kstr,"<Document", 9)) {
     175                kobj=(KML_Object*)new KML_Document();
     176                kobj->Read(fid,kstr);
     177                kmlobj    ->AddObject((Object*)kobj);
     178        }
     179
     180        else if (!strncmp(kstr,"<GroundOverlay",14)) {
     181                kobj=(KML_Object*)new KML_GroundOverlay();
     182                kobj->Read(fid,kstr);
     183                kmlobj    ->AddObject((Object*)kobj);
     184        }
     185
     186        else if (!strncmp(kstr,"<LatLonBox",10)) {
     187                kobj=(KML_Object*)new KML_LatLonBox();
     188                kobj->Read(fid,kstr);
     189                kmlobj    ->AddObject((Object*)kobj);
     190        }
     191
     192        else if (!strncmp(kstr,"<Icon", 5)) {
     193                kobj=(KML_Object*)new KML_Icon();
     194                kobj->Read(fid,kstr);
     195                kmlobj    ->AddObject((Object*)kobj);
     196        }
     197
     198        else if (!strncmp(kstr,"<Point", 6)) {
     199                kobj=(KML_Object*)new KML_Point();
     200                kobj->Read(fid,kstr);
     201                kmlobj    ->AddObject((Object*)kobj);
     202        }
     203
     204        else if (!strncmp(kstr,"<LineString",11)) {
     205                kobj=(KML_Object*)new KML_LineString();
     206                kobj->Read(fid,kstr);
     207                kmlobj    ->AddObject((Object*)kobj);
     208        }
     209
     210        else if (!strncmp(kstr,"<LinearRing",11)) {
     211                kobj=(KML_Object*)new KML_LinearRing();
     212                kobj->Read(fid,kstr);
     213                kmlobj    ->AddObject((Object*)kobj);
     214        }
     215
     216        else if (!strncmp(kstr,"<Polygon", 8)) {
     217                kobj=(KML_Object*)new KML_Polygon();
     218                kobj->Read(fid,kstr);
     219                kmlobj    ->AddObject((Object*)kobj);
     220        }
     221
     222        else if (!strncmp(kstr,"<MultiGeometry",14)) {
     223                kobj=(KML_Object*)new KML_MultiGeometry();
     224                kobj->Read(fid,kstr);
     225                kmlobj    ->AddObject((Object*)kobj);
     226        }
     227
     228//      else if (!strncmp(kstr,"<IconStyle",10)) {
     229//              kobj=(KML_Object*)new KML_IconStyle();
     230//              kobj->Read(fid,kstr);
     231//              kmlobj    ->AddObject((Object*)kobj);
     232//      }
     233
     234//      else if (!strncmp(kstr,"<LabelStyle",11)) {
     235//              kobj=(KML_Object*)new KML_LabelStyle();
     236//              kobj->Read(fid,kstr);
     237//              kmlobj    ->AddObject((Object*)kobj);
     238//      }
     239
     240        else if (!strncmp(kstr,"<LineStyle",10)) {
     241                kobj=(KML_Object*)new KML_LineStyle();
     242                kobj->Read(fid,kstr);
     243                kmlobj    ->AddObject((Object*)kobj);
     244        }
     245
     246        else if (!strncmp(kstr,"<PolyStyle",10)) {
     247                kobj=(KML_Object*)new KML_PolyStyle();
     248                kobj->Read(fid,kstr);
     249                kmlobj    ->AddObject((Object*)kobj);
     250        }
     251
     252//      else if (!strncmp(kstr,"<BalloonStyle",13)) {
     253//              kobj=(KML_Object*)new KML_BalloonStyle();
     254//              kobj->Read(fid,kstr);
     255//              kmlobj    ->AddObject((Object*)kobj);
     256//      }
     257
     258//      else if (!strncmp(kstr,"<ListStyle",10)) {
     259//              kobj=(KML_Object*)new KML_ListStyle();
     260//              kobj->Read(fid,kstr);
     261//              kmlobj    ->AddObject((Object*)kobj);
     262//      }
     263
    118264        else if (!strncmp(kstr,"<",1)) {
    119265                _printf_(true,"KML_Object::Read -- Unrecognized opening tag %s.\n",kstr);
     266<<<<<<< .working
    120267                KMLFileTagSkip(kstr,
    121268                                           fid);
     269=======
     270//              KMLFileTagSkip(kstr,
     271//                                         fid);
     272                kobj=(KML_Object*)new KML_Unknown();
     273                kobj->Read(fid,kstr);
     274                kmlobj    ->AddObject((Object*)kobj);
     275>>>>>>> .merge-right.r11410
    122276        }
    123277

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

    r11202 r11417  
    2121                DataSet* attrib;
    2222                DataSet* commnt;
     23<<<<<<< .working
     24=======
     25                DataSet* kmlobj;
     26>>>>>>> .merge-right.r11410
    2327
    2428                /*KML_Object constructors, destructors {{{1*/

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

    r11295 r11417  
    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"     );
    5153
    5254        return;
    … …  
    6668void  KML_Unknown::DeepEcho(const char* indent){
    6769
    68         bool  flag=true;
    69 
    70         _printf_(flag,"%sKML_Unknown_%s:\n",indent,name);
     70        char*        valuei;
     71        char*        vtoken;
     72        char         nl[]={'\n','\0'};
     73        bool         flag=true;
     74
     75        _printf_(flag,"%sKML_Unknown %s:\n",indent,name);
    7176        KML_Object::DeepEcho(indent);
    7277
    73         if (value     )
    74                 _printf_(flag,"%s         value: \"%s\"\n"     ,indent,value);
     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);
    7593
    7694        return;
    … …  
    8098void  KML_Unknown::Write(FILE* filout,const char* indent){
    8199
     100<<<<<<< .working
     101=======
     102        char*        valuei;
     103        char*        vtoken;
     104        char         nl[]={'\n','\0'};
     105
     106>>>>>>> .merge-right.r11410
    82107        fprintf(filout,"%s<%s",indent,name);
    83108        WriteAttrib(filout," ");
    … …  
    85110        WriteCommnt(filout,indent);
    86111
     112<<<<<<< .working
    87113        KML_Object::Write(filout,indent);
    88114
    89115        if (value     )
    90116                fprintf(filout,"%s  %s\n",indent,value);
    91 
     117=======
     118        if (value     ) {
     119                valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char));
     120                memcpy(valuei,value,(strlen(value)+1)*sizeof(char));
     121       
     122                vtoken=strtok(valuei,nl);
     123                fprintf(filout,"%s  %s\n",indent,vtoken);
     124   
     125                while (vtoken=strtok(NULL,nl))
     126                        fprintf(filout,"%s  %s\n",indent,vtoken);
     127>>>>>>> .merge-right.r11410
     128
     129<<<<<<< .working
     130=======
     131                xfree((void**)&valuei);
     132        }
     133
     134        KML_Object::Write(filout,indent);
     135
     136>>>>>>> .merge-right.r11410
    92137        fprintf(filout,"%s</%s>\n",indent,name);
    93138
    … …  
    101146        int          ncom=0;
    102147        char**       pcom=NULL;
     148        char         nl[]={'\n','\0'};
    103149
    104150/*  get object attributes and check for solo tag  */
    … …  
    120166                        _error_("KML_Unknown::Read -- Unexpected closing tag %s.\n",kstri);
    121167
     168<<<<<<< .working
    122169                else if (strncmp(kstri,"<",1))
    123170                        KMLFileTokenParse( value     ,NULL,0,
    124171                                                          NULL,
    125172                                                          fid);
     173=======
     174                else if (strncmp(kstri,"<",1)) {
     175                        if (value) {
     176                                value=(char *) xrealloc(value,(strlen(value)+1+strlen(kstri)+1)*sizeof(char));
     177                                strcat(value,nl);
     178                                strcat(value,kstri);
     179                        }
     180                        else {
     181                                value=(char *) xmalloc((strlen(kstri)+1)*sizeof(char));
     182                                memcpy(value,kstri,(strlen(kstri)+1)*sizeof(char));
     183                        }
     184                }
     185>>>>>>> .merge-right.r11410
    126186
    127187                else if (!strncmp(kstri,"<",1))

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

    r10576 r11417  
    361361}
    362362/*}}}*/
     363/*FUNCTION Icefront::CreateJacobianMatrix{{{1*/
     364void  Icefront::CreateJacobianMatrix(Mat Jff){
     365        this->CreateKMatrix(Jff,NULL);
     366}
     367/*}}}1*/
    363368/*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/
    364369void  Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){
    … …  
    373378}
    374379/*}}}*/
     380/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/
     381void  Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
     382        this->PenaltyCreateKMatrix(Jff,NULL,kmax);
     383}
     384/*}}}1*/
    375385/*FUNCTION Icefront::InAnalysis{{{1*/
    376386bool Icefront::InAnalysis(int in_analysis_type){

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

    r10576 r11417  
    7676                void  CreateKMatrix(Mat Kff, Mat Kfs);
    7777                void  CreatePVector(Vec pf);
     78<<<<<<< .working
     79=======
     80                void  CreateJacobianMatrix(Mat Jff);
     81>>>>>>> .merge-right.r11410
    7882                void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
    7983                void  PenaltyCreatePVector(Vec pf, double kmax);
     84                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
    8085                bool  InAnalysis(int analysis_type);
    8186                /*}}}*/

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

    r11247 r11417  
    214214}
    215215/*}}}1*/
     216<<<<<<< .working
     217=======
     218/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
     219void  Penpair::CreateJacobianMatrix(Mat Jff){
     220        this->CreateKMatrix(Jff,NULL);
     221}
     222/*}}}1*/
     223>>>>>>> .merge-right.r11410
    216224/*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/
    217225void  Penpair::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){

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

    r11373 r11417  
    606606        /*Return: */
    607607        *pviscosity_complement=viscosity_complement;
     608}
     609/*}}}*/
     610/*FUNCTION Matice::GetViscosity2dDerivativeEpsSquare{{{1*/
     611void  Matice::GetViscosity2dDerivativeEpsSquare(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;
    608639}
    609640/*}}}*/

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

    r11373 r11417  
    6969                void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
    7070                void   GetViscosityZComplement(double* pviscosity_complement, double* pepsilon);
     71                void   GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
    7172                double GetB();
    7273                double GetBbar();

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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 a string",enum_type,EnumToStringx(enum_type));};
     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));};
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",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));}
    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

    r11202 r11417  
    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 arrayl",enum_type,EnumToStringx(enum_type));}
     60                void  GetParameterValue(char*** pstringarray,int* pM){_error_("String param of enum %i (%s) cannot return a string array",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

    r9320 r11417  
    1212#include <stdio.h>
    1313#include <string.h>
     14#include <math.h>
    1415#include "./objects.h"
    1516#include "../Container/Container.h"

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

    r9320 r11417  
    165165
    166166        verbositylevel = (int)level;
     167        return verbositylevel;
    167168
    168169#else

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

    r9320 r11417  
    6464
    6565        verbositylevel = (int)level;
     66        return verbositylevel;
    6667
    6768#else

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

    r9320 r11417  
    124124******************************************************************************************************************************/
    125125
    126 int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
     126void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
    127127       
    128128        int i,counter;
    … …  
    379379                                   IsInPoly
    380380******************************************************************************************************************************/
    381 int 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 }
     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//}
    400400
    401401/******************************************************************************************************************************

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

    r9320 r11417  
    2626int IsNeighbor(int el1,int el2,double* index);
    2727int IsOnRift(int el,int nriftsegs,int* riftsegments);
    28 int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments);
     28void 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

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

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

    r11297 r11417  
    2525
    2626        /*TAO*/
     27<<<<<<< .working
    2728        int       ierr,numberofvertices;
    2829        AppCtx    user;
    … …  
    3132        Vec       XL = NULL;
    3233        Vec       XU = NULL;
     34=======
     35        int        ierr;
     36        int        num_controls,solution_type;
     37        int        nsteps,maxiter;
     38        AppCtx     user;
     39        TaoSolver  tao;
     40        double    *dummy          = NULL;
     41        int       *control_list   = NULL;
     42        Vec        X              = NULL;
     43        Vec        XL             = NULL;
     44        Vec        XU             = NULL;
     45>>>>>>> .merge-right.r11410
    3346
    3447        /*Initialize TAO*/
    … …  
    3851        if(ierr) _error_("Could not initialize Tao");
    3952
     53<<<<<<< .working
     54=======
     55        /*Recover some parameters*/
     56        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     57        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     58        femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum);
     59        femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
     60        femmodel->parameters->FindParam(&dummy,NULL,NULL,InversionMaxiterPerStepEnum);
     61        maxiter=nsteps*(int)dummy[0]; xfree((void**)&dummy);
     62
     63>>>>>>> .merge-right.r11410
    4064        /*Initialize TAO*/
    4165        TaoCreate(PETSC_COMM_WORLD,&tao);
    … …  
    4569
    4670        /*Prepare all TAO parameters*/
    47         TaoSetMaximumFunctionEvaluations(tao,50);
    48         TaoSetMaximumIterations(tao,10);
     71        TaoSetMaximumFunctionEvaluations(tao,maxiter);
     72        TaoSetMaximumIterations(tao,nsteps);
    4973        TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28);
    5074
    … …  
    6690        InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
    6791        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
     92
     93        /*Finalize*/
     94        _printf_(VerboseControl(),"%s\n","   preparing final solution");
     95        femmodel->parameters->SetParam(false,InversionIscontrolEnum); //needed to turn control result output in solutioncore
     96        void (*solutioncore)(FemModel*)=NULL;
     97        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     98        solutioncore(femmodel);
    6899
    69100        /*Clean up and return*/
    … …  
    90121        VecFree(&gradient);
    91122        CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     123<<<<<<< .working
     124=======
     125
     126        /*Compute gradient*/
     127        Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     128        VecCopy(gradient,G); VecFree(&gradient);
     129        VecScale(G,-1.);
     130
     131        /*Clean-up and return*/
     132        xfree((void**)&cost_functions);
     133        xfree((void**)&cost_functionsd);
     134>>>>>>> .merge-right.r11410
    92135        return 0;
    93136}

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

    r11245 r11417  
    2424        /*parameters: */
    2525        double finaltime,dt,yts;
    26         bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline;
     26        bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy;
    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);
    5354        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5455        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    … …  
    9192                        _printf_(VerboseSolution(),"   computing temperatures:\n");
    9293                        #ifdef _HAVE_THERMAL_
    93                         thermal_core_step(femmodel,step,time);
     94                        if(isenthalpy==0){
     95                                thermal_core_step(femmodel,step,time);
     96                        }
     97                        else{
     98                                enthalpy_core_step(femmodel,step,time);
     99                        }
    94100                        #else
    95101                        _error_("ISSM was not compiled with thermal capabilities. Exiting");
    … …  
    135141                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
    136142                        if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
    137                         InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,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);
    138146                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time);
    139147                        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

    r11293 r11417  
    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);
    8790                        break;
    8891                }

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

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

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

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

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

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

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

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

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

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

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

    r11219 r11417  
    66classdef solver
    77    properties (SetAccess=public)
    8                  options={NoneAnalysisEnum,mumpsoptions};
     8                 options=cell(0,0);
    99         end
    1010         methods
    … …  
    2626                 function obj = setdefaultparameters(obj) % {{{
    2727
     28                         %MUMPS is the default solver
     29                         obj.options={'NoneAnalysis',mumpsoptions};
     30
    2831                 end % }}}
    29                  function obj=addoptions(obj,analysis,solveroptions) % {{{1
     32                 function obj = addoptions(obj,analysis,solveroptions) % {{{1
     33
     34                         %Convert analysis from enum to string
     35                         analysis=EnumToString(analysis);
     36
    3037                         %first, find out if analysis has already been supplied
    3138                         found=false;
    3239                         for i=1:size(obj.options,1),
    3340                                 inanalysis=obj.options{i,1};
    34                                  if inanalysis==analysis,
     41                                 if strcmp(inanalysis,analysis),
    3542                                         found=true;
    36                                          obj.options{i,1}=analysis;
    37                                          obj.options{i,2}=solveroptions;
     43                                         obj.options{i,1} = analysis;
     44                                         obj.options{i,2} = solveroptions;
    3845                                         break;
    3946                                 end
    4047                         end
     48
    4149                         if ~found,
    42                                  obj.options{end+1,1}=analysis;
    43                                  obj.options{end,2}=solveroptions;
     50                                 obj.options{end+1,1}= analysis;
     51                                 obj.options{end,2}  = solveroptions;
    4452                         end
    4553                 end
    4654                 %}}}
    4755                 function checkconsistency(obj,md,solution,analyses) % {{{
    48 
     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
    4961                 end % }}}
    5062                 function PetscFile(solver,filename) % {{{
    … …  
    7082
    7183                                 %first write analysis:
    72                                  fprintf(fid,'\n+%s\n',EnumToString(analysis)); %append a + to recognize it's an analysis enum
     84                                 fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
    7385
    7486                                 %now, write options
    … …  
    126138                                end
    127139
    128                                 disp(sprintf('   %s -> ''%s''',EnumToString(analysis),string));
     140                                disp(sprintf('   %s -> ''%s''',analysis,string));
    129141                        end
    130142                 end

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

    r11265 r11417  
    1212                penalty_lock      = 0;
    1313                penalty_factor    = 0;
     14                isenthalpy        = 0;
    1415        end
    1516        methods
    … …  
    4243                        %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
    4344                        obj.penalty_factor=3;
     45
     46                        %Should we use cold ice (default) or enthalpy formulation
     47                        obj.isenthalpy=0;
    4448                end % }}}
    4549                function checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5054                        checkfield(md,'thermal.stabilization','numel',1,'values',[0 1 2]);
    5155                        checkfield(md,'thermal.spctemperature','forcing',1);
    52                         if ismember(EnthalpyAnalysisEnum,analyses),
     56                        if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy),
    5357                        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]);
    5459                        end
    5560                end % }}}
    … …  
    6267                        fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
    6368                        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)');
    6470
    6571                end % }}}
    … …  
    7177                        WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
    7278                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     79                        WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
    7380                end % }}}
    7481        end

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

    r11368 r11417  
    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;
     31if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
     32if ~isnan(md.inversion.min_parameters), md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
     33if ~isnan(md.inversion.max_parameters), md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
    3134if ~isnan(md.surfaceforcings.mass_balance),
    3235        md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers);
    … …  
    4952if ~isnan(md.flowequation.element_equation)
    5053        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);
    5158end     
    5259
    … …  
    5562md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers);
    5663md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers);
     64md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers);
    5765md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers);
    5866md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
    … …  
    6068%Extrusion of Neumann BC
    6169if ~isnan(md.diagnostic.icefront),
    62         numberofneumann2d=size(md.diagnostic.icefront,1)/md.mesh.numberoflayers;
     70        numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1);
    6371        md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
    6472end
    … …  
    8997md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1);
    9098md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1);
     99md.mask.elementonwater=project2d(md,md.mask.elementonwater,1);
     100md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1);
    91101
    92102%Initialize with the 2d mesh

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

    r10538 r11417  
    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
    1419
    1520        %initialize md.results if not a structure yet

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

    r10683 r11417  
    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/');
     13%setenv('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

    r9762 r11417  
    1010if strcmpi(string,'diagnostic'),
    1111        if md.mesh.dimension==2,
    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
     12                md.initialization.vx=md.results.DiagnosticSolution.Vx;
     13                md.initialization.vy=md.results.DiagnosticSolution.Vy;
    2214        else
    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
     15                md.initialization.vx=md.results.DiagnosticSolution.Vx;
     16                md.initialization.vy=md.results.DiagnosticSolution.Vy;
     17                md.initialization.vz=md.results.DiagnosticSolution.Vz;
    3018        end
    31         md.initialization.vel=PatchToVec(md.results.DiagnosticSolution.Vel);
     19        md.initialization.vel=md.results.DiagnosticSolution.Vel;
    3220
    3321        if isfield(md.results.DiagnosticSolution,'Pressure'),
    34                 md.initialization.pressure=PatchToVec(md.results.DiagnosticSolution.Pressure);
     22                md.initialization.pressure=md.results.DiagnosticSolution.Pressure;
    3523        end
    3624        if md.rifts.numrifts,
    … …  
    4230                for control_parameters=md.inversion.control_parameters
    4331                        %Will need to be updated... good luck ;)
    44                         md.(EnumToModelField(control_parameters))=PatchToVec(md.results.DiagnosticSolution.(EnumToString(control_parameters)));
     32                        md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters));
    4533                end
    4634        end
    … …  
    5947        for i=1:length(results),
    6048                if ~isempty(md.results.TransientSolution(i).Vel),
    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);
     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;
    6755                        results2(count).time=md.results.TransientSolution(i).time;
    6856                        results2(count).step=md.results.TransientSolution(i).step;
    … …  
    7664        clear results,results2;
    7765elseif strcmpi(string,'steadystate'),
    78         md.initialization.vx=PatchToVec(md.results.SteadystateSolution.Vx);
    79         md.initialization.vy=PatchToVec(md.results.SteadystateSolution.Vy);
     66        md.initialization.vx=md.results.SteadystateSolution.Vx;
     67        md.initialization.vy=md.results.SteadystateSolution.Vy;
    8068        if isfield(md.results.SteadystateSolution,'Vz'),
    81                 md.initialization.vz=PatchToVec(md.results.SteadystateSolution.Vz);
     69                md.initialization.vz=md.results.SteadystateSolution.Vz;
    8270        end
    8371
    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);
     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;
    8876
    8977        if md.inversion.iscontrol==1,
    9078                for control_parameters=md.inversion.control_parameters
    91                         md.(EnumToModelField(control_parameters))=PatchToVec(md.results.SteadystateSolution.(EnumToString(control_parameters)));
     79                        md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters));
    9280                end
    9381        end
    9482
    9583elseif strcmpi(string,'thermal'),
    96         md.initialization.temperature=PatchToVec(md.results.ThermalSolution.Temperature);
    97         md.basalforcings.melting_rate=PatchToVec(md.results.ThermalSolution.BasalMeltingRate);
     84        md.initialization.temperature=md.results.ThermalSolution.Temperature;
     85        md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate;
    9886elseif strcmpi(string,'hydrology'),
    99         md.initialization.watercolumn=PatchToVec(md.results.HydrologySolution.Watercolumn);
     87        md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn;
    10088
    10189else

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

    r9668 r11417  
    11function dresp=expandresponses(md,responses)
     2%EXPANDRESPONSES - expand responses
    23
    34fnames=fieldnames(responses);

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

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

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

    r10649 r11417  
    1515        ismacayealpattyn=femmodel.parameters.FlowequationIsmacayealpattyn;
    1616        isstokes=femmodel.parameters.FlowequationIsstokes;
     17        isnewton=femmodel.parameters.DiagnosticIsnewton;
    1718        dakota_analysis=femmodel.parameters.QmuIsdakota;
    1819        control_analysis=femmodel.parameters.InversionIscontrol;
    … …  
    5354                issmprintf(VerboseSolution,'\n%s',['   computing horizontal velocities']);
    5455                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    55                 femmodel=solver_nonlinear(femmodel,modify_loads);
     56                if isnewton,
     57                        femmodel=solver_newton(femmodel);
     58                else
     59                        femmodel=solver_nonlinear(femmodel,modify_loads);
     60                end
    5661        end
    5762       

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

    r10999 r11417  
    1919        isthermal=femmodel.parameters.TransientIsthermal;
    2020        isgroundingline=femmodel.parameters.TransientIsgroundingline;
     21        isenthalpy=femmodel.parameters.ThermalIsenthalpy;
    2122        groundinglinemigration=femmodel.parameters.GroundinglineMigration;
    2223
    … …  
    5758                if (isthermal & dim==3)
    5859                        issmprintf(VerboseSolution,'\n%s',['   computing temperature']);
    59                         femmodel=thermal_core_step(femmodel);
     60                        if (isenthalpy==0),
     61                                femmodel=thermal_core_step(femmodel);
     62                        else
     63                                femmodel=enthalpy_core_step(femmodel);
     64                        end
    6065                end
    6166
    … …  
    9095                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time);
    9196                        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
    9299                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time);
    93100                        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

    r11295 r11417  
    2222                                CostFunction \
    2323                                CreateNodalConstraints\
     24                                CreateJacobianMatrix\
    2425                                Echo\
    2526                                ElementConnectivity\
    … …  
    327328                          SystemMatrices/SystemMatrices.h
    328329
     330
     331CreateJacobianMatrix_SOURCES = CreateJacobianMatrix/CreateJacobianMatrix.cpp\
     332                                                                 CreateJacobianMatrix/CreateJacobianMatrix.h
     333
    329334SurfaceArea_SOURCES = SurfaceArea/SurfaceArea.cpp\
    330335                                                                 SurfaceArea/SurfaceArea.h
Note: See TracChangeset for help on using the changeset viewer.