Index: /issm/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 11417) @@ -162,4 +162,5 @@ ThermalSpctemperatureEnum, ThermalStabilizationEnum, + ThermalIsenthalpyEnum, ThicknessEnum, TimesteppingCflCoefficientEnum, Index: /issm/trunk-jpl-damage/src/c/EnumDefinitions/Synchronize.sh =================================================================== --- /issm/trunk-jpl-damage/src/c/EnumDefinitions/Synchronize.sh (revision 11416) +++ /issm/trunk-jpl-damage/src/c/EnumDefinitions/Synchronize.sh (revision 11417) @@ -9,4 +9,7 @@ rm $ISSM_TIER/src/c/modules/EnumToStringx/EnumToStringx.cpp rm $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp + +#Get number of enums +NUMENUMS=$(wc -l temp | awk '{printf("%s",$1);}'); #Take care of EnumToModelField.m first (easy) @@ -95,20 +98,33 @@ int StringToEnumx(const char* name){ + int stage=1; + END + #core -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 -#Footer +i1=1; +i2=120; +for (( i=1 ; i<=100 ; i++ )); do + echo " if(stage==$i){" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp + awk -v i1=$i1 -v i2=$i2 '{if(NR>=i1 && NR<=i2) print $0 }' 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 + echo " else stage=$(($i+1));" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp + echo " }" >> $ISSM_TIER//src/c/modules/StringToEnumx/StringToEnumx.cpp + + if [ $i2 -ge $NUMENUMS ]; then break; fi + let i1=$i1+120 + let i2=$i2+120 +done + +#footer cat <> $ISSM_TIER/src/c/modules/StringToEnumx/StringToEnumx.cpp - else _error_("Enum %s not found",name); - + /*If we reach this point, the string provided has not been found*/ + _error_("Enum %s not found",name); } END #}}} -#get number of lines in temp -NUMBEROFLINES=$(wc -l temp | awk '{printf("%s",$1);}'); - # go through the lines of temp -for (( i=1 ; i<=$NUMBEROFLINES ; i++ )); do +for (( i=1 ; i<=$NUMENUMS ; i++ )); do #Get name and enum of the line i @@ -123,13 +139,13 @@ then printf "\r " - printf "\r $i/$NUMBEROFLINES Adding "$NAME"..." + printf "\r $i/$NUMENUMS Adding "$NAME"..." else if [ $i -lt 100 ] then printf "\r " - printf "\r $i/$NUMBEROFLINES Adding "$NAME"..." + printf "\r $i/$NUMENUMS Adding "$NAME"..." else printf "\r " - printf "\r$i/$NUMBEROFLINES Adding "$NAME"..." + printf "\r$i/$NUMENUMS Adding "$NAME"..." fi fi @@ -153,5 +169,4 @@ done - #clean up{{{ rm temp Index: /issm/trunk-jpl-damage/src/c/modules/BamgTriangulatex/BamgTriangulatex.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/BamgTriangulatex/BamgTriangulatex.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/BamgTriangulatex/BamgTriangulatex.cpp (revision 11417) @@ -17,4 +17,5 @@ Th.WriteIndex(pindex,pnels); //delete &Th; + return 0; } Index: /issm/trunk-jpl-damage/src/c/modules/Chacox/Chacox.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/Chacox/Chacox.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/Chacox/Chacox.cpp (revision 11417) @@ -191,4 +191,6 @@ - #endif //ifdef _HAVE_CHACO_ + #else //ifdef _HAVE_CHACO_ + return (0); + #endif } Index: /issm/trunk-jpl-damage/src/c/modules/Chacox/input_parse.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/Chacox/input_parse.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/Chacox/input_parse.cpp (revision 11417) @@ -262,4 +262,6 @@ return(0); - #endif //#ifdef _HAVE_CHACO_ + #else //#ifdef _HAVE_CHACO_ + return(0); + #endif } Index: /issm/trunk-jpl-damage/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp (revision 11417) @@ -69,3 +69,5 @@ } + return NULL; + } Index: /issm/trunk-jpl-damage/src/c/modules/ContourToNodesx/ContourToNodesx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/ContourToNodesx/ContourToNodesx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/ContourToNodesx/ContourToNodesx.cpp (revision 11417) @@ -36,4 +36,5 @@ /*Assign output pointers: */ *pflags=flags; - + + return 1; } Index: /issm/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11417) @@ -166,4 +166,5 @@ case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; case ThermalStabilizationEnum : return "ThermalStabilization"; + case ThermalIsenthalpyEnum : return "ThermalIsenthalpy"; case ThicknessEnum : return "Thickness"; case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient"; Index: /issm/trunk-jpl-damage/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp (revision 11417) @@ -126,4 +126,6 @@ } if (debug && my_thread==0) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + return NULL; } Index: /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 11417) @@ -78,4 +78,5 @@ parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum)); parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum)); + parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum)); parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum)); Index: /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 11417) @@ -20,5 +20,5 @@ int i,analysis_type,dim,verbose; - bool isthermal,isprognostic,isdiagnostic,isgroundingline; + bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy; /*output: */ @@ -39,4 +39,5 @@ iomodel->Constant(&verbose,VerboseEnum); iomodel->Constant(&isthermal,TransientIsthermalEnum); + iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum); iomodel->Constant(&isprognostic,TransientIsprognosticEnum); iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum); @@ -52,6 +53,11 @@ if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && dim==2) continue; if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && dim==2) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && dim==2) continue; if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isthermal==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isthermal==false) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isthermal==false) continue; + if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isenthalpy==true) continue; + if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isenthalpy==true) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isenthalpy==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==PrognosticAnalysisEnum && isprognostic==false && isgroundingline==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==DiagnosticHorizAnalysisEnum && isdiagnostic==false) continue; Index: /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp (revision 11417) @@ -36,3 +36,5 @@ /*Assign output pointers: */ *pflags=flags; + + return 1; } Index: /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp (revision 11417) @@ -74,3 +74,5 @@ /*Free ressources:*/ xfree((void**)&already); + + return NULL; } Index: /issm/trunk-jpl-damage/src/c/modules/Scotchx/Scotchx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/Scotchx/Scotchx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/Scotchx/Scotchx.cpp (revision 11417) @@ -339,4 +339,6 @@ return (0); -#endif //#ifdef _HAVE_SCOTCH_ +#else //#ifdef _HAVE_SCOTCH_ + return(0); +#endif } Index: /issm/trunk-jpl-damage/src/c/modules/Shp2Kmlx/Shp2Kmlx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/Shp2Kmlx/Shp2Kmlx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/Shp2Kmlx/Shp2Kmlx.cpp (revision 11417) @@ -24,5 +24,7 @@ sgn,cm,sp)); - #endif //ifdef _HAVE_SHAPELIB_ + #else //ifdef _HAVE_SHAPELIB_ + return 0; + #endif } @@ -607,5 +609,7 @@ return(iret); - #endif //ifdef _HAVE_SHAPELIB_ + #else //ifdef _HAVE_SHAPELIB_ + return 0; + #endif } Index: /issm/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11417) @@ -14,4 +14,5 @@ int StringToEnumx(const char* name){ +<<<<<<< .working if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum; else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum; @@ -443,4 +444,449 @@ else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else _error_("Enum %s not found",name); +======= + int stage=1; +>>>>>>> .merge-right.r11410 + if(stage==1){ + if (strcmp(name,"AutodiffForward")==0) return AutodiffForwardEnum; + else if (strcmp(name,"AutodiffIsautodiff")==0) return AutodiffIsautodiffEnum; + else if (strcmp(name,"AutodiffReverse")==0) return AutodiffReverseEnum; + else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum; + else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum; + else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum; + else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; + else if (strcmp(name,"BasalforcingsMeltingRateCorrection")==0) return BasalforcingsMeltingRateCorrectionEnum; + else if (strcmp(name,"BasalforcingsMeltingRate")==0) return BasalforcingsMeltingRateEnum; + else if (strcmp(name,"Bathymetry")==0) return BathymetryEnum; + else if (strcmp(name,"Bed")==0) return BedEnum; + else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum; + else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum; + else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum; + else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum; + else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum; + else if (strcmp(name,"DiagnosticIsnewton")==0) return DiagnosticIsnewtonEnum; + else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum; + else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum; + else if (strcmp(name,"DiagnosticPenaltyFactor")==0) return DiagnosticPenaltyFactorEnum; + else if (strcmp(name,"DiagnosticReferential")==0) return DiagnosticReferentialEnum; + else if (strcmp(name,"DiagnosticReltol")==0) return DiagnosticReltolEnum; + else if (strcmp(name,"DiagnosticRequestedOutputs")==0) return DiagnosticRequestedOutputsEnum; + else if (strcmp(name,"DiagnosticRestol")==0) return DiagnosticRestolEnum; + else if (strcmp(name,"DiagnosticRiftPenaltyLock")==0) return DiagnosticRiftPenaltyLockEnum; + else if (strcmp(name,"DiagnosticRiftPenaltyThreshold")==0) return DiagnosticRiftPenaltyThresholdEnum; + else if (strcmp(name,"DiagnosticShelfDampening")==0) return DiagnosticShelfDampeningEnum; + else if (strcmp(name,"DiagnosticSpcvx")==0) return DiagnosticSpcvxEnum; + else if (strcmp(name,"DiagnosticSpcvy")==0) return DiagnosticSpcvyEnum; + else if (strcmp(name,"DiagnosticSpcvz")==0) return DiagnosticSpcvzEnum; + else if (strcmp(name,"DiagnosticStokesreconditioning")==0) return DiagnosticStokesreconditioningEnum; + else if (strcmp(name,"DiagnosticVertexPairing")==0) return DiagnosticVertexPairingEnum; + else if (strcmp(name,"DiagnosticViscosityOvershoot")==0) return DiagnosticViscosityOvershootEnum; + else if (strcmp(name,"FlowequationBordermacayeal")==0) return FlowequationBordermacayealEnum; + else if (strcmp(name,"FlowequationBorderpattyn")==0) return FlowequationBorderpattynEnum; + else if (strcmp(name,"FlowequationBorderstokes")==0) return FlowequationBorderstokesEnum; + else if (strcmp(name,"FlowequationElementEquation")==0) return FlowequationElementEquationEnum; + else if (strcmp(name,"FlowequationIshutter")==0) return FlowequationIshutterEnum; + else if (strcmp(name,"FlowequationIsmacayealpattyn")==0) return FlowequationIsmacayealpattynEnum; + else if (strcmp(name,"FlowequationIsstokes")==0) return FlowequationIsstokesEnum; + else if (strcmp(name,"FlowequationVertexEquation")==0) return FlowequationVertexEquationEnum; + else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum; + else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; + else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; + else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; + else if (strcmp(name,"HydrologyCR")==0) return HydrologyCREnum; + else if (strcmp(name,"HydrologyKn")==0) return HydrologyKnEnum; + else if (strcmp(name,"HydrologyN")==0) return HydrologyNEnum; + else if (strcmp(name,"HydrologyP")==0) return HydrologyPEnum; + else if (strcmp(name,"HydrologyQ")==0) return HydrologyQEnum; + else if (strcmp(name,"HydrologySpcwatercolumn")==0) return HydrologySpcwatercolumnEnum; + else if (strcmp(name,"HydrologyStabilization")==0) return HydrologyStabilizationEnum; + else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; + else if (strcmp(name,"InversionCostFunction")==0) return InversionCostFunctionEnum; + else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum; + else if (strcmp(name,"InversionGradientOnly")==0) return InversionGradientOnlyEnum; + else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum; + else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum; + else if (strcmp(name,"InversionTao")==0) return InversionTaoEnum; + else if (strcmp(name,"InversionMaxParameters")==0) return InversionMaxParametersEnum; + else if (strcmp(name,"InversionMaxiterPerStep")==0) return InversionMaxiterPerStepEnum; + else if (strcmp(name,"InversionMinParameters")==0) return InversionMinParametersEnum; + else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum; + else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum; + else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum; + else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; + else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; + else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; + else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; + else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; + else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum; + else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum; + else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum; + else if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum; + else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum; + else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum; + else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; + else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum; + else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; + else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; + else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum; + else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; + else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum; + else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; + else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; + else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum; + else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; + else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; + else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; + else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum; + else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum; + else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum; + else if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum; + else if (strcmp(name,"MeshElementonbed")==0) return MeshElementonbedEnum; + else if (strcmp(name,"MeshElementonsurface")==0) return MeshElementonsurfaceEnum; + else if (strcmp(name,"MeshElements2d")==0) return MeshElements2dEnum; + else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; + else if (strcmp(name,"MeshLowerelements")==0) return MeshLowerelementsEnum; + else if (strcmp(name,"MeshNumberofedges")==0) return MeshNumberofedgesEnum; + else if (strcmp(name,"MeshNumberofelements2d")==0) return MeshNumberofelements2dEnum; + else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum; + else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum; + else if (strcmp(name,"MeshNumberofvertices2d")==0) return MeshNumberofvertices2dEnum; + else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum; + else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum; + else if (strcmp(name,"MeshVertexonbed")==0) return MeshVertexonbedEnum; + else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; + else if (strcmp(name,"MeshX")==0) return MeshXEnum; + else if (strcmp(name,"MeshY")==0) return MeshYEnum; + else if (strcmp(name,"MeshZ")==0) return MeshZEnum; + else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum; + else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum; + else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum; + else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum; + else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum; + else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum; + else stage=2; + } + if(stage==2){ + if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum; + else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; + else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum; + else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum; + else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum; + else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum; + else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum; + else if (strcmp(name,"QmuVariabledescriptors")==0) return QmuVariabledescriptorsEnum; + else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum; + else if (strcmp(name,"RiftsRiftstruct")==0) return RiftsRiftstructEnum; + else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; + else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum; + else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; + else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum; + else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; + else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; + else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; + else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; + else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; + else if (strcmp(name,"Surface")==0) return SurfaceEnum; + else if (strcmp(name,"SurfaceforcingsAblationRate")==0) return SurfaceforcingsAblationRateEnum; + else if (strcmp(name,"SurfaceforcingsAccumulationRate")==0) return SurfaceforcingsAccumulationRateEnum; + else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; + else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; + else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; + else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; + else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum; + else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; + else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum; + else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum; + else if (strcmp(name,"Thickness")==0) return ThicknessEnum; + else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; + else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum; + else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum; + else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum; + else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum; + else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; + else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum; + else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; + else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; + else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; + else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; + else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum; + else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum; + else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; + else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum; + else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; + else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; + else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; + else if (strcmp(name,"BedSlopeAnalysis")==0) return BedSlopeAnalysisEnum; + else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; + else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum; + else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum; + else if (strcmp(name,"DiagnosticHorizAnalysis")==0) return DiagnosticHorizAnalysisEnum; + else if (strcmp(name,"DiagnosticHutterAnalysis")==0) return DiagnosticHutterAnalysisEnum; + else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum; + else if (strcmp(name,"DiagnosticVertAnalysis")==0) return DiagnosticVertAnalysisEnum; + else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; + else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum; + else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum; + else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum; + else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum; + else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; + else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; + else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum; + else if (strcmp(name,"PrognosticAnalysis")==0) return PrognosticAnalysisEnum; + else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum; + else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; + else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum; + else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; + else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum; + else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum; + else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; + else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; + else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + else if (strcmp(name,"Approximation")==0) return ApproximationEnum; + else if (strcmp(name,"HutterApproximation")==0) return HutterApproximationEnum; + else if (strcmp(name,"MacAyealApproximation")==0) return MacAyealApproximationEnum; + else if (strcmp(name,"MacAyealPattynApproximation")==0) return MacAyealPattynApproximationEnum; + else if (strcmp(name,"MacAyealStokesApproximation")==0) return MacAyealStokesApproximationEnum; + else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; + else if (strcmp(name,"PattynApproximation")==0) return PattynApproximationEnum; + else if (strcmp(name,"PattynStokesApproximation")==0) return PattynStokesApproximationEnum; + else if (strcmp(name,"StokesApproximation")==0) return StokesApproximationEnum; + else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; + else if (strcmp(name,"Loads")==0) return LoadsEnum; + else if (strcmp(name,"Materials")==0) return MaterialsEnum; + else if (strcmp(name,"Nodes")==0) return NodesEnum; + else if (strcmp(name,"Parameters")==0) return ParametersEnum; + else if (strcmp(name,"Vertices")==0) return VerticesEnum; + else if (strcmp(name,"Results")==0) return ResultsEnum; + else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; + else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; + else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; + else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum; + else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; + else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; + else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; + else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; + else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; + else if (strcmp(name,"Element")==0) return ElementEnum; + else if (strcmp(name,"ElementResult")==0) return ElementResultEnum; + else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; + else if (strcmp(name,"FileParam")==0) return FileParamEnum; + else if (strcmp(name,"Hook")==0) return HookEnum; + else if (strcmp(name,"Icefront")==0) return IcefrontEnum; + else if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"IntInput")==0) return IntInputEnum; + else if (strcmp(name,"IntParam")==0) return IntParamEnum; + else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; + else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum; + else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum; + else if (strcmp(name,"Matice")==0) return MaticeEnum; + else if (strcmp(name,"Matpar")==0) return MatparEnum; + else if (strcmp(name,"Node")==0) return NodeEnum; + else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + else if (strcmp(name,"Param")==0) return ParamEnum; + else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"Pengrid")==0) return PengridEnum; + else if (strcmp(name,"Penpair")==0) return PenpairEnum; + else if (strcmp(name,"Penta")==0) return PentaEnum; + else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum; + else if (strcmp(name,"PetscMatParam")==0) return PetscMatParamEnum; + else if (strcmp(name,"PetscVecParam")==0) return PetscVecParamEnum; + else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + else if (strcmp(name,"Segment")==0) return SegmentEnum; + else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; + else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; + else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; + else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; + else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum; + else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; + else if (strcmp(name,"StringParam")==0) return StringParamEnum; + else if (strcmp(name,"Tria")==0) return TriaEnum; + else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum; + else if (strcmp(name,"Vertex")==0) return VertexEnum; + else if (strcmp(name,"Air")==0) return AirEnum; + else if (strcmp(name,"Ice")==0) return IceEnum; + else if (strcmp(name,"Melange")==0) return MelangeEnum; + else if (strcmp(name,"Water")==0) return WaterEnum; + else if (strcmp(name,"Closed")==0) return ClosedEnum; + else if (strcmp(name,"Free")==0) return FreeEnum; + else if (strcmp(name,"Open")==0) return OpenEnum; + else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; + else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; + else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; + else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; + else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; + else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; + else if (strcmp(name,"Boundary")==0) return BoundaryEnum; + else if (strcmp(name,"Constant")==0) return ConstantEnum; + else if (strcmp(name,"Converged")==0) return ConvergedEnum; + else if (strcmp(name,"ExtToIu")==0) return ExtToIuEnum; + else if (strcmp(name,"Fill")==0) return FillEnum; + else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; + else if (strcmp(name,"Friction")==0) return FrictionEnum; + else if (strcmp(name,"GroundinglineMeltingRate")==0) return GroundinglineMeltingRateEnum; + else if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"IuToExt")==0) return IuToExtEnum; + else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; + else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum; + else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; + else if (strcmp(name,"Misfit")==0) return MisfitEnum; + else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum; + else if (strcmp(name,"Pressure")==0) return PressureEnum; + else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum; + else if (strcmp(name,"QmuPressure")==0) return QmuPressureEnum; + else if (strcmp(name,"QmuVx")==0) return QmuVxEnum; + else if (strcmp(name,"QmuVy")==0) return QmuVyEnum; + else if (strcmp(name,"QmuVz")==0) return QmuVzEnum; + else if (strcmp(name,"QmuThickness")==0) return QmuThicknessEnum; + else if (strcmp(name,"QmuBed")==0) return QmuBedEnum; + else if (strcmp(name,"QmuSurface")==0) return QmuSurfaceEnum; + else if (strcmp(name,"QmuMelting")==0) return QmuMeltingEnum; + else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum; + else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum; + else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; + else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; + else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; + else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; + else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; + else if (strcmp(name,"Temperature")==0) return TemperatureEnum; + else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum; + else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; + else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; + else if (strcmp(name,"Type")==0) return TypeEnum; + else if (strcmp(name,"Vel")==0) return VelEnum; + else if (strcmp(name,"Velocity")==0) return VelocityEnum; + else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + else if (strcmp(name,"Vx")==0) return VxEnum; + else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; + else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else if (strcmp(name,"Vy")==0) return VyEnum; + else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; + else if (strcmp(name,"Vz")==0) return VzEnum; + else if (strcmp(name,"VzMacAyeal")==0) return VzMacAyealEnum; + else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum; + else if (strcmp(name,"VzPicard")==0) return VzPicardEnum; + else if (strcmp(name,"VzStokes")==0) return VzStokesEnum; + else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; + else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; + else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; + else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; + else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; + else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; + else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum; + else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; + else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; + else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; + else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum; + else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; + else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; + else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum; + else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; + else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; + else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; + else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; + else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; + else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; + else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; + else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; + else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; + else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; + else if (strcmp(name,"P0")==0) return P0Enum; + else if (strcmp(name,"P1")==0) return P1Enum; + else if (strcmp(name,"P1DG")==0) return P1DGEnum; + else if (strcmp(name,"BoolElementResult")==0) return BoolElementResultEnum; + else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; + else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum; + else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; + else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; + else if (strcmp(name,"J")==0) return JEnum; + else if (strcmp(name,"Patch")==0) return PatchEnum; + else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum; + else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum; + else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum; + else if (strcmp(name,"PetscVecExternalResult")==0) return PetscVecExternalResultEnum; + else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; + else if (strcmp(name,"Time")==0) return TimeEnum; + else if (strcmp(name,"TriaP1ElementResult")==0) return TriaP1ElementResultEnum; + else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; + else if (strcmp(name,"MinVel")==0) return MinVelEnum; + else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; + else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; + else if (strcmp(name,"MinVy")==0) return MinVyEnum; + else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; + else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; + else if (strcmp(name,"MinVz")==0) return MinVzEnum; + else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; + else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; + else if (strcmp(name,"Relative")==0) return RelativeEnum; + else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; + else if (strcmp(name,"Incremental")==0) return IncrementalEnum; + else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum; + else if (strcmp(name,"None")==0) return NoneEnum; + else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; + else if (strcmp(name,"StokesSolver")==0) return StokesSolverEnum; + else if (strcmp(name,"Adjoint")==0) return AdjointEnum; + else if (strcmp(name,"Colinear")==0) return ColinearEnum; + else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum; + else if (strcmp(name,"Fset")==0) return FsetEnum; + else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; + else if (strcmp(name,"Gradient2")==0) return Gradient2Enum; + else if (strcmp(name,"Gradient3")==0) return Gradient3Enum; + else if (strcmp(name,"Gradient")==0) return GradientEnum; + else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; + else if (strcmp(name,"Gset")==0) return GsetEnum; + else if (strcmp(name,"Index")==0) return IndexEnum; + else if (strcmp(name,"Indexed")==0) return IndexedEnum; + else if (strcmp(name,"Intersect")==0) return IntersectEnum; + else if (strcmp(name,"Nodal")==0) return NodalEnum; + else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; + else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; + else if (strcmp(name,"PetscOptionsAnalyses")==0) return PetscOptionsAnalysesEnum; + else if (strcmp(name,"PetscOptionsStrings")==0) return PetscOptionsStringsEnum; + else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; + else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; + else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum; + else if (strcmp(name,"Regular")==0) return RegularEnum; + else if (strcmp(name,"Scaled")==0) return ScaledEnum; + else if (strcmp(name,"Separate")==0) return SeparateEnum; + else if (strcmp(name,"Sset")==0) return SsetEnum; + else if (strcmp(name,"Verbose")==0) return VerboseEnum; + else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; + else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; + else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; + else if (strcmp(name,"XY")==0) return XYEnum; + else if (strcmp(name,"XYZP")==0) return XYZPEnum; + else if (strcmp(name,"Option")==0) return OptionEnum; + else if (strcmp(name,"OptionCell")==0) return OptionCellEnum; + else if (strcmp(name,"OptionChar")==0) return OptionCharEnum; + else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; + else if (strcmp(name,"OptionDouble")==0) return OptionDoubleEnum; + else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum; + else if (strcmp(name,"Paterson")==0) return PatersonEnum; + else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; + else stage=5; + } + /*If we reach this point, the string provided has not been found*/ + _error_("Enum %s not found",name); } Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.cpp (revision 11417) @@ -769,4 +769,40 @@ } /*}}}*/ +<<<<<<< .working +======= +/*FUNCTION Penta::CreateJacobianMatrix{{{1*/ +void Penta::CreateJacobianMatrix(Mat Jff){ + + /*retrieve parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging {{{2*/ + _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ +#ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + Ke=CreateJacobianDiagnosticHoriz(); + break; +#endif + default: + _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Jff); + delete Ke; + } +} +/*}}}*/ +>>>>>>> .merge-right.r11410 /*FUNCTION Penta::DeepEcho{{{1*/ void Penta::DeepEcho(void){ @@ -1099,6 +1135,8 @@ /*}}}*/ /*FUNCTION Penta::GetStabilizationParameter {{{1*/ -double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){ +double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){ /*Compute stabilization parameter*/ + /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/ + /*kappa=enthalpydiffusionparameter for enthalpy model*/ double normu; @@ -1106,6 +1144,6 @@ normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); - if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){ - tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity)); + if(normu*diameter/(3*2*kappa)<1){ + tau_parameter=pow(diameter,2)/(3*2*2*kappa); } else tau_parameter=diameter/(2*normu); @@ -3283,11 +3321,7 @@ GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&u, gauss); - vy_input->GetInputValue(&v, gauss); - vz_input->GetInputValue(&w, gauss); - vxm_input->GetInputValue(&um,gauss); - vym_input->GetInputValue(&vm,gauss); - vzm_input->GetInputValue(&wm,gauss); - vx=u-um; vy=v-vm; vz=w-wm; + vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; + vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; + vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; D_scalar_advec=gauss->weight*Jdet; @@ -3321,7 +3355,7 @@ vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14; h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.)); - 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); - 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); - 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); + K[0][0]=h/(2*vel)*vx*vx; K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz; + K[1][0]=h/(2*vel)*vy*vx; K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz; + K[2][0]=h/(2*vel)*vz*vx; K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz; D_scalar_stab=gauss->weight*Jdet; if(dt) D_scalar_stab=D_scalar_stab*dt; @@ -3337,5 +3371,5 @@ else if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity); + tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); for(i=0;iGetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); + kappa=thermalconductivity/(rho_ice*heatcapacity); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); @@ -3499,5 +3534,5 @@ GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); - D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity)); + D_scalar_conduct=gauss->weight*Jdet*kappa; if(dt) D_scalar_conduct=D_scalar_conduct*dt; @@ -3568,5 +3603,5 @@ else if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity); + tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); for(i=0;iparameters->FindParam(&dt,TimesteppingTimeStepEnum); this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* enthalpy_input=NULL; - if (dt) enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(inputs); - if (stabilization==2) diameter=MinEdgeLength(xyz_list); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* enthalpy_input=NULL; + Input* enthalpypicard_input=NULL; + if(dt){ + enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + } + if (stabilization==2){ + diameter=MinEdgeLength(xyz_list); + enthalpypicard_input=inputs->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input); + } /* Start looping on the number of gaussian points: */ @@ -3715,5 +3758,5 @@ GetPhi(&phi, &epsilon[0], viscosity); - scalar_def=phi/(rho_ice)*Jdet*gauss->weight; + scalar_def=phi/rho_ice*Jdet*gauss->weight; if(dt) scalar_def=scalar_def*dt; @@ -3733,6 +3776,8 @@ vy_input->GetInputValue(&v, gauss); vz_input->GetInputValue(&w, gauss); - - tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity); + pressure_input->GetInputValue(&pressure, gauss); + enthalpypicard_input->GetInputValue(&enthalpypicard, gauss); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); + tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); @@ -3819,6 +3864,260 @@ double rho_ice,heatcapacity,geothermalflux_value; double basalfriction,alpha2,vx,vy; + double scalar,enthalpy,enthalpyup; + double pressure,pressureup; + double basis[NUMVERTICES]; + Friction* friction=NULL; + GaussPenta* gauss=NULL; + GaussPenta* gaussup=NULL; + + /* Geothermal flux on ice sheet base and basal friction */ + if (!IsOnBed() || IsFloating()) return NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); + + /*Build frictoin element, needed later: */ + friction=new Friction("3d",inputs,matpar,analysis_type); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + gauss=new GaussPenta(0,1,2,2); + gaussup=new GaussPenta(3,4,5,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gaussup->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&basis[0], gauss); + + enthalpy_input->GetInputValue(&enthalpy,gauss); + pressure_input->GetInputValue(&pressure,gauss); +// if(enthalpy>matpar->PureIceEnthalpy(pressure)){ +// enthalpy_input->GetInputValue(&enthalpyup,gaussup); +// pressure_input->GetInputValue(&pressureup,gaussup); +// if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){ +// //do nothing, don't add heatflux +// } +// else{ +// //need to change spcenthalpy according to Aschwanden +// //NEED TO UPDATE +// } +// } +// else{ + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + if(dt) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*basis[i]; +// } + } + + /*Clean up and return*/ + delete gauss; + delete gaussup; + delete friction; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorMelting {{{1*/ +ElementVector* Penta::CreatePVectorMelting(void){ + return NULL; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermal {{{1*/ +ElementVector* Penta::CreatePVectorThermal(void){ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorThermalVolume(); + ElementVector* pe2=CreatePVectorThermalSheet(); + ElementVector* pe3=CreatePVectorThermalShelf(); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/ +ElementVector* Penta::CreatePVectorThermalVolume(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries*/ + int i,j,ig,found=0; + int friction_type,stabilization; + double Jdet,phi,dt; + double rho_ice,heatcapacity; + double thermalconductivity,kappa; + double viscosity,temperature; + double tau_parameter,diameter; + double u,v,w; + double scalar_def,scalar_transient; + double temperature_list[NUMVERTICES]; + double xyz_list[NUMVERTICES][3]; + double L[numdof]; + double dbasis[3][6]; + double epsilon[6]; + GaussPenta *gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + thermalconductivity=matpar->GetThermalConductivity(); + kappa=thermalconductivity/(rho_ice*heatcapacity); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* temperature_input=NULL; + if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs); + if (stabilization==2) diameter=MinEdgeLength(xyz_list); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(2,3); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&L[0], gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + GetPhi(&phi, &epsilon[0], viscosity); + + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; + if(dt) scalar_def=scalar_def*dt; + + for(i=0;ivalues[i]+=scalar_def*L[i]; + + /* Build transient now */ + if(dt){ + temperature_input->GetInputValue(&temperature, gauss); + scalar_transient=temperature*Jdet*gauss->weight; + for(i=0;ivalues[i]+=scalar_transient*L[i]; + } + + if(stabilization==2){ + GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); + + vx_input->GetInputValue(&u, gauss); + vy_input->GetInputValue(&v, gauss); + vz_input->GetInputValue(&w, gauss); + + tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); + + for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); + if(dt){ + for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); + } + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/ +ElementVector* Penta::CreatePVectorThermalShelf(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries */ + int i,j,ig; + double Jdet2d; + double mixed_layer_capacity,thermal_exchange_velocity; + double rho_ice,rho_water,pressure,dt,scalar_ocean; + double heatcapacity,t_pmp; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double basis[NUMVERTICES]; + GaussPenta* gauss=NULL; + + /* Ice/ocean heat exchange flux on ice shelf base */ + if (!IsOnBed() || !IsFloating()) return NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iGetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + gauss=new GaussPenta(0,1,2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&basis[0], gauss); + + pressure_input->GetInputValue(&pressure,gauss); + t_pmp=matpar->TMeltingPoint(pressure); + + scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); + if(dt) scalar_ocean=dt*scalar_ocean; + + for(i=0;ivalues[i]+=scalar_ocean*basis[i]; + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/ +ElementVector* Penta::CreatePVectorThermalSheet(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries */ + int i,j,ig; + int analysis_type; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]={0.0}; + double Jdet2d,dt; + double rho_ice,heatcapacity,geothermalflux_value; + double basalfriction,alpha2,vx,vy; + double basis[NUMVERTICES]; double scalar; - double basis[NUMVERTICES]; Friction* friction=NULL; GaussPenta* gauss=NULL; @@ -3854,245 +4153,14 @@ GetNodalFunctionsP1(&basis[0], gauss); - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); - if(dt) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - delete friction; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorMelting {{{1*/ -ElementVector* Penta::CreatePVectorMelting(void){ - return NULL; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermal {{{1*/ -ElementVector* Penta::CreatePVectorThermal(void){ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorThermalVolume(); - ElementVector* pe2=CreatePVectorThermalSheet(); - ElementVector* pe3=CreatePVectorThermalShelf(); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/ -ElementVector* Penta::CreatePVectorThermalVolume(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries*/ - int i,j,ig,found=0; - int friction_type,stabilization; - double Jdet,phi,dt; - double rho_ice,heatcapacity; - double thermalconductivity; - double viscosity,temperature; - double tau_parameter,diameter; - double u,v,w; - double scalar_def,scalar_transient; - double temperature_list[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double L[numdof]; - double dbasis[3][6]; - double epsilon[6]; - GaussPenta *gauss=NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* temperature_input=NULL; - if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs); - if (stabilization==2) diameter=MinEdgeLength(xyz_list); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,3); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&L[0], gauss); - - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - GetPhi(&phi, &epsilon[0], viscosity); - - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; - if(dt) scalar_def=scalar_def*dt; - - for(i=0;ivalues[i]+=scalar_def*L[i]; - - /* Build transient now */ - if(dt){ - temperature_input->GetInputValue(&temperature, gauss); - scalar_transient=temperature*Jdet*gauss->weight; - for(i=0;ivalues[i]+=scalar_transient*L[i]; - } - - if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - - vx_input->GetInputValue(&u, gauss); - vy_input->GetInputValue(&v, gauss); - vz_input->GetInputValue(&w, gauss); - - tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity); - - for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); - if(dt){ - for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); - } - } - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/ -ElementVector* Penta::CreatePVectorThermalShelf(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j,ig; - double Jdet2d; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,pressure,dt,scalar_ocean; - double heatcapacity,t_pmp; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]; - double basis[NUMVERTICES]; - GaussPenta* gauss=NULL; - - /* Ice/ocean heat exchange flux on ice shelf base */ - if (!IsOnBed() || !IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - for(i=0;iGetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - gauss=new GaussPenta(0,1,2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - pressure_input->GetInputValue(&pressure,gauss); - t_pmp=matpar->TMeltingPoint(pressure); - - scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); - if(dt) scalar_ocean=dt*scalar_ocean; - - for(i=0;ivalues[i]+=scalar_ocean*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/ -ElementVector* Penta::CreatePVectorThermalSheet(void){ - - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - - /*Intermediaries */ - int i,j,ig; - int analysis_type; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[NUMVERTICES2D][3]={0.0}; - double Jdet2d,dt; - double rho_ice,heatcapacity,geothermalflux_value; - double basalfriction,alpha2,vx,vy; - double scalar; - double basis[NUMVERTICES]; - Friction* friction=NULL; - GaussPenta* gauss=NULL; - - /* Geothermal flux on ice sheet base and basal friction */ - if (!IsOnBed() || IsFloating()) return NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - - /*Build frictoin element, needed later: */ - friction=new Friction("3d",inputs,matpar,analysis_type); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - gauss=new GaussPenta(0,1,2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); - if(dt) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); + if(dt) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*basis[i]; } @@ -4271,5 +4339,13 @@ if(converged){ /*Convert enthalpy into temperature and water fraction*/ +<<<<<<< .working for(i=0;iEnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]); +======= + for(i=0;iEnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]); + if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector"); + //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector"); + } +>>>>>>> .merge-right.r11410 this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values)); @@ -7319,4 +7395,188 @@ } /*}}}*/ +/*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/ +ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){ + + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + + switch(approximation){ + case MacAyealApproximationEnum: + return CreateJacobianDiagnosticMacayeal2d(); + case PattynApproximationEnum: + return CreateJacobianDiagnosticPattyn(); + case StokesApproximationEnum: + return CreateJacobianDiagnosticStokes(); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation %s not supported yet",EnumToStringx(approximation)); + } +} +/*}}}*/ +/*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{1*/ +ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){ + + /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build + the stiffness matrix. */ + if (!IsOnBed()) return NULL; + + /*Depth Averaging B*/ + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal(); + delete tria->matice; delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/ +ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double xyz_list[NUMVERTICES][3]; + double Jdet; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double mu_prime; + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double eps1[3],eps2[3]; + double phi[NUMVERTICES]; + double dphi[3][NUMVERTICES]; + GaussPenta *gauss=NULL; + + /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn(); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); + + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; + + for(i=0;i<6;i++){ + for(j=0;j<6;j++){ + eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; + eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; + eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; + eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; + + Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{1*/ +ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){ + + /*Constants*/ + const int numdof=NDOF4*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double xyz_list[NUMVERTICES][3]; + double Jdet; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double eps3dotdphii,eps3dotdphij; + double mu_prime; + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double eps1[3],eps2[3],eps3[3]; + double phi[NUMVERTICES]; + double dphi[3][NUMVERTICES]; + GaussPenta *gauss=NULL; + + /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticStokes(); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); + + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; + + for(i=0;i<6;i++){ + for(j=0;j<6;j++){ + eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; + eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; + eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; + eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; + eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i]; + eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j]; + + Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; + + Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; + + Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Penta.h (revision 11417) @@ -179,5 +179,5 @@ void GetPhi(double* phi, double* epsilon, double viscosity); void GetSolutionFromInputsEnthalpy(Vec solutiong); - double GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity); + double GetStabilizationParameter(double u, double v, double w, double diameter, double kappa); void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input); @@ -231,4 +231,11 @@ ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); ElementMatrix* CreateKMatrixDiagnosticVertSurface(void); +<<<<<<< .working +======= + ElementMatrix* CreateJacobianDiagnosticHoriz(void); + ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void); + ElementMatrix* CreateJacobianDiagnosticPattyn(void); +>>>>>>> .merge-right.r11410 + ElementMatrix* CreateJacobianDiagnosticStokes(void); void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); void InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong); Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11417) @@ -886,4 +886,37 @@ delete gauss; return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreateJacobianMatrix{{{1*/ +void Tria::CreateJacobianMatrix(Mat Jff){ + + /*retrieve parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging {{{2*/ + _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ +#ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + Ke=CreateJacobianDiagnosticMacayeal(); + break; +#endif + default: + _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Jff); + delete Ke; + } } /*}}}*/ @@ -3031,4 +3064,70 @@ } /*}}}*/ +/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/ +ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double xyz_list[NUMVERTICES][3]; + double Jdet,thickness; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double mu_prime; + double epsilon[3];/* epsilon=[exx,eyy,exy];*/ + double eps1[2],eps2[2]; + double phi[NUMVERTICES]; + double dphi[2][NUMVERTICES]; + GaussTria *gauss=NULL; + + /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss); + + thickness_input->GetInputValue(&thickness, gauss); + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; + + for(i=0;i<3;i++){ + for(j=0;j<3;j++){ + eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]; + eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]; + eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]; + eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]; + + Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; + Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; + Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; + Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/ void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ @@ -3332,6 +3431,6 @@ /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ - /*If on water, grad = 0: */ - if(IsOnWater()) return; + /*If on water, grad = 0: (WHY???) */ + if(IsOnWater()) return; /*First deal with ∂/∂alpha(KU-F)*/ Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11417) @@ -84,4 +84,8 @@ void CreateKMatrix(Mat Kff, Mat Kfs,Vec df); void CreatePVector(Vec pf); +<<<<<<< .working +======= + void CreateJacobianMatrix(Mat Jff); +>>>>>>> .merge-right.r11410 int GetNodeIndex(Node* node); int Sid(); @@ -206,4 +210,5 @@ ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorDiagnosticHutter(void); + ElementMatrix* CreateJacobianDiagnosticMacayeal(void); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsDiagnosticHutter(Vec solution); Index: /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.cpp (revision 11417) @@ -24,5 +24,5 @@ KML_File::KML_File(){ - kmlobj =new DataSet; + ; } @@ -31,8 +31,5 @@ KML_File::~KML_File(){ - if (kmlobj) { - delete kmlobj; - kmlobj =NULL; - } + ; } @@ -47,6 +44,4 @@ _printf_(flag,"KML_File:\n"); KML_Object::Echo(); - - _printf_(flag," kmlobj: (size=%d)\n" ,kmlobj->Size()); return; @@ -66,6 +61,4 @@ void KML_File::DeepEcho(const char* indent){ - int i; - char indent2[81]; bool flag=true; @@ -73,4 +66,5 @@ KML_Object::DeepEcho(indent); +<<<<<<< .working /* loop over the kml objects for the file */ @@ -88,4 +82,6 @@ _printf_(flag,"%s kmlobj: [empty]\n" ,indent); +======= +>>>>>>> .merge-right.r11410 return; } @@ -93,7 +89,4 @@ /*FUNCTION KML_File::Write {{{1*/ void KML_File::Write(FILE* filout,const char* indent){ - - int i; - char indent2[81]; fprintf(filout,"%sSize(); i++) - ((KML_Object *)kmlobj->GetObjectByOffset(i))->Write(filout,indent2); fprintf(filout,"%s\n",indent); @@ -145,106 +129,4 @@ _error_("KML_File::Read -- Unexpected field \"%s\".\n",kstri); - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - -// else if (!strncmp(kstri,"Read(fid,kstri); -// kmlobj ->AddObject((Object*)kobj); -// } - -// else if (!strncmp(kstri,"Read(fid,kstri); -// kmlobj ->AddObject((Object*)kobj); -// } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - - else if (!strncmp(kstri,"Read(fid,kstri); - kmlobj ->AddObject((Object*)kobj); - } - -// else if (!strncmp(kstri,"Read(fid,kstri); -// kmlobj ->AddObject((Object*)kobj); -// } - -// else if (!strncmp(kstri,"Read(fid,kstri); -// kmlobj ->AddObject((Object*)kobj); -// } - else if (!strncmp(kstri,"<",1)) KML_Object::Read(fid,kstri); Index: /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/KML/KML_File.h (revision 11417) @@ -19,6 +19,4 @@ public: - - DataSet* kmlobj; /*KML_File constructors, destructors {{{1*/ Index: /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.cpp (revision 11417) @@ -26,4 +26,8 @@ attrib =new DataSet; commnt =new DataSet; +<<<<<<< .working +======= + kmlobj =new DataSet; +>>>>>>> .merge-right.r11410 } @@ -40,4 +44,11 @@ commnt =NULL; } +<<<<<<< .working +======= + if (kmlobj) { + delete kmlobj; + kmlobj =NULL; + } +>>>>>>> .merge-right.r11410 } @@ -52,4 +63,8 @@ _printf_(flag," attrib: (size=%d)\n" ,attrib->Size()); _printf_(flag," commnt: (size=%d)\n" ,commnt->Size()); +<<<<<<< .working +======= + _printf_(flag," kmlobj: (size=%d)\n" ,kmlobj->Size()); +>>>>>>> .merge-right.r11410 return; @@ -70,4 +85,5 @@ int i; + char indent2[81]; bool flag=true; @@ -90,4 +106,21 @@ _printf_(flag,"%s commnt: [empty]\n" ,indent); +<<<<<<< .working +======= +/* loop over the unknown objects for the object */ + + memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char)); + strcat(indent2," "); + + if (kmlobj->Size()) + for (i=0; iSize(); i++) { + _printf_(flag,"%s kmlobj: -------- begin [%d] --------\n" ,indent,i); + ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->DeepEcho(indent2); + _printf_(flag,"%s kmlobj: -------- end [%d] --------\n" ,indent,i); + } + else + _printf_(flag,"%s kmlobj: [empty]\n" ,indent); + +>>>>>>> .merge-right.r11410 return; } @@ -101,4 +134,15 @@ ; +<<<<<<< .working +======= + memcpy(indent2,indent,(strlen(indent)+1)*sizeof(char)); + strcat(indent2," "); + + if (kmlobj->Size()) + for (i=0; iSize(); i++) { + ((KML_Unknown *)kmlobj->GetObjectByOffset(i))->Write(filout,indent2); + } + +>>>>>>> .merge-right.r11410 return; } @@ -116,8 +160,118 @@ _error_("KML_Object::Read -- Unexpected field \"%s\".\n",kstr); + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + +// else if (!strncmp(kstr,"Read(fid,kstr); +// kmlobj ->AddObject((Object*)kobj); +// } + +// else if (!strncmp(kstr,"Read(fid,kstr); +// kmlobj ->AddObject((Object*)kobj); +// } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + + else if (!strncmp(kstr,"Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); + } + +// else if (!strncmp(kstr,"Read(fid,kstr); +// kmlobj ->AddObject((Object*)kobj); +// } + +// else if (!strncmp(kstr,"Read(fid,kstr); +// kmlobj ->AddObject((Object*)kobj); +// } + else if (!strncmp(kstr,"<",1)) { _printf_(true,"KML_Object::Read -- Unrecognized opening tag %s.\n",kstr); +<<<<<<< .working KMLFileTagSkip(kstr, fid); +======= +// KMLFileTagSkip(kstr, +// fid); + kobj=(KML_Object*)new KML_Unknown(); + kobj->Read(fid,kstr); + kmlobj ->AddObject((Object*)kobj); +>>>>>>> .merge-right.r11410 } Index: /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/KML/KML_Object.h (revision 11417) @@ -21,4 +21,8 @@ DataSet* attrib; DataSet* commnt; +<<<<<<< .working +======= + DataSet* kmlobj; +>>>>>>> .merge-right.r11410 /*KML_Object constructors, destructors {{{1*/ Index: /issm/trunk-jpl-damage/src/c/objects/KML/KML_Unknown.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/KML/KML_Unknown.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/KML/KML_Unknown.cpp (revision 11417) @@ -44,9 +44,11 @@ bool flag=true; - _printf_(flag,"KML_Unknown_%s:\n",name); + _printf_(flag,"KML_Unknown %s:\n",name); KML_Object::Echo(); if (value ) _printf_(flag," value: \"%s\"\n" ,value); + else + _printf_(flag," value: [none]\n" ); return; @@ -66,11 +68,27 @@ void KML_Unknown::DeepEcho(const char* indent){ - bool flag=true; - - _printf_(flag,"%sKML_Unknown_%s:\n",indent,name); + char* valuei; + char* vtoken; + char nl[]={'\n','\0'}; + bool flag=true; + + _printf_(flag,"%sKML_Unknown %s:\n",indent,name); KML_Object::DeepEcho(indent); - if (value ) - _printf_(flag,"%s value: \"%s\"\n" ,indent,value); + if (value ) { + valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); + + vtoken=strtok(valuei,nl); + _printf_(flag,"%s value: \"%s" ,indent,vtoken); + + while (vtoken=strtok(NULL,nl)) + _printf_(flag,"\n%s %s" ,indent,vtoken); + _printf_(flag,"\"\n"); + + xfree((void**)&valuei); + } + else + _printf_(flag,"%s value: [none]\n" ,indent); return; @@ -80,4 +98,11 @@ void KML_Unknown::Write(FILE* filout,const char* indent){ +<<<<<<< .working +======= + char* valuei; + char* vtoken; + char nl[]={'\n','\0'}; + +>>>>>>> .merge-right.r11410 fprintf(filout,"%s<%s",indent,name); WriteAttrib(filout," "); @@ -85,9 +110,29 @@ WriteCommnt(filout,indent); +<<<<<<< .working KML_Object::Write(filout,indent); if (value ) fprintf(filout,"%s %s\n",indent,value); - +======= + if (value ) { + valuei=(char *) xmalloc((strlen(value)+1)*sizeof(char)); + memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); + + vtoken=strtok(valuei,nl); + fprintf(filout,"%s %s\n",indent,vtoken); + + while (vtoken=strtok(NULL,nl)) + fprintf(filout,"%s %s\n",indent,vtoken); +>>>>>>> .merge-right.r11410 + +<<<<<<< .working +======= + xfree((void**)&valuei); + } + + KML_Object::Write(filout,indent); + +>>>>>>> .merge-right.r11410 fprintf(filout,"%s\n",indent,name); @@ -101,4 +146,5 @@ int ncom=0; char** pcom=NULL; + char nl[]={'\n','\0'}; /* get object attributes and check for solo tag */ @@ -120,8 +166,22 @@ _error_("KML_Unknown::Read -- Unexpected closing tag %s.\n",kstri); +<<<<<<< .working else if (strncmp(kstri,"<",1)) KMLFileTokenParse( value ,NULL,0, NULL, fid); +======= + else if (strncmp(kstri,"<",1)) { + if (value) { + value=(char *) xrealloc(value,(strlen(value)+1+strlen(kstri)+1)*sizeof(char)); + strcat(value,nl); + strcat(value,kstri); + } + else { + value=(char *) xmalloc((strlen(kstri)+1)*sizeof(char)); + memcpy(value,kstri,(strlen(kstri)+1)*sizeof(char)); + } + } +>>>>>>> .merge-right.r11410 else if (!strncmp(kstri,"<",1)) Index: /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.cpp (revision 11417) @@ -361,4 +361,9 @@ } /*}}}*/ +/*FUNCTION Icefront::CreateJacobianMatrix{{{1*/ +void Icefront::CreateJacobianMatrix(Mat Jff){ + this->CreateKMatrix(Jff,NULL); +} +/*}}}1*/ /*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/ void Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){ @@ -373,4 +378,9 @@ } /*}}}*/ +/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/ +void Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){ + this->PenaltyCreateKMatrix(Jff,NULL,kmax); +} +/*}}}1*/ /*FUNCTION Icefront::InAnalysis{{{1*/ bool Icefront::InAnalysis(int in_analysis_type){ Index: /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Loads/Icefront.h (revision 11417) @@ -76,6 +76,11 @@ void CreateKMatrix(Mat Kff, Mat Kfs); void CreatePVector(Vec pf); +<<<<<<< .working +======= + void CreateJacobianMatrix(Mat Jff); +>>>>>>> .merge-right.r11410 void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); void PenaltyCreatePVector(Vec pf, double kmax); + void PenaltyCreateJacobianMatrix(Mat Jff,double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl-damage/src/c/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Loads/Penpair.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Loads/Penpair.cpp (revision 11417) @@ -214,4 +214,12 @@ } /*}}}1*/ +<<<<<<< .working +======= +/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/ +void Penpair::CreateJacobianMatrix(Mat Jff){ + this->CreateKMatrix(Jff,NULL); +} +/*}}}1*/ +>>>>>>> .merge-right.r11410 /*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/ void Penpair::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){ Index: /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11417) @@ -606,4 +606,35 @@ /*Return: */ *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosity2dDerivativeEpsSquare{{{1*/ +void Matice::GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* epsilon){ + + /*output: */ + double mu_prime; + double mu,n,eff2; + + /*input strain rate: */ + double exx,eyy,exy,exz; + + /*Get visocisty and n*/ + GetViscosity2d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + mu_prime=0.5*pow((double)10,(double)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + exy=epsilon[2]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; } /*}}}*/ Index: /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11417) @@ -69,4 +69,5 @@ void GetViscosityComplement(double* pviscosity_complement, double* pepsilon); void GetViscosityZComplement(double* pviscosity_complement, double* pepsilon); + void GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon); double GetB(); double GetBbar(); Index: /issm/trunk-jpl-damage/src/c/objects/Params/BoolParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/BoolParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/BoolParam.h (revision 11417) @@ -57,5 +57,5 @@ void GetParameterValue(double* pdouble){_error_("Bool param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("Bool param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("Bool param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("Bool param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("Bool param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatArrayParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatArrayParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatArrayParam.h (revision 11417) @@ -60,5 +60,5 @@ void GetParameterValue(double* pdouble){_error_("DoubleMatArray param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("DoubleMatArray param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("DoubleMatArray param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/DoubleMatParam.h (revision 11417) @@ -59,5 +59,5 @@ void GetParameterValue(double* pdouble){_error_("DoubleMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("DoubleMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM,int* pN); Index: /issm/trunk-jpl-damage/src/c/objects/Params/DoubleParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/DoubleParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/DoubleParam.h (revision 11417) @@ -58,5 +58,5 @@ void GetParameterValue(double* pdouble){*pdouble=value;} void GetParameterValue(char** pstring){_error_("Double param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("Double param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("Double param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM); void GetParameterValue(double** pdoublearray,int* pM, int* pN); Index: /issm/trunk-jpl-damage/src/c/objects/Params/DoubleVecParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/DoubleVecParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/DoubleVecParam.h (revision 11417) @@ -55,8 +55,8 @@ void GetParameterValue(int* pinteger){_error_("DoubleVec param of enum %i (%s) cannot return an integer",enum_type,EnumToStringx(enum_type));} void GetParameterValue(int** pintarray,int* pM); - void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("DoubleVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}; + 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));}; void GetParameterValue(double* pdouble){_error_("DoubleVec param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("DoubleVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleVec param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("DoubleVec param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM); void GetParameterValue(double** pdoublearray,int* pM, int* pN); Index: /issm/trunk-jpl-damage/src/c/objects/Params/FileParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/FileParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/FileParam.h (revision 11417) @@ -57,5 +57,5 @@ void GetParameterValue(double* pdouble){_error_("FileParam of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("FileParam of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("FileParam of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("FileParam of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("FileParam of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM, int* pN){_error_("FileParam of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/IntMatParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/IntMatParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/IntMatParam.h (revision 11417) @@ -59,5 +59,5 @@ void GetParameterValue(double* pdouble){_error_("IntMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("IntMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("IntMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));}; Index: /issm/trunk-jpl-damage/src/c/objects/Params/IntParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/IntParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/IntParam.h (revision 11417) @@ -58,5 +58,5 @@ void GetParameterValue(double* pdouble){_error_("Int param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("Int param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("Int param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("Int param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("Int param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/PetscMatParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/PetscMatParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/PetscMatParam.h (revision 11417) @@ -58,5 +58,5 @@ void GetParameterValue(double* pdouble){_error_("PetscMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("PetscMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("PetscMat param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/PetscVecParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/PetscVecParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/PetscVecParam.h (revision 11417) @@ -58,5 +58,5 @@ void GetParameterValue(double* pdouble){_error_("PetscVec param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring){_error_("PetscVec param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} - void GetParameterValue(char*** pstringarray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("PetscVec param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Params/StringParam.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Params/StringParam.h (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Params/StringParam.h (revision 11417) @@ -58,5 +58,5 @@ void GetParameterValue(double* pdouble){_error_("String param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} void GetParameterValue(char** pstring); - void GetParameterValue(char*** pstringarray,int* pM){_error_("String param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToStringx(enum_type));} + void GetParameterValue(char*** pstringarray,int* pM){_error_("String param of enum %i (%s) cannot return a string array",enum_type,EnumToStringx(enum_type));} void GetParameterValue(double** pdoublearray,int* pM){_error_("String param of enum %i (%s) cannot return a double array",enum_type,EnumToStringx(enum_type));} 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));} Index: /issm/trunk-jpl-damage/src/c/objects/Patch.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Patch.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/objects/Patch.cpp (revision 11417) @@ -12,4 +12,5 @@ #include #include +#include #include "./objects.h" #include "../Container/Container.h" Index: /issm/trunk-jpl-damage/src/c/shared/Numerics/Synchronize.sh =================================================================== --- /issm/trunk-jpl-damage/src/c/shared/Numerics/Synchronize.sh (revision 11416) +++ /issm/trunk-jpl-damage/src/c/shared/Numerics/Synchronize.sh (revision 11417) @@ -165,4 +165,5 @@ verbositylevel = (int)level; + return verbositylevel; #else Index: /issm/trunk-jpl-damage/src/c/shared/Numerics/Verbosity.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/shared/Numerics/Verbosity.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/shared/Numerics/Verbosity.cpp (revision 11417) @@ -64,4 +64,5 @@ verbositylevel = (int)level; + return verbositylevel; #else Index: /issm/trunk-jpl-damage/src/c/shared/TriMesh/TriMeshUtils.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/shared/TriMesh/TriMeshUtils.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/shared/TriMesh/TriMeshUtils.cpp (revision 11417) @@ -124,5 +124,5 @@ ******************************************************************************************************************************/ -int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){ +void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){ int i,counter; @@ -379,23 +379,23 @@ IsInPoly ******************************************************************************************************************************/ -int IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){ - - int i; - double x0,y0; - - /*Go through all nodes of the mesh:*/ - for (i=0;i>>>>>> .merge-right.r11410 /*Initialize TAO*/ @@ -38,4 +51,15 @@ if(ierr) _error_("Could not initialize Tao"); +<<<<<<< .working +======= + /*Recover some parameters*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum); + femmodel->parameters->FindParam(&dummy,NULL,NULL,InversionMaxiterPerStepEnum); + maxiter=nsteps*(int)dummy[0]; xfree((void**)&dummy); + +>>>>>>> .merge-right.r11410 /*Initialize TAO*/ TaoCreate(PETSC_COMM_WORLD,&tao); @@ -45,6 +69,6 @@ /*Prepare all TAO parameters*/ - TaoSetMaximumFunctionEvaluations(tao,50); - TaoSetMaximumIterations(tao,10); + TaoSetMaximumFunctionEvaluations(tao,maxiter); + TaoSetMaximumIterations(tao,nsteps); TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28); @@ -66,4 +90,11 @@ InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum); + + /*Finalize*/ + _printf_(VerboseControl(),"%s\n"," preparing final solution"); + femmodel->parameters->SetParam(false,InversionIscontrolEnum); //needed to turn control result output in solutioncore + void (*solutioncore)(FemModel*)=NULL; + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + solutioncore(femmodel); /*Clean up and return*/ @@ -90,4 +121,16 @@ VecFree(&gradient); CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); +<<<<<<< .working +======= + + /*Compute gradient*/ + Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + VecCopy(gradient,G); VecFree(&gradient); + VecScale(G,-1.); + + /*Clean-up and return*/ + xfree((void**)&cost_functions); + xfree((void**)&cost_functionsd); +>>>>>>> .merge-right.r11410 return 0; } Index: /issm/trunk-jpl-damage/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/solutions/transient_core.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/solutions/transient_core.cpp (revision 11417) @@ -24,5 +24,5 @@ /*parameters: */ double finaltime,dt,yts; - bool control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline; + bool control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy; bool dakota_analysis=false; bool time_adapt=false; @@ -51,4 +51,5 @@ femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); @@ -91,5 +92,10 @@ _printf_(VerboseSolution()," computing temperatures:\n"); #ifdef _HAVE_THERMAL_ - thermal_core_step(femmodel,step,time); + if(isenthalpy==0){ + thermal_core_step(femmodel,step,time); + } + else{ + enthalpy_core_step(femmodel,step,time); + } #else _error_("ISSM was not compiled with thermal capabilities. Exiting"); @@ -135,5 +141,7 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum,step,time); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,step,time); + if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum,step,time); Index: /issm/trunk-jpl-damage/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/solvers/solver_nonlinear.cpp (revision 11416) +++ /issm/trunk-jpl-damage/src/c/solvers/solver_nonlinear.cpp (revision 11417) @@ -85,4 +85,7 @@ if(count>=max_nonlinear_iterations){ _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + converged=true; + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); break; } Index: /issm/trunk-jpl-damage/src/m/classes/clusters/castor.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/clusters/castor.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/clusters/castor.m (revision 11417) @@ -10,11 +10,11 @@ % {{{1 name='castor' - login='larour'; - np =128; %number of processors + login='username'; + np =128; port=0; queue='shortc'; time=180; - codepath='/workp/edw/larour/issm-2.0/bin' - executionpath='/workp/edw/larour/Testing/Execution' + codepath='/workp/edw/issm-2.0/bin' + executionpath='/workp/edw/Testing/Execution' %}}} end Index: /issm/trunk-jpl-damage/src/m/classes/clusters/cosmos.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/clusters/cosmos.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/clusters/cosmos.m (revision 11417) @@ -10,11 +10,11 @@ % {{{1 name='cosmos' - login='larour'; + login='username'; np=128; port=0; queue='shortq'; time=3*60; - codepath='/work00/edw/larour/issm-2.0/bin'; - executionpath='/work00/edw/larour/Execution'; + codepath='/work00/edw/issm-2.0/bin'; + executionpath='/work00/edw/Execution'; %}}} end Index: /issm/trunk-jpl-damage/src/m/classes/clusters/gemini.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/clusters/gemini.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/clusters/gemini.m (revision 11417) @@ -10,11 +10,11 @@ % {{{1 name='gemini' - login='larour'; + login='username'; np=50; port=0; queue='debug'; time=60; - codepath='/workg/edw/larour/issm-2.0/bin' - executionpath='/workg/edw/larour/Testing/Execution' + codepath='/workg/edw/issm-2.0/bin' + executionpath='/workg/edw/Testing/Execution' %}}} end Index: /issm/trunk-jpl-damage/src/m/classes/clusters/greenplanet.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/clusters/greenplanet.m (revision 11417) +++ /issm/trunk-jpl-damage/src/m/classes/clusters/greenplanet.m (revision 11417) @@ -0,0 +1,202 @@ +%PFE class definition +% +% Usage: +% cluster=greenplanet(); +% cluster=greenplanet('np',3); +% cluster=greenplanet('np',3,'login','username'); + +classdef greenplanet + properties (SetAccess=public) + % {{{1 + name='greenplanet' + login=''; + numnodes=20; + cpuspernode=8; + port=8000; + queue='rignot'; + codepath=''; + executionpath=''; + interactive=0; + end + properties (SetAccess=private) + np=20*8; + % }}} + end + methods + function cluster=greenplanet(varargin) % {{{1 + + %initialize cluster using default settings if provided + if (exist('greenplanet_settings')==2), greenplanet_settings; end + + %use provided options to change fields + options=pairoptions(varargin{:}); + for i=1:size(options.list,1), + fieldname=options.list{i,1}; + fieldvalue=options.list{i,2}; + if ismember(fieldname,properties('greenplanet')), + cluster.(fieldname)=fieldvalue; + else + disp(['''' fieldname ''' is not a property of cluster greenplanet']); + end + end + end + %}}} + function disp(cluster) % {{{1 + % display the object + disp(sprintf('class ''%s'' object ''%s'' = ',class(cluster),inputname(1))); + disp(sprintf(' name: %s',cluster.name)); + disp(sprintf(' login: %s',cluster.login)); + disp(sprintf(' port: %i',cluster.port)); + disp(sprintf(' numnodes: %i',cluster.numnodes)); + disp(sprintf(' cpuspernode: %i',cluster.cpuspernode)); + disp(sprintf(' np: %i',cluster.cpuspernode*cluster.numnodes)); + disp(sprintf(' queue: %s',cluster.queue)); + disp(sprintf(' codepath: %s',cluster.codepath)); + disp(sprintf(' executionpath: %s',cluster.executionpath)); + disp(sprintf(' interactive: %i',cluster.interactive)); + end + %}}} + function checkconsistency(cluster,md,solution,analyses) % {{{1 + + available_queues={'rignot','default'}; + queue_requirements_time=[Inf Inf]; + queue_requirements_np=[80 80]; + + QueueRequirements(available_queues,queue_requirements_time,queue_requirements_np,cluster.queue,cluster.np,1) + + %Miscelaneous + if isempty(cluster.login), checkmessage('login empty'); end + if isempty(cluster.codepath), checkmessage('codepath empty'); end + if isempty(cluster.executionpath), checkmessage('executionpath empty'); end + + end + %}}} + function BuildQueueScript(cluster,md) % {{{1 + + %retrieve parameters + modelname=md.miscellaneous.name; + solution=md.private.solution; + isvalgrind=md.debug.valgrind; + + %compute number of processors + cluster.np=cluster.numnodes*cluster.cpuspernode; + + %open file for writing: + fid=fopen([modelname '.queue'],'w'); + + fprintf(fid,'#PBS -S /bin/bash\n'); + fprintf(fid,'#PBS -N %s\n',modelname); + fprintf(fid,'#PBS -q %s \n',cluster.queue); + fprintf(fid,'#PBS -l nodes=%i:ppn=%i\n',cluster.numnodes,cluster.cpuspernode); + fprintf(fid,'#PBS -m bea\n'); + fprintf(fid,'#PBS -M mmorligh@uci.edu\n'); + fprintf(fid,'#PBS -o %s.outlog \n',modelname); + fprintf(fid,'#PBS -e %s.errlog \n\n',modelname); + + fprintf(fid,'cd %s/%s\n\n',cluster.executionpath,md.private.runtimename); + fprintf(fid,'/sopt/mpi/mvapich-1.1/intel/bin/mpiexec -mpich-p4-no-shmem -np %i %s/issm.exe %s %s %s\n',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); + + if ~md.settings.io_gather, + %concatenate the output files: + fprintf(fid,'cat %s.outbin.* > %s.outbin',modelname,modelname); + end + + %close file + fclose(fid); + + %in interactive mode, create a run file, and errlog and outlog file + if cluster.interactive, + fid=fopen([modelname '.run'],'w'); + fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s\n',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); + + if ~md.settings.io_gather, + %concatenate the output files: + fprintf(fid,'cat %s.outbin.* > %s.outbin',modelname,modelname); + end + fclose(fid); + fid=fopen([modelname '.errlog'],'w'); + fclose(fid); + fid=fopen([modelname '.outlog'],'w'); + fclose(fid); + end + end %}}} + function LaunchQueueJob(cluster,md,options)% {{{1 + + %lauch command, to be executed via ssh + if ~cluster.interactive, + launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' md.private.runtimename ' && mkdir ' md.private.runtimename ... + ' && cd ' md.private.runtimename ' && mv ../' md.private.runtimename '.tar.gz ./ && tar -zxf ' md.private.runtimename '.tar.gz && qsub ' md.miscellaneous.name '.queue ']; + else + launchcommand=['cd ' cluster.executionpath '/Interactive' num2str(cluster.interactive) ' && tar -zxf ' md.private.runtimename '.tar.gz']; + end + + if ~strcmpi(options.batch,'yes'), + + %compress the files into one zip. + compressstring=['tar -zcf ' md.private.runtimename '.tar.gz ' md.miscellaneous.name '.bin ' md.miscellaneous.name '.queue ' md.miscellaneous.name '.petsc ']; + if md.qmu.isdakota, + compressstring=[compressstring md.miscellaneous.name '.qmu.in ']; + end + if cluster.interactive, + compressstring=[compressstring md.miscellaneous.name '.run ' md.miscellaneous.name '.errlog ' md.miscellaneous.name '.outlog ']; + end + system(compressstring); + + disp('uploading input file and queueing script'); + if cluster.interactive, + directory=[cluster.executionpath '/Interactive' num2str(cluster.interactive)]; + else + directory=cluster.executionpath; + end + + issmscpout(cluster.name,directory,cluster.login,cluster.port,{[md.private.runtimename '.tar.gz']}); + + disp('launching solution sequence on remote cluster'); + issmssh(cluster.name,cluster.login,cluster.port,launchcommand); + + else + disp('batch mode requested: not launching job interactively'); + disp('launch solution sequence on remote cluster by hand'); + end + end + %}}} + function Download(cluster,md)% {{{1 + + %some check + if isempty(md.private.runtimename), + if ~cluster.interactive, + error('greenplanet Download error message: supply runtime name for results to be loaded!'); + end + end + + %Figure out the directory where all the files are in: + if ~cluster.interactive, + directory=[cluster.executionpath '/' md.private.runtimename '/']; + else + directory=[cluster.executionpath '/Interactive' num2str(cluster.interactive) '/']; + end + + %What packages are we picking up from remote cluster + if ~cluster.interactive, + packages={[md.miscellaneous.name '.outlog'],[md.miscellaneous.name '.errlog']}; + else + packages={}; + end + if md.qmu.isdakota, + packages{end+1}=[md.miscellaneous.name '.qmu.err']; + packages{end+1}=[md.miscellaneous.name '.qmu.out']; + if isfield(md.qmu.params,'tabular_graphics_data'), + if md.qmu.params.tabular_graphics_data==true, + packages{end+1}='dakota_tabular.dat'; + end + end + else + packages{end+1}=[md.miscellaneous.name '.outbin']; + end + + %copy files from cluster to present directory + issmscpin(cluster.name, cluster.login, cluster.port, directory, packages); + + end %}}} + end +end Index: /issm/trunk-jpl-damage/src/m/classes/clusters/pollux.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/clusters/pollux.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/clusters/pollux.m (revision 11417) @@ -10,11 +10,11 @@ % {{{1 name='pollux' - login='larour'; + login='username'; np=128; port=0; queue='shortp'; time=180; - codepath='/workc/edw/larour/issm-2.0/bin' - executionpath='/workc/edw/larour/Testing/Execution' + codepath='/workc/edw/issm-2.0/bin' + executionpath='/workc/edw/Testing/Execution' %}}} end Index: /issm/trunk-jpl-damage/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/initialization.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/initialization.m (revision 11417) @@ -62,5 +62,5 @@ checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]); end - if ismember(EnthalpyAnalysisEnum,analyses), + if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy), checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]); end Index: /issm/trunk-jpl-damage/src/m/classes/solver.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/solver.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/solver.m (revision 11417) @@ -6,5 +6,5 @@ classdef solver properties (SetAccess=public) - options={NoneAnalysisEnum,mumpsoptions}; + options=cell(0,0); end methods @@ -26,25 +26,37 @@ function obj = setdefaultparameters(obj) % {{{ + %MUMPS is the default solver + obj.options={'NoneAnalysis',mumpsoptions}; + end % }}} - function obj=addoptions(obj,analysis,solveroptions) % {{{1 + function obj = addoptions(obj,analysis,solveroptions) % {{{1 + + %Convert analysis from enum to string + analysis=EnumToString(analysis); + %first, find out if analysis has already been supplied found=false; for i=1:size(obj.options,1), inanalysis=obj.options{i,1}; - if inanalysis==analysis, + if strcmp(inanalysis,analysis), found=true; - obj.options{i,1}=analysis; - obj.options{i,2}=solveroptions; + obj.options{i,1} = analysis; + obj.options{i,2} = solveroptions; break; end end + if ~found, - obj.options{end+1,1}=analysis; - obj.options{end,2}=solveroptions; + obj.options{end+1,1}= analysis; + obj.options{end,2} = solveroptions; end end %}}} function checkconsistency(obj,md,solution,analyses) % {{{ - + for i=1:size(obj.options,1), + if ~ischar(obj.options{i,1}), + checkmessage('solver is not well formatted: Analyses are not strings'); + end + end end % }}} function PetscFile(solver,filename) % {{{ @@ -70,5 +82,5 @@ %first write analysis: - fprintf(fid,'\n+%s\n',EnumToString(analysis)); %append a + to recognize it's an analysis enum + fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum %now, write options @@ -126,5 +138,5 @@ end - disp(sprintf(' %s -> ''%s''',EnumToString(analysis),string)); + disp(sprintf(' %s -> ''%s''',analysis,string)); end end Index: /issm/trunk-jpl-damage/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl-damage/src/m/classes/thermal.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/classes/thermal.m (revision 11417) @@ -12,4 +12,5 @@ penalty_lock = 0; penalty_factor = 0; + isenthalpy = 0; end methods @@ -42,4 +43,7 @@ %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor obj.penalty_factor=3; + + %Should we use cold ice (default) or enthalpy formulation + obj.isenthalpy=0; end % }}} function checkconsistency(obj,md,solution,analyses) % {{{ @@ -50,6 +54,7 @@ checkfield(md,'thermal.stabilization','numel',1,'values',[0 1 2]); checkfield(md,'thermal.spctemperature','forcing',1); - if ismember(EnthalpyAnalysisEnum,analyses), + if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy), 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'); + checkfield(md,'thermal.isenthalpy','numel',1,'values',[0 1]); end end % }}} @@ -62,4 +67,5 @@ fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)'); + fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); end % }}} @@ -71,4 +77,5 @@ WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean'); end % }}} end Index: /issm/trunk-jpl-damage/src/m/enum/ThermalIsenthalpyEnum.m =================================================================== --- /issm/trunk-jpl-damage/src/m/enum/ThermalIsenthalpyEnum.m (revision 11417) +++ /issm/trunk-jpl-damage/src/m/enum/ThermalIsenthalpyEnum.m (revision 11417) @@ -0,0 +1,11 @@ +function macro=ThermalIsenthalpyEnum() +%THERMALISENTHALPYENUM - Enum of ThermalIsenthalpy +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalIsenthalpyEnum() + +macro=StringToEnum('ThermalIsenthalpy'); Index: /issm/trunk-jpl-damage/src/m/model/collapse.m =================================================================== --- /issm/trunk-jpl-damage/src/m/model/collapse.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/model/collapse.m (revision 11417) @@ -29,4 +29,7 @@ if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; +if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; +if ~isnan(md.inversion.min_parameters), md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; +if ~isnan(md.inversion.max_parameters), md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; if ~isnan(md.surfaceforcings.mass_balance), md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); @@ -49,4 +52,8 @@ if ~isnan(md.flowequation.element_equation) md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); + md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); + md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); + md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); + md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); end @@ -55,4 +62,5 @@ md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); +md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); @@ -60,5 +68,5 @@ %Extrusion of Neumann BC if ~isnan(md.diagnostic.icefront), - numberofneumann2d=size(md.diagnostic.icefront,1)/md.mesh.numberoflayers; + numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer end @@ -89,4 +97,6 @@ md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); +md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); +md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); %Initialize with the 2d mesh Index: /issm/trunk-jpl-damage/src/m/model/loadresultsfromdisk.m =================================================================== --- /issm/trunk-jpl-damage/src/m/model/loadresultsfromdisk.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/model/loadresultsfromdisk.m (revision 11417) @@ -12,4 +12,9 @@ if ~md.qmu.isdakota, + + %Check that file exists + if ~exist(filename,'file'), + error(['binary file ' filename ' not found.']); + end %initialize md.results if not a structure yet Index: /issm/trunk-jpl-damage/src/m/model/radarpower.m =================================================================== --- /issm/trunk-jpl-damage/src/m/model/radarpower.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/model/radarpower.m (revision 11417) @@ -11,5 +11,5 @@ %If gdal does not work, uncomment the following line -setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/'); +%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/'); %Parse inputs if nargin==1, Index: /issm/trunk-jpl-damage/src/m/model/tres.m =================================================================== --- /issm/trunk-jpl-damage/src/m/model/tres.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/model/tres.m (revision 11417) @@ -10,27 +10,15 @@ if strcmpi(string,'diagnostic'), if md.mesh.dimension==2, - if isfield(md.results.DiagnosticSolution,'VxAverage'), - md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.VxAverage); - else - md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.Vx); - end - if isfield(md.results.DiagnosticSolution,'VyAverage'), - md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.VyAverage); - else - md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.Vy); - end + md.initialization.vx=md.results.DiagnosticSolution.Vx; + md.initialization.vy=md.results.DiagnosticSolution.Vy; else - md.initialization.vx=PatchToVec(md.results.DiagnosticSolution.Vx); - md.initialization.vy=PatchToVec(md.results.DiagnosticSolution.Vy); - if isfield(md.results.DiagnosticSolution,'Vz'), - md.initialization.vz=PatchToVec(md.results.DiagnosticSolution.Vz); - else - md.initialization.vz=zeros(md.mesh.numberofvertices,1); - end + md.initialization.vx=md.results.DiagnosticSolution.Vx; + md.initialization.vy=md.results.DiagnosticSolution.Vy; + md.initialization.vz=md.results.DiagnosticSolution.Vz; end - md.initialization.vel=PatchToVec(md.results.DiagnosticSolution.Vel); + md.initialization.vel=md.results.DiagnosticSolution.Vel; if isfield(md.results.DiagnosticSolution,'Pressure'), - md.initialization.pressure=PatchToVec(md.results.DiagnosticSolution.Pressure); + md.initialization.pressure=md.results.DiagnosticSolution.Pressure; end if md.rifts.numrifts, @@ -42,5 +30,5 @@ for control_parameters=md.inversion.control_parameters %Will need to be updated... good luck ;) - md.(EnumToModelField(control_parameters))=PatchToVec(md.results.DiagnosticSolution.(EnumToString(control_parameters))); + md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters)); end end @@ -59,10 +47,10 @@ for i=1:length(results), if ~isempty(md.results.TransientSolution(i).Vel), - results2(count).Vel=PatchToVec(md.results.TransientSolution(i).Vel); - results2(count).Surface=PatchToVec(md.results.TransientSolution(i).Surface); - results2(count).Thickness=PatchToVec(md.results.TransientSolution(i).Thickness); - results2(count).Bed=PatchToVec(md.results.TransientSolution(i).Bed); - results2(count).Vx=PatchToVec(md.results.TransientSolution(i).Vx); - results2(count).Vy=PatchToVec(md.results.TransientSolution(i).Vy); + results2(count).Vel=md.results.TransientSolution(i).Vel; + results2(count).Surface=md.results.TransientSolution(i).Surface; + results2(count).Thickness=md.results.TransientSolution(i).Thickness; + results2(count).Bed=md.results.TransientSolution(i).Bed; + results2(count).Vx=md.results.TransientSolution(i).Vx; + results2(count).Vy=md.results.TransientSolution(i).Vy; results2(count).time=md.results.TransientSolution(i).time; results2(count).step=md.results.TransientSolution(i).step; @@ -76,26 +64,26 @@ clear results,results2; elseif strcmpi(string,'steadystate'), - md.initialization.vx=PatchToVec(md.results.SteadystateSolution.Vx); - md.initialization.vy=PatchToVec(md.results.SteadystateSolution.Vy); + md.initialization.vx=md.results.SteadystateSolution.Vx; + md.initialization.vy=md.results.SteadystateSolution.Vy; if isfield(md.results.SteadystateSolution,'Vz'), - md.initialization.vz=PatchToVec(md.results.SteadystateSolution.Vz); + md.initialization.vz=md.results.SteadystateSolution.Vz; end - md.initialization.vel=PatchToVec(md.results.SteadystateSolution.Vel); - md.initialization.pressure=PatchToVec(md.results.SteadystateSolution.Pressure); - md.initialization.temperature=PatchToVec(md.results.SteadystateSolution.Temperature); - md.basalforcings.melting_rate=PatchToVec(md.results.SteadystateSolution.BasalforcingsMeltingRate); + md.initialization.vel=md.results.SteadystateSolution.Vel; + md.initialization.pressure=md.results.SteadystateSolution.Pressure; + md.initialization.temperature=md.results.SteadystateSolution.Temperature; + md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate; if md.inversion.iscontrol==1, for control_parameters=md.inversion.control_parameters - md.(EnumToModelField(control_parameters))=PatchToVec(md.results.SteadystateSolution.(EnumToString(control_parameters))); + md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters)); end end elseif strcmpi(string,'thermal'), - md.initialization.temperature=PatchToVec(md.results.ThermalSolution.Temperature); - md.basalforcings.melting_rate=PatchToVec(md.results.ThermalSolution.BasalMeltingRate); + md.initialization.temperature=md.results.ThermalSolution.Temperature; + md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate; elseif strcmpi(string,'hydrology'), - md.initialization.watercolumn=PatchToVec(md.results.HydrologySolution.Watercolumn); + md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn; else Index: /issm/trunk-jpl-damage/src/m/qmu/expandresponses.m =================================================================== --- /issm/trunk-jpl-damage/src/m/qmu/expandresponses.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/qmu/expandresponses.m (revision 11417) @@ -1,3 +1,4 @@ function dresp=expandresponses(md,responses) +%EXPANDRESPONSES - expand responses fnames=fieldnames(responses); Index: /issm/trunk-jpl-damage/src/m/solutions/AnalysisConfiguration.m =================================================================== --- /issm/trunk-jpl-damage/src/m/solutions/AnalysisConfiguration.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/solutions/AnalysisConfiguration.m (revision 11417) @@ -42,6 +42,6 @@ case TransientSolutionEnum, - numanalyses=8; - analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum]; + numanalyses=9; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum]; case FlaimSolutionEnum, Index: /issm/trunk-jpl-damage/src/m/solutions/diagnostic_core.m =================================================================== --- /issm/trunk-jpl-damage/src/m/solutions/diagnostic_core.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/solutions/diagnostic_core.m (revision 11417) @@ -15,4 +15,5 @@ ismacayealpattyn=femmodel.parameters.FlowequationIsmacayealpattyn; isstokes=femmodel.parameters.FlowequationIsstokes; + isnewton=femmodel.parameters.DiagnosticIsnewton; dakota_analysis=femmodel.parameters.QmuIsdakota; control_analysis=femmodel.parameters.InversionIscontrol; @@ -53,5 +54,9 @@ issmprintf(VerboseSolution,'\n%s',[' computing horizontal velocities']); femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum); - femmodel=solver_nonlinear(femmodel,modify_loads); + if isnewton, + femmodel=solver_newton(femmodel); + else + femmodel=solver_nonlinear(femmodel,modify_loads); + end end Index: /issm/trunk-jpl-damage/src/m/solutions/transient_core.m =================================================================== --- /issm/trunk-jpl-damage/src/m/solutions/transient_core.m (revision 11416) +++ /issm/trunk-jpl-damage/src/m/solutions/transient_core.m (revision 11417) @@ -19,4 +19,5 @@ isthermal=femmodel.parameters.TransientIsthermal; isgroundingline=femmodel.parameters.TransientIsgroundingline; + isenthalpy=femmodel.parameters.ThermalIsenthalpy; groundinglinemigration=femmodel.parameters.GroundinglineMigration; @@ -57,5 +58,9 @@ if (isthermal & dim==3) issmprintf(VerboseSolution,'\n%s',[' computing temperature']); - femmodel=thermal_core_step(femmodel); + if (isenthalpy==0), + femmodel=thermal_core_step(femmodel); + else + femmodel=enthalpy_core_step(femmodel); + end end @@ -90,4 +95,6 @@ femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time); if (dim==3 & isthermal), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,step,time);end + if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,WaterfractionEnum,step,time);end + if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,EnthalpyEnum,step,time);end femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time); femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceforcingsMassBalanceEnum,step,time); Index: /issm/trunk-jpl-damage/src/m/solvers/solver_newton.m =================================================================== --- /issm/trunk-jpl-damage/src/m/solvers/solver_newton.m (revision 11417) +++ /issm/trunk-jpl-damage/src/m/solvers/solver_newton.m (revision 11417) @@ -0,0 +1,60 @@ +function femmodel=solver_newton(femmodel) +%SOLVER_NEWTON - core solver of diagnostic run +% +% Usage: +% [femmodel]=solver_newton(femmodel) + + %Branch on partitioning schema requested + maxiter=femmodel.parameters.DiagnosticMaxiter; + configuration_type=femmodel.parameters.ConfigurationType; + [femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters); + + %initialize solution vector + converged=0; count=1; + + %Start non-linear iteration using input velocity: + ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters); + uf=Reducevectorgtof( ug, femmodel.nodes,femmodel.parameters); + + %Update the solution to make sure that vx and vxold are similar + [femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum); + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug); + + while(~converged), + + %save pointer to old velocity + old_ug=ug; + old_uf=uf; + + %Solver forward model + [K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); + ys=CreateNodalConstraints(femmodel.nodes,configuration_type); + p_f=Reduceload( p_f, K_fs, ys); + issmprintf(VerboseSolver(),'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters); + ug=Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug); + + %Check convergence + converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters); + if(converged==1) break; end + if(count>maxiter), + issmprintf(true,'%s%i%s',' maximum number of iterations ',maxiter,' exceeded'); + break; + end + + %Prepare next iteration using Newton's method + [K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters); + ys=CreateNodalConstraints(femmodel.nodes,configuration_type); + p_f=Reduceload( p_f, K_fs, ys); + pJf=p_f-K_ff*uf; + Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax); + duf=Solver(Jff,pJf,[],[],femmodel.parameters); + uf=uf+duf; + ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters); + [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug); + + %increase count + count=count+1; + end +end Index: /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp =================================================================== --- /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp (revision 11417) +++ /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.cpp (revision 11417) @@ -0,0 +1,66 @@ +/*\file CreateJacobianMatrix.c + *\brief: build system matrices (stiffness matrix, loads vector) + */ + +#include "./CreateJacobianMatrix.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*input datasets: */ + Elements *elements = NULL; + Nodes *nodes = NULL; + Vertices *vertices = NULL; + Loads *loads = NULL; + Materials *materials = NULL; + Parameters *parameters = NULL; + double kmax; + + /* output datasets: */ + Mat Jff = NULL; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&CreateJacobianMatrixUsage); + + /*Input datasets: */ + FetchMatlabData((DataSet**)&elements,ELEMENTS); + FetchMatlabData((DataSet**)&nodes,NODES); + FetchMatlabData((DataSet**)&vertices,VERTICES); + FetchMatlabData((DataSet**)&loads,LOADS); + FetchMatlabData((DataSet**)&materials,MATERIALS); + FetchMatlabData(¶meters,PARAMETERS); + FetchMatlabData(&kmax,KMAX); + + /*configure: */ + elements-> Configure(elements,loads, nodes,vertices, materials,parameters); + nodes-> Configure(elements,loads, nodes,vertices, materials,parameters); + loads-> Configure(elements, loads, nodes,vertices, materials,parameters); + materials-> Configure(elements, loads, nodes,vertices, materials,parameters); + + /*!Generate internal degree of freedom numbers: */ + CreateJacobianMatrixx(&Jff,elements,nodes,vertices,loads,materials,parameters,kmax); + + /*write output datasets: */ + WriteMatlabData(JFF,Jff); + + /*Free ressources: */ + delete elements; + delete nodes; + delete vertices; + delete loads; + delete materials; + delete parameters; + MatFree(&Jff); + + /*end module: */ + MODULEEND(); +} + +void CreateJacobianMatrixUsage(void) +{ + _printf_(true,"\n"); + _printf_(true," usage: [Jff] = %s(elements,nodes,vertices,loads,materials,parameters,kmax);\n",__FUNCT__); + _printf_(true,"\n"); +} Index: /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.h =================================================================== --- /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.h (revision 11417) +++ /issm/trunk-jpl-damage/src/mex/CreateJacobianMatrix/CreateJacobianMatrix.h (revision 11417) @@ -0,0 +1,37 @@ +/* + CreateJacobianMatrix.h +*/ + +#ifndef _CREATEJACOBIANMATRIX_H +#define _CREATEJACOBIANMATRIX_H + +/* local prototypes: */ +void CreateJacobianMatrixUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/Container/Container.h" +#include "../../c/shared/shared.h" +#include "../../c/EnumDefinitions/EnumDefinitions.h" + +#undef __FUNCT__ +#define __FUNCT__ "CreateJacobianMatrix" + +/* serial input macros: */ +#define ELEMENTS (mxArray *)prhs[0] +#define NODES (mxArray *)prhs[1] +#define VERTICES (mxArray *)prhs[2] +#define LOADS (mxArray *)prhs[3] +#define MATERIALS (mxArray *)prhs[4] +#define PARAMETERS (mxArray *)prhs[5] +#define KMAX (mxArray *)prhs[6] + +/* serial output macros: */ +#define JFF (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 7 + +#endif /* _CREATEJACOBIANMATRIX_H*/ Index: /issm/trunk-jpl-damage/src/mex/Makefile.am =================================================================== --- /issm/trunk-jpl-damage/src/mex/Makefile.am (revision 11416) +++ /issm/trunk-jpl-damage/src/mex/Makefile.am (revision 11417) @@ -22,4 +22,5 @@ CostFunction \ CreateNodalConstraints\ + CreateJacobianMatrix\ Echo\ ElementConnectivity\ @@ -327,4 +328,8 @@ SystemMatrices/SystemMatrices.h + +CreateJacobianMatrix_SOURCES = CreateJacobianMatrix/CreateJacobianMatrix.cpp\ + CreateJacobianMatrix/CreateJacobianMatrix.h + SurfaceArea_SOURCES = SurfaceArea/SurfaceArea.cpp\ SurfaceArea/SurfaceArea.h