Index: /issm/trunk-jpl/src/c/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/CMakeLists.txt (revision 16469) +++ /issm/trunk-jpl/src/c/CMakeLists.txt (revision 16470) @@ -96,5 +96,4 @@ ./classes/Vertex.cpp ./classes/Hook.cpp - ./classes/Patch.cpp ./classes/ElementResults/DoubleElementResult.cpp ./classes/ElementResults/TriaP1ElementResult.cpp @@ -279,5 +278,4 @@ ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp ./analyses/thermal_core.cpp - ./analyses/enthalpy_core.cpp ./solutionsequences/solutionsequence_thermal_nonlinear.cpp) #}}} Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16469) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16470) @@ -63,6 +63,4 @@ ./classes/Hook.h\ ./classes/Hook.cpp\ - ./classes/Patch.h\ - ./classes/Patch.cpp\ ./classes/ElementResults/ElementResultLocal.h\ ./classes/ElementResults/DoubleElementResult.h\ @@ -426,5 +424,4 @@ ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\ ./analyses/thermal_core.cpp\ - ./analyses/enthalpy_core.cpp\ ./solutionsequences/solutionsequence_thermal_nonlinear.cpp #}}} Index: /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 16470) @@ -49,14 +49,9 @@ case ThermalSolutionEnum: - numanalyses=2; + numanalyses=3; analyses=xNew(numanalyses); analyses[0]=ThermalAnalysisEnum; analyses[1]=MeltingAnalysisEnum; - break; - - case EnthalpySolutionEnum: - numanalyses=1; - analyses=xNew(numanalyses); - analyses[0]=EnthalpyAnalysisEnum; + analyses[2]=EnthalpyAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 16470) @@ -40,11 +40,4 @@ #ifdef _HAVE_THERMAL_ solutioncore=&thermal_core; - #else - _error_("ISSM was not compiled with thermal capabilities. Exiting"); - #endif - break; - case EnthalpySolutionEnum: - #ifdef _HAVE_THERMAL_ - solutioncore=&enthalpy_core; #else _error_("ISSM was not compiled with thermal capabilities. Exiting"); Index: /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp (revision 16470) @@ -34,5 +34,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,AdjointEnum); + const char* outputs [] = {"Adjoint"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1); } } Index: /issm/trunk-jpl/src/c/analyses/adjointstressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/adjointstressbalance_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/adjointstressbalance_core.cpp (revision 16470) @@ -35,11 +35,13 @@ /*Save results*/ - if(save_results){ + if(save_results || true){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,AdjointxEnum); - InputToResultx(femmodel,AdjointyEnum); if (isFS){ - InputToResultx(femmodel,AdjointzEnum); - InputToResultx(femmodel,AdjointpEnum); + const char* outputs [] = {"Adjointx","Adjointy","Adjointz","Adjointp"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],4); + } + else{ + const char* outputs [] = {"Adjointx","Adjointy"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2); } } Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 16470) @@ -23,5 +23,4 @@ void hydrology_core(FemModel* femmodel); void thermal_core(FemModel* femmodel); -void enthalpy_core(FemModel* femmodel); void surfaceslope_core(FemModel* femmodel); void bedslope_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 16470) @@ -26,5 +26,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,ThicknessEnum); + const char* outputs [] = {"Thickness"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1); } Index: /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp (revision 16470) @@ -31,7 +31,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,SurfaceSlopeXEnum); - InputToResultx(femmodel,SurfaceSlopeYEnum); - InputToResultx(femmodel,VelEnum); + const char* outputs [] = {"SurfaceSlopeX","SurfaceSlopeY","Vel"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3); } Index: /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 16470) @@ -32,6 +32,12 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,BedSlopeXEnum); - if(meshtype!=Mesh2DverticalEnum) InputToResultx(femmodel,BedSlopeYEnum); + if(meshtype!=Mesh2DverticalEnum){ + const char* outputs [] = {"BedSlopeX"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1); + } + else{ + const char* outputs [] = {"BedSlopeX","BedSlopeY"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2); + } } Index: sm/trunk-jpl/src/c/analyses/enthalpy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 16469) +++ (revision ) @@ -1,35 +1,0 @@ -/*!\file: enthalpy_core.cpp - * \brief: core of the enthalpy solution - */ - -#include "./analyses.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../modules/modules.h" -#include "../shared/shared.h" -#include "../solutionsequences/solutionsequences.h" - -void enthalpy_core(FemModel* femmodel){ - - /*intermediary*/ - bool save_results; - - //first recover parameters common to all solutions - femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - - if(VerboseSolution()) _printf0_(" computing enthalpy\n"); - femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum); - solutionsequence_nonlinear(femmodel,true); - - /*transfer enthalpy to enthalpy picard for the next step: */ - InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum); - - if(save_results){ - if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,TemperatureEnum); - InputToResultx(femmodel,EnthalpyEnum); - InputToResultx(femmodel,WaterfractionEnum); - InputToResultx(femmodel,WatercolumnEnum); - InputToResultx(femmodel,BasalforcingsMeltingRateEnum); - } -} Index: /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 16470) @@ -68,7 +68,6 @@ if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ if(VerboseSolution()) _printf0_(" saving results \n"); - InputToResultx(femmodel,WatercolumnEnum); - InputToResultx(femmodel,HydrologyWaterVxEnum); - InputToResultx(femmodel,HydrologyWaterVyEnum); + const char* outputs [] = {"Watercolumn","HydrologyWaterVx","HydrologyWaterVy"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3); /*unload results*/ @@ -92,6 +91,10 @@ InputToResultx(femmodel,SedimentHeadResidualEnum); if(isefficientlayer){ - InputToResultx(femmodel,EplHeadEnum); - InputToResultx(femmodel,HydrologydcMaskEplactiveEnum); + const char* outputs [] = {"SedimentHead","SedimentHeadResidual","EplHead","HydrologydcMaskEplactive"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],4); + } + else{ + const char* outputs [] = {"SedimentHead","SedimentHeadResidual"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2); } /*unload results*/ Index: /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp (revision 16470) @@ -32,6 +32,6 @@ femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum); - femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum); Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 16470) @@ -33,13 +33,13 @@ IssmDouble reltol; int numoutputs = 0; - char** requested_outputs = NULL; + char** requested_outputs = NULL; /* recover parameters:*/ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum); - femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum); femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); + femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum); @@ -51,13 +51,7 @@ if(VerboseSolution()) _printf0_(" computing temperature and velocity for step: " << step << "\n"); #ifdef _HAVE_THERMAL_ - if(isenthalpy==0){ - thermal_core(femmodel); - femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); - GetSolutionFromInputsx(&tg,femmodel); - } - else{ - enthalpy_core(femmodel); - GetSolutionFromInputsx(&tg,femmodel); - } + thermal_core(femmodel); + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/ + GetSolutionFromInputsx(&tg,femmodel); #else _error_("ISSM was not compiled with thermal capabilities. Exiting"); @@ -85,16 +79,5 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,VxEnum); - InputToResultx(femmodel,VyEnum); - InputToResultx(femmodel,VzEnum); - InputToResultx(femmodel,VelEnum); - InputToResultx(femmodel,PressureEnum); - InputToResultx(femmodel,TemperatureEnum); - if(isenthalpy) InputToResultx(femmodel,WaterfractionEnum); - if(isenthalpy) InputToResultx(femmodel,EnthalpyEnum); - if(isenthalpy) InputToResultx(femmodel,WatercolumnEnum); - //if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum); - InputToResultx(femmodel,BasalforcingsMeltingRateEnum); - femmodel->RequestedOutputsx(requested_outputs,numoutputs); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 16470) @@ -14,7 +14,9 @@ /*parameters: */ bool save_results; + int meshtype; /*Recover some parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); if(VerboseSolution()) _printf0_("computing slope...\n"); @@ -28,6 +30,12 @@ if(save_results){ if(VerboseSolution()) _printf0_("saving results:\n"); - InputToResultx(femmodel,SurfaceSlopeXEnum); - InputToResultx(femmodel,SurfaceSlopeYEnum); + if(meshtype!=Mesh2DverticalEnum){ + const char* outputs [] = {"SurfaceSlopeX"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1); + } + else{ + const char* outputs [] = {"SurfaceSlopeX","SurfaceSlopeY"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2); + } } Index: /issm/trunk-jpl/src/c/analyses/thermal_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 16470) @@ -13,12 +13,16 @@ /*intermediary*/ - bool save_results; - bool dakota_analysis = false; - int solution_type; + bool save_results,isenthalpy; + bool dakota_analysis; + int solution_type,numoutputs; + char** requested_outputs = NULL; - //first recover parameters common to all solutions + /*first recover parameters common to all solutions*/ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + femmodel->parameters->FindParam(&numoutputs,ThermalNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,ThermalRequestedOutputsEnum); if(dakota_analysis && solution_type!=TransientSolutionEnum){ @@ -33,19 +37,32 @@ } - if(VerboseSolution()) _printf0_(" setting basal Dirichlet boundary conditions\n"); - femmodel->UpdateBasalConstraintsEnthalpyx(); + if(isenthalpy){ - if(VerboseSolution()) _printf0_(" computing temperatures\n"); - femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); - solutionsequence_thermal_nonlinear(femmodel); + if(VerboseSolution()) _printf0_(" computing enthalpy\n"); + femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum); + solutionsequence_nonlinear(femmodel,true); - if(VerboseSolution()) _printf0_(" computing melting\n"); - femmodel->SetCurrentConfiguration(MeltingAnalysisEnum); - solutionsequence_linear(femmodel); + /*transfer enthalpy to enthalpy picard for the next step: */ + InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum); + + /*Post process*/ + PostprocessingEnthalpyx(femmodel); + } + else{ + if(VerboseSolution()) _printf0_(" setting basal Dirichlet boundary conditions\n"); + femmodel->UpdateBasalConstraintsEnthalpyx(); + + if(VerboseSolution()) _printf0_(" computing temperatures\n"); + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + solutionsequence_thermal_nonlinear(femmodel); + + if(VerboseSolution()) _printf0_(" computing melting\n"); + femmodel->SetCurrentConfiguration(MeltingAnalysisEnum); + solutionsequence_linear(femmodel); + } if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - InputToResultx(femmodel,TemperatureEnum); - InputToResultx(femmodel,BasalforcingsMeltingRateEnum); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } } Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16470) @@ -22,5 +22,5 @@ int i; IssmDouble starttime,finaltime,dt,yts; - bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isenthalpy,isdelta18o,isgia; + bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isdelta18o,isgia; bool save_results,dakota_analysis; bool time_adapt=false; @@ -49,5 +49,4 @@ femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); - femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); @@ -116,11 +115,5 @@ if(VerboseSolution()) _printf0_(" computing temperatures\n"); #ifdef _HAVE_THERMAL_ - if(isenthalpy==0){ - thermal_core(femmodel); - } - else{ - enthalpy_core(femmodel); - PostprocessingEnthalpyx(femmodel); - } + thermal_core(femmodel); #else _error_("ISSM was not compiled with thermal capabilities. Exiting"); @@ -152,7 +145,6 @@ #endif if(save_results){ - InputToResultx(femmodel,SurfaceEnum); - InputToResultx(femmodel,BedEnum); - InputToResultx(femmodel,MaskGroundediceLevelsetEnum); + const char* outputs [] = {"Surface","Bed","MaskGroundediceLevelset"}; + femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3); } } Index: /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.cpp (revision 16470) @@ -80,15 +80,4 @@ } /*}}}*/ -/*FUNCTION BoolElementResult::PatchFill{{{*/ -void BoolElementResult::PatchFill(int row, Patch* patch){ - - /*Here, we fill the result information into the patch object. First, let's remember what is in a row - * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values - * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */ - IssmDouble IssmDoublevalue=this->value?1:0; - patch->fillresultinfo(row,this->enum_type,this->step,this->time,P0Enum,&IssmDoublevalue,1); - -} -/*}}}*/ /*FUNCTION BoolElementResult::GetVectorFromResults{{{*/ void BoolElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){ Index: /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.h (revision 16470) @@ -37,5 +37,4 @@ int GetStep(void){return step;}; int NumberOfNodalValues(void); - void PatchFill(int row, Patch* patch); /*}}}*/ /*BoolElementResult management: {{{*/ Index: /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.cpp (revision 16470) @@ -80,14 +80,4 @@ } /*}}}*/ -/*FUNCTION DoubleElementResult::PatchFill{{{*/ -void DoubleElementResult::PatchFill(int row, Patch* patch){ - - /*Here, we fill the result information into the patch object. First, let's remember what is in a row - * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values - * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */ - patch->fillresultinfo(row,this->enum_type,this->step,this->time,P0Enum,&this->value,1); - -} -/*}}}*/ /*FUNCTION DoubleElementResult::GetVectorFromResults{{{1*/ void DoubleElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){ Index: /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.h (revision 16470) @@ -37,5 +37,4 @@ int GetStep(void){return step;}; int NumberOfNodalValues(void); - void PatchFill(int row, Patch* patch); /*DoubleElementResult management*/ Index: /issm/trunk-jpl/src/c/classes/ElementResults/ElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/ElementResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/ElementResult.h (revision 16470) @@ -8,5 +8,4 @@ /*Headers:*/ #include "../../datastructures/datastructures.h" -class Patch; class Parameters; @@ -19,5 +18,4 @@ virtual int GetStep(void) = 0; virtual int NumberOfNodalValues(void) = 0; - virtual void PatchFill(int row, Patch *patch)=0; virtual int InstanceEnum() = 0; virtual void GetVectorFromResults(Vector *vector,int*doflist,int*connectivitylist,int numdof)=0; Index: /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.cpp (revision 16470) @@ -83,14 +83,4 @@ } /*}}}*/ -/*FUNCTION PentaP1ElementResult::PatchFill{{{*/ -void PentaP1ElementResult::PatchFill(int row, Patch* patch){ - - /*Here, we fill the result information into the patch object. First, let's remember what is in a row - * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values - * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */ - patch->fillresultinfo(row,this->enum_type,this->step,this->time,P1Enum,this->values,6); - -} -/*}}}*/ /*FUNCTION PentaP1ElementResult::GetVectorFromResults{{{*/ void PentaP1ElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){ Index: /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.h (revision 16470) @@ -36,5 +36,4 @@ int GetStep(void){return step;}; int NumberOfNodalValues(void); - void PatchFill(int row, Patch* patch); /*}}}*/ /*PentaP1ElementResult management: {{{*/ Index: /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp (revision 16470) @@ -82,14 +82,4 @@ } /*}}}*/ -/*FUNCTION TriaP1ElementResult::PatchFill{{{*/ -void TriaP1ElementResult::PatchFill(int row, Patch* patch){ - - /*Here, we fill the result information into the patch object. First, let's remember what is in a row - * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values - * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */ - patch->fillresultinfo(row,this->enum_type,this->step,this->time,P1Enum,this->values,3); - -} -/*}}}*/ /*FUNCTION TriaP1ElementResult::GetVectorFromResults{{{*/ void TriaP1ElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){ Index: /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.h (revision 16470) @@ -35,5 +35,4 @@ int GetStep(void){return step;}; int NumberOfNodalValues(void); - void PatchFill(int row, Patch* patch); /*}}}*/ /*TriaP1ElementResult management: {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16470) @@ -15,5 +15,4 @@ class DataSet; class Parameters; -class Patch; class Elements; class Loads; @@ -60,6 +59,4 @@ virtual void ComputeBasalStress(Vector* sigma_b)=0; virtual void ComputeStrainRate(Vector* eps)=0; - virtual void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes)=0; - virtual void PatchFill(int* pcount, Patch* patch)=0; virtual void ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results)=0; virtual void ResultInterpolation(int* pinterpolation,int output_enum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp (revision 16470) @@ -16,5 +16,4 @@ #include "../ExternalResults/Results.h" #include "../ExternalResults/GenericExternalResult.h" -#include "../Patch.h" #include "../../toolkits/toolkits.h" #include "../../shared/shared.h" @@ -61,78 +60,4 @@ } /*}}}*/ -/*FUNCTION Elements::ResultsToPatch{{{*/ -Patch* Elements::ResultsToPatch(void){ - - /*output: */ - Patch* patch=NULL; - - /*intermediary: */ - int i; - int count; - int numrows; - int numvertices; - int numnodes; - int max_numvertices; - int max_numnodes; - int element_numvertices; - int element_numrows; - int element_numnodes; - - /*We are going to extract from the results within the elements, the desired results , and create a table - * of patch information, that will hold, for each element that computed the result that - * we desire, the enum_type of the result, the step and time, the id of the element, the interpolation type, the vertices ids, and the values - * at the nodes (could be different from the vertices). This will be used for visualization purposes. - * For example, we could build the following patch table, for velocities: - * - 1. on a Beam element, Vx, at step 1, time .5, element id 1, interpolation type P0 (constant), vertices ids 1 and 2, one constant value 4.5 - VxEnum 1 .5 1 P0 1 2 4.5 NaN NaN - 2. on a Tria element, Vz, at step 2, time .8, element id 2, interpolation type P1 (linear), vertices ids 1 3 and 4, with values at 3 nodes 4.5, 3.2, 2.5 - VzEnum 2 .8 2 P1 1 3 4 4.5 3.2 2.5 - * ... etc ... - * - * So what do we need to build the table: the maximum number of vertices included in the table, - * and the maximum number of nodal values, as well as the number of rows. Once we have that, - * we ask the elements to fill their own row in the table, by looping on the elememnts. - * - * We will use the Patch object, which will store all of the information needed, and will be able - * to output itself to disk on its own. See object/Patch.h for format of this object.*/ - - /*First, determine maximum number of vertices, nodes, and number of results: */ - numrows=0; - numvertices=0; - numnodes=0; - - for(i=0;iSize();i++){ - - Element* element=dynamic_cast(this->GetObjectByOffset(i)); - element->PatchSize(&element_numrows,&element_numvertices,&element_numnodes); - - numrows+=element_numrows; - if(element_numvertices>numvertices)numvertices=element_numvertices; - if(element_numnodes>numnodes)numnodes=element_numnodes; - } - - /*Synchronize across cluster, so as to not end up with different sizes for each patch on each cpu: */ - ISSM_MPI_Reduce (&numvertices,&max_numvertices,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&max_numvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - numvertices=max_numvertices; - - ISSM_MPI_Reduce (&numnodes,&max_numnodes,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&max_numnodes,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - numnodes=max_numnodes; - - /*Ok, initialize Patch object: */ - patch=new Patch(numrows,numvertices,numnodes); - - /*Now, go through elements, and fill the Patch object: */ - count=0; - for(i=0;iSize();i++){ - Element* element=dynamic_cast(this->GetObjectByOffset(i)); - element->PatchFill(&count,patch); - } - - return patch; -} -/*}}}*/ /*FUNCTION Elements::SetCurrentConfiguration{{{*/ void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){ @@ -156,5 +81,4 @@ int num_procs; - Patch *patch = NULL; int *resultsenums = NULL; int *resultssizes = NULL; @@ -164,7 +88,6 @@ Vector *vector = NULL; bool io_gather; - bool results_as_patches; - int numberofvertices ,numberofelements; - int numberofresults ,vectorsize; + int numberofvertices,numberofelements; + int numberofresults ,vectorsize; int rank; int minrank; @@ -176,92 +99,66 @@ /*Recover parameters: */ parameters->FindParam(&io_gather,SettingsIoGatherEnum); - parameters->FindParam(&results_as_patches,SettingsResultsAsPatchesEnum); parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum); parameters->FindParam(&numberofelements,MeshNumberofelementsEnum); - if(!results_as_patches){ - /*No patch here, we prepare vectors*/ - - /*Get rank of first cpu that has results*/ - if(this->Size()) rank=my_rank; - else rank=num_procs; - ISSM_MPI_Allreduce (&rank,&minrank,1,ISSM_MPI_INT,ISSM_MPI_MIN,IssmComm::GetComm()); - - /*see what the first element of this partition has in stock (this is common to all partitions)*/ - if(my_rank==minrank){ - if(this->Size()==0) _error_("Cannot write results because there is no element??"); - Element* element=dynamic_cast(this->GetObjectByOffset(0)); - element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults); + /*Get rank of first cpu that has results*/ + if(this->Size()) rank=my_rank; + else rank=num_procs; + ISSM_MPI_Allreduce (&rank,&minrank,1,ISSM_MPI_INT,ISSM_MPI_MIN,IssmComm::GetComm()); + + /*see what the first element of this partition has in stock (this is common to all partitions)*/ + if(my_rank==minrank){ + if(this->Size()==0) _error_("Cannot write results because there is no element??"); + Element* element=dynamic_cast(this->GetObjectByOffset(0)); + element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults); + } + ISSM_MPI_Bcast(&numberofresults,1,ISSM_MPI_INT,minrank,IssmComm::GetComm()); + + /*Get out if there is no results. Otherwise broadcast info*/ + if(!numberofresults) return; + if(my_rank!=minrank){ + resultsenums=xNew(numberofresults); + resultssizes=xNew(numberofresults); + resultstimes=xNew(numberofresults); + resultssteps=xNew(numberofresults); + } + ISSM_MPI_Bcast(resultsenums,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); + ISSM_MPI_Bcast(resultssizes,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); + ISSM_MPI_Bcast(resultstimes,numberofresults,ISSM_MPI_DOUBLE,minrank,IssmComm::GetComm()); + ISSM_MPI_Bcast(resultssteps,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); + + /*Loop over all results and get nodal vector*/ + for(int i=0;i(vectorsize); + + for(int j=0;jSize();j++){ + Element* element=dynamic_cast(this->GetObjectByOffset(j)); + element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]); } - ISSM_MPI_Bcast(&numberofresults,1,ISSM_MPI_INT,minrank,IssmComm::GetComm()); - - /*Get out if there is no results. Otherwise broadcast info*/ - if(!numberofresults) return; - if(my_rank!=minrank){ - resultsenums=xNew(numberofresults); - resultssizes=xNew(numberofresults); - resultstimes=xNew(numberofresults); - resultssteps=xNew(numberofresults); + vector->Assemble(); + + /*Serialize and add to results*/ + vector_serial=vector->ToMPISerial(); + results->DeleteResult(resultsenums[i],resultssteps[i]); + if(my_rank==0){ + /*No need to add this vector for all cpus*/ +#ifdef _HAVE_ADOLC_ + IssmPDouble* vector_serial_passive=xNew(vectorsize); + for(int k=0;k(vector_serial[k]); + results->AddResult(new GenericExternalResult(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i])); + xDelete(vector_serial_passive); +#else + results->AddResult(new GenericExternalResult(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i])); +#endif } - ISSM_MPI_Bcast(resultsenums,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); - ISSM_MPI_Bcast(resultssizes,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); - ISSM_MPI_Bcast(resultstimes,numberofresults,ISSM_MPI_DOUBLE,minrank,IssmComm::GetComm()); - ISSM_MPI_Bcast(resultssteps,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm()); - - /*Loop over all results and get nodal vector*/ - for(int i=0;i(vectorsize); - - for(int j=0;jSize();j++){ - Element* element=dynamic_cast(this->GetObjectByOffset(j)); - element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]); - } - vector->Assemble(); - - /*Serialize and add to results*/ - vector_serial=vector->ToMPISerial(); - results->DeleteResult(resultsenums[i],resultssteps[i]); - if(my_rank==0){ - /*No need to add this vector for all cpus*/ - #ifdef _HAVE_ADOLC_ - IssmPDouble* vector_serial_passive=xNew(vectorsize); - for(int k=0;k(vector_serial[k]); - results->AddResult(new GenericExternalResult(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i])); - xDelete(vector_serial_passive); - #else - results->AddResult(new GenericExternalResult(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i])); - #endif - } - - /*clean up*/ - delete vector; - xDelete(vector_serial); - } - } - else{ - /*create patch object out of all results in this dataset: */ - patch=this->ResultsToPatch(); - - /*Gather onto master cpu 0, if needed: */ - if(io_gather)patch->Gather(); - - /*create result object and add to results dataset:*/ - if (patch->maxvertices && patch->maxnodes){ - results->AddObject(new GenericExternalResult(results->Size()+1,PatchVerticesEnum,patch->maxvertices,1,0)); - results->AddObject(new GenericExternalResult(results->Size()+1,PatchNodesEnum, patch->maxnodes,1,0)); - #ifdef _HAVE_ADOLC_ - IssmPDouble* values_passive=xNew(patch->numrows*patch->numcols); - for(int k=0;k<(patch->numrows*patch->numcols);k++)values_passive[k]=reCast(patch->values[k]); - results->AddObject(new GenericExternalResult(results->Size()+1,PatchEnum, values_passive,patch->numrows,patch->numcols,1,0)); - xDelete(values_passive); - #else - results->AddObject(new GenericExternalResult(results->Size()+1,PatchEnum, patch->values,patch->numrows,patch->numcols,1,0)); - #endif - } + + /*clean up*/ + delete vector; + xDelete(vector_serial); } @@ -271,5 +168,4 @@ xDelete(resultssteps); xDelete(resultstimes); - delete patch; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Elements.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Elements.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Elements.h (revision 16470) @@ -9,5 +9,4 @@ class Loads; class Nodes; -class Patch; class Results; @@ -30,5 +29,4 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); void ToResults(Results* results,Parameters* parameters); - Patch* ResultsToPatch(void); int NumberOfElements(void); void InputDuplicate(int input_enum,int output_enum); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16470) @@ -2882,55 +2882,4 @@ } /*}}}*/ -/*FUNCTION Penta::PatchFill{{{*/ -void Penta::PatchFill(int* pcount, Patch* patch){ - - int i,count; - int vertices_ids[6]; - - /*recover pointer: */ - count=*pcount; - - /*will be needed later: */ - for(i=0;i<6;i++) vertices_ids[i]=vertices[i]->Id(); //vertices id start at column 3 of the patch. - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - - /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand - *it to the result object, to fill the rest: */ - patch->fillelementinfo(count,this->sid+1,vertices_ids,6); - elementresult->PatchFill(count,patch); - - /*increment counter: */ - count++; - } - - /*Assign output pointers:*/ - *pcount=count; -}/*}}}*/ -/*FUNCTION Penta::PatchSize{{{*/ -void Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - int i; - int numrows = 0; - int numnodes = 0; - int temp_numnodes = 0; - - /*Go through all the results objects, and update the counters: */ - for (i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - /*first, we have one more result: */ - numrows++; - /*now, how many vertices and how many nodal values for this result? :*/ - temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object. - if(temp_numnodes>numnodes)numnodes=temp_numnodes; - } - - /*Assign output pointers:*/ - *pnumrows=numrows; - *pnumvertices=NUMVERTICES; - *pnumnodes=numnodes; -} -/*}}}*/ /*FUNCTION Penta::PositiveDegreeDay{{{*/ void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){ @@ -3436,5 +3385,4 @@ bool dakota_analysis; bool isFS; - IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; int numnodes; int* penta_node_ids = NULL; @@ -3444,9 +3392,4 @@ iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); iomodel->Constant(&isFS,FlowequationIsFSEnum); - iomodel->Constant(&beta,MaterialsBetaEnum); - iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum); - iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum); - iomodel->Constant(&meltingpoint,MaterialsMeltingpointEnum); - iomodel->Constant(&latentheat,MaterialsLatentheatEnum); /*Checks if debuging*/ @@ -3636,23 +3579,4 @@ } break; - - - case EnthalpyAnalysisEnum: - /*Initialize mesh velocity*/ - for(i=0;i<6;i++)nodeinputs[i]=0; - if (iomodel->Data(TemperatureEnum) && iomodel->Data(WaterfractionEnum) && iomodel->Data(PressureEnum)) { - for(i=0;i<6;i++){ - if(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1] < meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]){ - nodeinputs[i]=heatcapacity*(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1]-referencetemperature); - } - else nodeinputs[i]=heatcapacity* - (meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]-referencetemperature) - +latentheat*iomodel->Data(WaterfractionEnum)[penta_vertex_ids[i]-1]; - } - this->inputs->AddInput(new PentaInput(EnthalpyEnum,nodeinputs,P1Enum)); - } - else _error_("temperature and waterfraction required for the enthalpy solution"); - break; - default: /*No update for other solution types*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16470) @@ -103,6 +103,4 @@ void ResultInterpolation(int* pinterpolation,int output_enum); void ResultToVector(Vector* vector,int output_enum); - void PatchFill(int* pcount, Patch* patch); - void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); void ResetCoordinateSystem(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16470) @@ -117,6 +117,4 @@ void RequestedOutput(int output_enum,int step,IssmDouble time){_error_("not implemented yet");}; void ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results){_error_("not implemented yet");}; - void PatchFill(int* pcount, Patch* patch){_error_("not implemented yet");}; - void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){_error_("not implemented yet");}; void ResultInterpolation(int* pinterpolation,int output_enum){_error_("not implemented");}; void ResultToVector(Vector* vector,int output_enum){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16470) @@ -2449,55 +2449,4 @@ if(found)*pvalue=value; return found; -} -/*}}}*/ -/*FUNCTION Tria::PatchFill{{{*/ -void Tria::PatchFill(int* prow, Patch* patch){ - - int i,row; - int vertices_ids[3]; - - /*recover pointer: */ - row=*prow; - - for(i=0;i<3;i++) vertices_ids[i]=vertices[i]->Id(); //vertices id start at column 3 of the patch. - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - - /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand - *it to the result object, to fill the rest: */ - patch->fillelementinfo(row,this->sid+1,vertices_ids,3); - elementresult->PatchFill(row,patch); - - /*increment rower: */ - row++; - } - - /*Assign output pointers:*/ - *prow=row; -} -/*}}}*/ -/*FUNCTION Tria::PatchSize{{{*/ -void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - int i; - int numrows = 0; - int numnodes = 0; - int temp_numnodes = 0; - - /*Go through all the results objects, and update the counters: */ - for (i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - /*first, we have one more result: */ - numrows++; - /*now, how many vertices and how many nodal values for this result? :*/ - temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object. - if(temp_numnodes>numnodes)numnodes=temp_numnodes; - } - - /*Assign output pointers:*/ - *pnumrows=numrows; - *pnumvertices=NUMVERTICES; - *pnumnodes=numnodes; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16470) @@ -110,6 +110,4 @@ void ResultInterpolation(int* pinterpolation,int output_enum); void ResultToVector(Vector* vector,int output_enum); - void PatchFill(int* pcount, Patch* patch); - void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); void ResetCoordinateSystem(void); void SmbGradients(); Index: /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 16470) @@ -93,5 +93,5 @@ } /*}}}*/ - GenericExternalResult(int in_id, char* in_enum_string,ResultType in_value,int in_step, IssmDouble in_time){ /*{{{*/ + GenericExternalResult(int in_id,const char* in_enum_string,ResultType in_value,int in_step, IssmDouble in_time){ /*{{{*/ id = in_id; enum_type = -1; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16470) @@ -533,5 +533,5 @@ IssmDouble time; IssmDouble double_result; - char *output_string = NULL; + const char *output_string = NULL; /*recover results*/ Index: sm/trunk-jpl/src/c/classes/Patch.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Patch.cpp (revision 16469) +++ (revision ) @@ -1,160 +1,0 @@ -/*!\file Patch.c - * \brief: implementation of the Patch object - */ - -/*Include files */ -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include -#include -#include -#include "./classes.h" -#include "../shared/shared.h" - -/*Object constructors and destructors:*/ -/*FUNCTION Patch::Patch() default constructor {{{*/ -Patch::Patch(){ - this->numrows = 0; - this->numcols = 0; - this->maxvertices = 0; - this->maxnodes = 0; - this->values = NULL; -} -/*}}}*/ -/*FUNCTION Patch::Patch(int numrows, int maxvertices, int maxnodes){{{*/ -Patch::Patch(int in_numrows, int in_maxvertices, int in_maxnodes){ - - this->numrows=in_numrows; - this->maxvertices=in_maxvertices; - this->maxnodes=in_maxnodes; - this->numcols=1 //enum_type - +1 //step - +1 //time - +1 //element id - +1 //interpolation type - +maxvertices //vertices - +maxnodes; //nodes - - //Allocate values and fill with NaN: - if (this->numcols*this->numrows==0){ - this->values=NULL; - } - else{ - this->values=xNew(this->numcols*this->numrows); - for(int i=0;inumrows;i++){ - for(int j=0;jnumcols;j++){ - this->values[i*this->numcols+j]=NAN; - } - } - } - -} -/*}}}*/ -/*FUNCTION Patch::~Patch(){{{*/ -Patch::~Patch(){ - xDelete(values); -} -/*}}}*/ - -/*Object methods*/ -/*FUNCTION Patch::fillelementinfo{{{*/ -void Patch::fillelementinfo(int count, int element_id, int* vertices_ids, int num_vertices){ - - int i; - IssmDouble* row=NULL; - - /*point to the start of the row: */ - row=this->values+count*this->numcols; - - /*Let's remember what is on a row: - enum_type step time element_id interpolation vertices_ids nodal_values - */ - row[3]=element_id; - for(i=0;ivalues+count*this->numcols; - - /*Let's remember what is on a row: - enum_type step time element_id interpolation vertices_ids nodal_values */ - row[0]=enum_type; - row[1]=(IssmDouble)step; - row[2]=time; - row[4]=interpolation; - for(i=0;imaxvertices+i]=nodal_values[i]; - } - -} -/*}}}*/ -/*FUNCTION Patch::Gather{{{*/ -void Patch::Gather(void){ - - int count; - int my_rank; - int num_procs; - int total_numrows; - int node_numrows; - IssmDouble *total_values = NULL; - ISSM_MPI_Status status; - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - num_procs=IssmComm::GetSize(); - - /*First, figure out total number of rows combining all the cpus: */ - ISSM_MPI_Reduce(&this->numrows,&total_numrows,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_numrows,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*return if patch empty*/ - if(total_numrows==0) return; - - /*Now, allocate buffer to holds all the values, on node 0: */ - if(my_rank==0)total_values=xNew(this->numcols*total_numrows); - - /*Start by copying node 0 values onto total_values: */ - if(my_rank==0){ - count=0; - xMemCpy(total_values+count,this->values,this->numcols*this->numrows); - count+=this->numrows*this->numcols; - } - - /*Now, ask other nodes to send their values: */ - for(int i=1;inumrows,1,ISSM_MPI_INT,0,1,IssmComm::GetComm()); - if (this->numrows)ISSM_MPI_Send(this->values,this->numrows*this->numcols,ISSM_MPI_DOUBLE,0,1,IssmComm::GetComm()); - } - if (my_rank==0){ - ISSM_MPI_Recv(&node_numrows,1,ISSM_MPI_INT,i,1,IssmComm::GetComm(),&status); - if (node_numrows)ISSM_MPI_Recv(total_values+count,node_numrows*this->numcols,ISSM_MPI_DOUBLE,i,1,IssmComm::GetComm(),&status); - count+=node_numrows*this->numcols; - } - } - - /*Now, node 0 has total_values, of size total_numrows*this->numcols. Update the fields in the patch, to reflect this new - * reality. For other cpus, no point in keeping their data anymore: */ - if(my_rank==0){ - this->numrows=total_numrows; - xDelete(this->values); - this->values=total_values; - } - else{ - this->numrows=0; - xDelete(this->values); - } -}/*}}}*/ Index: sm/trunk-jpl/src/c/classes/Patch.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Patch.h (revision 16469) +++ (revision ) @@ -1,44 +1,0 @@ -/*!\file: Patch.h - * \brief prototypes for Patch.h - */ - -/*A Patch object is used to store all results held in elements, in a serial way, and to dump them to disk. - - Each row of the Patch object is made of the following information: - - the result enum_type, - - the step and time, - - the id of the element, - - the interpolation type, - - the vertices ids, - - and the values at the nodes (could be different from the vertices). - - For ex: - 1. on a Beam element, Vx, at step 1, time .5, element id 1, interpolation type P0 (constant), vertices ids 1 and 2, one constant value 4.5 - VxEnum 1 .5 1 P0 1 2 4.5 NaN NaN - 2. on a Tria element, Vz, at step 2, time .8, element id 2, interpolation type P1 (linear), vertices ids 1 3 and 4, with values at 3 nodes 4.5, 3.2, 2.5 - VzEnum 2 .8 2 P1 1 3 4 4.5 3.2 2.5 -*/ - -#ifndef _PATCH_H_ -#define _PATCH_H_ - -/*Headers:*/ -class Patch{ - - public: - - int numrows; //number of results held in Patch object - int numcols; //number of columns - int maxvertices; //largest amount of vertices on a given element, determined by geometry - int maxnodes; //the largest amout of nodes on a given element, determined by interpolation - IssmDouble *values; //result values - - Patch(); - Patch(int numrows, int maxvertices, int maxnodes); - ~Patch(); - void fillelementinfo(int row, int element_id, int* vertices_ids, int num_vertices); - void fillresultinfo(int row,int enum_type,int step, IssmDouble time, int interpolation, IssmDouble* nodal_values, int num_nodes); - void Gather(void); - -}; -#endif //ifndef _PATCH_H_ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 16469) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 16470) @@ -117,5 +117,4 @@ #include "./DofIndexing.h" #include "./IoModel.h" -#include "./Patch.h" #include "./Update.h" #include "./FemModel.h" Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16470) @@ -75,5 +75,4 @@ parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum)); parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum)); - parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsAsPatchesEnum)); parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum)); parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum)); @@ -152,4 +151,9 @@ iomodel->DeleteData(&requestedoutputs,numoutputs,MasstransportRequestedOutputsEnum); + iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum); + parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum); + /*Deal with mass flux segments: {{{*/ iomodel->FetchData(&qmu_mass_flux_present,QmuMassFluxSegmentsPresentEnum); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 16470) @@ -12,6 +12,13 @@ void UpdateElementsEnthalpy(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + bool dakota_analysis; + bool isenthalpy; + /*Now, is the model 3d? otherwise, do nothing: */ if(iomodel->meshtype==Mesh2DhorizontalEnum)return; + + /*Is enthalpy requested?*/ + iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum); + if(!isenthalpy) return; /*Fetch data needed: */ @@ -28,5 +35,4 @@ } - bool dakota_analysis; iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); @@ -47,4 +53,5 @@ iomodel->FetchDataToInput(elements,TemperatureEnum); iomodel->FetchDataToInput(elements,WaterfractionEnum); + iomodel->FetchDataToInput(elements,EnthalpyEnum); iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum); iomodel->FetchDataToInput(elements,WatercolumnEnum); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp (revision 16470) @@ -9,5 +9,5 @@ #include "../ModelProcessorx.h" -void CreateNodesThermal(Nodes** pnodes, IoModel* iomodel){ +void CreateNodesThermal(Nodes** pnodes, IoModel* iomodel){ if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); Index: /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 16470) @@ -8,10 +8,9 @@ #include "../../classes/classes.h" -IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,char* output_string){ +IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,const char* output_string){ - int i; - Definition* definition=NULL; - DataSet* output_definitions=NULL; - IssmDouble return_value; + Definition *definition = NULL; + DataSet *output_definitions = NULL; + IssmDouble return_value; /*Ok, go find the output definitions dataset in the parameters, where our responses are hiding: */ @@ -19,5 +18,5 @@ /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_string: */ - for (i=0;iSize();i++){ + for(int i=0;iSize();i++){ definition=dynamic_cast(output_definitions->GetObjectByOffset(i)); Index: /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.h (revision 16469) +++ /issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.h (revision 16470) @@ -8,5 +8,5 @@ /* local prototypes: */ -IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,char* output_string); +IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,const char* output_string); #endif /* _OUTPUTDEFINITIONSRESPONSEXX_H */ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16469) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16470) @@ -50,8 +50,8 @@ StressbalanceIsnewtonEnum, StressbalanceMaxiterEnum, - StressbalanceNumRequestedOutputsEnum, StressbalancePenaltyFactorEnum, StressbalanceReferentialEnum, StressbalanceReltolEnum, + StressbalanceNumRequestedOutputsEnum, StressbalanceRequestedOutputsEnum, StressbalanceRestolEnum, @@ -228,5 +228,4 @@ SettingsLowmemEnum, SettingsOutputFrequencyEnum, - SettingsResultsAsPatchesEnum, SettingsWaitonlockEnum, SurfaceforcingsDelta18oEnum, @@ -258,5 +257,5 @@ SurfaceforcingsBNegEnum, ThermalIsenthalpyEnum, - ThermalIsdynamicbasalspcEnum, + ThermalIsdynamicbasalspcEnum, ThermalMaxiterEnum, ThermalPenaltyFactorEnum, @@ -265,4 +264,6 @@ ThermalSpctemperatureEnum, ThermalStabilizationEnum, + ThermalNumRequestedOutputsEnum, + ThermalRequestedOutputsEnum, GiaMantleViscosityEnum, GiaLithosphereThicknessEnum, @@ -308,5 +309,4 @@ StressbalanceVerticalAnalysisEnum, EnthalpyAnalysisEnum, - EnthalpySolutionEnum, FlaimAnalysisEnum, FlaimSolutionEnum, @@ -539,7 +539,4 @@ IntExternalResultEnum, JEnum, - PatchEnum, - PatchNodesEnum, - PatchVerticesEnum, PentaP1ElementResultEnum, StringExternalResultEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16470) @@ -58,8 +58,8 @@ case StressbalanceIsnewtonEnum : return "StressbalanceIsnewton"; case StressbalanceMaxiterEnum : return "StressbalanceMaxiter"; - case StressbalanceNumRequestedOutputsEnum : return "StressbalanceNumRequestedOutputs"; case StressbalancePenaltyFactorEnum : return "StressbalancePenaltyFactor"; case StressbalanceReferentialEnum : return "StressbalanceReferential"; case StressbalanceReltolEnum : return "StressbalanceReltol"; + case StressbalanceNumRequestedOutputsEnum : return "StressbalanceNumRequestedOutputs"; case StressbalanceRequestedOutputsEnum : return "StressbalanceRequestedOutputs"; case StressbalanceRestolEnum : return "StressbalanceRestol"; @@ -236,5 +236,4 @@ case SettingsLowmemEnum : return "SettingsLowmem"; case SettingsOutputFrequencyEnum : return "SettingsOutputFrequency"; - case SettingsResultsAsPatchesEnum : return "SettingsResultsAsPatches"; case SettingsWaitonlockEnum : return "SettingsWaitonlock"; case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o"; @@ -273,4 +272,6 @@ case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; case ThermalStabilizationEnum : return "ThermalStabilization"; + case ThermalNumRequestedOutputsEnum : return "ThermalNumRequestedOutputs"; + case ThermalRequestedOutputsEnum : return "ThermalRequestedOutputs"; case GiaMantleViscosityEnum : return "GiaMantleViscosity"; case GiaLithosphereThicknessEnum : return "GiaLithosphereThickness"; @@ -314,5 +315,4 @@ case StressbalanceVerticalAnalysisEnum : return "StressbalanceVerticalAnalysis"; case EnthalpyAnalysisEnum : return "EnthalpyAnalysis"; - case EnthalpySolutionEnum : return "EnthalpySolution"; case FlaimAnalysisEnum : return "FlaimAnalysis"; case FlaimSolutionEnum : return "FlaimSolution"; @@ -529,7 +529,4 @@ case IntExternalResultEnum : return "IntExternalResult"; case JEnum : return "J"; - case PatchEnum : return "Patch"; - case PatchNodesEnum : return "PatchNodes"; - case PatchVerticesEnum : return "PatchVertices"; case PentaP1ElementResultEnum : return "PentaP1ElementResult"; case StringExternalResultEnum : return "StringExternalResult"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16469) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16470) @@ -58,8 +58,8 @@ else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum; else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum; - else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum; else if (strcmp(name,"StressbalanceReferential")==0) return StressbalanceReferentialEnum; else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum; + else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum; else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum; @@ -239,5 +239,4 @@ 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,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum; @@ -260,9 +259,9 @@ else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; + else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; - else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum; + if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum; @@ -279,4 +278,6 @@ else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum; + else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum; + else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum; else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; @@ -320,5 +321,4 @@ else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; 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; @@ -541,7 +541,4 @@ else 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,"StringExternalResult")==0) return StringExternalResultEnum; Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 16469) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 16470) @@ -56,5 +56,5 @@ md = checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]); end - if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy) | solution==EnthalpySolutionEnum(), + if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy) md = checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'initialization.watercolumn' ,'>=',0,'size',[md.mesh.numberofvertices 1]); @@ -106,4 +106,12 @@ WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum); WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum); + + if md.thermal.isenthalpy, + tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; + pos = find(md.initialization.temperature>tpmp); + enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); + enthalpy(pos) = md.materials.heatcapacity*tpmp(pos) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction(pos); + WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum()); + end end % }}} end Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 16469) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 16470) @@ -68,5 +68,5 @@ md = checkfield(md,'initialization.vz','NaN',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices]) - if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy) or solution==EnthalpySolutionEnum(): + if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy): md = checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices]) md = checkfield(md,'initialization.watercolumn' ,'>=',0,'size',[md.mesh.numberofvertices]) @@ -87,3 +87,11 @@ WriteData(fid,'data',self.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum()) WriteData(fid,'data',self.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum()) + + if md.thermal.isenthalpy: + tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; + pos = numpy.nonzero(md.initialization.temperature > tpmp)[0] + enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); + enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,1) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,1) + WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum()); + # }}} Index: /issm/trunk-jpl/src/m/classes/masstransport.m =================================================================== --- /issm/trunk-jpl/src/m/classes/masstransport.m (revision 16469) +++ /issm/trunk-jpl/src/m/classes/masstransport.m (revision 16470) @@ -60,5 +60,5 @@ end end % }}} - function list=defaultoutputs(self,md) % {{{ + function list = defaultoutputs(self,md) % {{{ list = {'Thickness','Surface','Bed'}; Index: /issm/trunk-jpl/src/m/classes/settings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/settings.m (revision 16469) +++ /issm/trunk-jpl/src/m/classes/settings.m (revision 16470) @@ -8,5 +8,4 @@ io_gather = 0; lowmem = 0; - results_as_patches = 0; output_frequency = 0; waitonlock = 0; @@ -38,7 +37,4 @@ obj.output_frequency=1; - %do not use patches by default (difficult to plot) - obj.results_as_patches=0; - %this option can be activated to load automatically the results %onto the model after a parallel run by waiting for the lock file @@ -55,5 +51,4 @@ md = checkfield(md,'settings.io_gather','numel',[1],'values',[0 1]); md = checkfield(md,'settings.lowmem','numel',[1],'values',[0 1]); - md = checkfield(md,'settings.results_as_patches','numel',[1],'values',[0 1]); md = checkfield(md,'settings.output_frequency','numel',[1],'>=',1); md = checkfield(md,'settings.waitonlock','numel',[1]); @@ -65,5 +60,4 @@ fielddisplay(obj,'io_gather','I/O gathering strategy for result outputs (default 1)'); fielddisplay(obj,'lowmem','is the memory limited ? (0 or 1)'); - fielddisplay(obj,'results_as_patches','provide results as patches for each element (0 or 1)'); fielddisplay(obj,'output_frequency','frequency at which results are saved in all solutions with multiple time_steps'); fielddisplay(obj,'waitonlock','maximum number of minutes to wait for batch results (NaN to deactivate)'); @@ -78,5 +72,4 @@ WriteData(fid,'object',obj,'fieldname','io_gather','format','Boolean'); WriteData(fid,'object',obj,'fieldname','lowmem','format','Boolean'); - WriteData(fid,'object',obj,'fieldname','results_as_patches','format','Boolean'); WriteData(fid,'object',obj,'fieldname','output_frequency','format','Integer'); if obj.waitonlock>0, Index: /issm/trunk-jpl/src/m/classes/settings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/settings.py (revision 16469) +++ /issm/trunk-jpl/src/m/classes/settings.py (revision 16470) @@ -15,5 +15,4 @@ self.io_gather = 0 self.lowmem = 0 - self.results_as_patches = 0 self.output_frequency = 0 self.waitonlock = 0 @@ -28,5 +27,4 @@ string="%s\n%s"%(string,fielddisplay(self,"io_gather","I/O gathering strategy for result outputs (default 1)")) string="%s\n%s"%(string,fielddisplay(self,"lowmem","is the memory limited ? (0 or 1)")) - string="%s\n%s"%(string,fielddisplay(self,"results_as_patches","provide results as patches for each element (0 or 1)")) string="%s\n%s"%(string,fielddisplay(self,"output_frequency","frequency at which results are saved in all solutions with multiple time_steps")) string="%s\n%s"%(string,fielddisplay(self,"waitonlock","maximum number of minutes to wait for batch results, or return 0")) @@ -44,7 +42,4 @@ self.output_frequency=1 - #do not use patches by default (difficult to plot) - self.results_as_patches=0 - #this option can be activated to load automatically the results #onto the model after a parallel run by waiting for the lock file @@ -58,5 +53,4 @@ md = checkfield(md,'settings.io_gather','numel',[1],'values',[0,1]) md = checkfield(md,'settings.lowmem','numel',[1],'values',[0,1]) - md = checkfield(md,'settings.results_as_patches','numel',[1],'values',[0,1]) md = checkfield(md,'settings.output_frequency','numel',[1],'>=',1) md = checkfield(md,'settings.waitonlock','numel',[1]) @@ -67,5 +61,4 @@ WriteData(fid,'object',self,'fieldname','io_gather','format','Boolean') WriteData(fid,'object',self,'fieldname','lowmem','format','Boolean') - WriteData(fid,'object',self,'fieldname','results_as_patches','format','Boolean') WriteData(fid,'object',self,'fieldname','output_frequency','format','Integer') if self.waitonlock>0: Index: /issm/trunk-jpl/src/m/classes/steadystate.m =================================================================== --- /issm/trunk-jpl/src/m/classes/steadystate.m (revision 16469) +++ /issm/trunk-jpl/src/m/classes/steadystate.m (revision 16470) @@ -25,4 +25,12 @@ %Relative tolerance for the steadystate convertgence obj.reltol=0.01; + + %default output + obj.requested_outputs={'default'}; + end % }}} + function list=defaultoutputs(self,md) % {{{ + + list = [md.stressbalance.defaultoutputs(md) md.thermal.defaultoutputs(md)]; + end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ @@ -51,5 +59,13 @@ WriteData(fid,'object',obj,'fieldname','reltol','format','Double'); WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer'); - WriteData(fid,'object',obj,'fieldname','requested_outputs','format','StringArray'); + + %process requested outputs + outputs = obj.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(obj,md)]; %add defaults + end + WriteData(fid,'data',outputs,'enum',SteadystateRequestedOutputsEnum,'format','StringArray'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/steadystate.py =================================================================== --- /issm/trunk-jpl/src/m/classes/steadystate.py (revision 16469) +++ /issm/trunk-jpl/src/m/classes/steadystate.py (revision 16470) @@ -29,4 +29,9 @@ return string #}}} + def defaultoutputs(self,md): # {{{ + + return md.stressbalance.defaultoutputs(md)+md.thermal.defaultoutputs(md) + + #}}} def setdefaultparameters(self): # {{{ @@ -37,4 +42,6 @@ self.reltol=0.01 + #default output + self.requested_outputs=['default'] return self #}}} @@ -58,4 +65,12 @@ WriteData(fid,'object',self,'fieldname','reltol','format','Double') WriteData(fid,'object',self,'fieldname','maxiter','format','Integer') - WriteData(fid,'object',self,'fieldname','requested_outputs','format','StringArray') + + #process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] + outputs =outputscopy + print outputs + WriteData(fid,'data',outputs,'enum',SteadystateRequestedOutputsEnum(),'format','StringArray') # }}} Index: /issm/trunk-jpl/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.m (revision 16469) +++ /issm/trunk-jpl/src/m/classes/thermal.m (revision 16470) @@ -14,4 +14,5 @@ isenthalpy = 0; isdynamicbasalspc = 0; + requested_outputs = {}; end methods @@ -23,4 +24,13 @@ error('constructor not supported'); end + end % }}} + function list = defaultoutputs(self,md) % {{{ + + if self.isenthalpy, + list = {'Enthalpy','Temperature','Waterfraction','Watercolumn','BasalforcingsMeltingRate'}; + else + list = {'Temperature','BasalforcingsMeltingRate'}; + end + end % }}} function obj = setdefaultparameters(obj) % {{{ @@ -43,4 +53,7 @@ %will basal boundary conditions be set dynamically obj.isdynamicbasalspc=0; + + %default output + obj.requested_outputs={'default'}; end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ @@ -51,5 +64,5 @@ md = checkfield(md,'thermal.stabilization','numel',[1],'values',[0 1 2]); md = checkfield(md,'thermal.spctemperature','forcing',1); - if (ismember(EnthalpyAnalysisEnum(),analyses) & (md.thermal.isenthalpy | solution==EnthalpySolutionEnum()) & md.mesh.dimension==3), + if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy & md.mesh.dimension==3), pos=find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN); replicate=repmat(md.geometry.surface-md.mesh.z,1,size(md.thermal.spctemperature,2)); @@ -61,4 +74,6 @@ end end + + md = checkfield(md,'thermal.requested_outputs','stringrow',1); end % }}} function disp(obj) % {{{ @@ -73,4 +88,5 @@ fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']); + fielddisplay(obj,'requested_outputs','additional outputs requested'); end % }}} @@ -84,4 +100,13 @@ WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean'); WriteData(fid,'object',obj,'fieldname','isdynamicbasalspc','format','Boolean'); + + %process requested outputs + outputs = obj.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(obj,md)]; %add defaults + end + WriteData(fid,'data',outputs,'enum',ThermalRequestedOutputsEnum(),'format','StringArray'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/thermal.py =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.py (revision 16469) +++ /issm/trunk-jpl/src/m/classes/thermal.py (revision 16470) @@ -22,4 +22,5 @@ self.isenthalpy = 0 self.isdynamicbasalspc = 0; + self.requested_outputs = [] #set defaults @@ -36,6 +37,15 @@ string="%s\n%s"%(string,fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)')) + string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) return string #}}} + def defaultoutputs(self,md): # {{{ + + if self.isenthalpy: + return ['Enthalpy','Temperature','Waterfraction','Watercolumn','BasalforcingsMeltingRate'] + else: + return ['Temperature','BasalforcingsMeltingRate'] + + #}}} def setdefaultparameters(self): # {{{ @@ -57,4 +67,7 @@ #will basal boundary conditions be set dynamically self.isdynamicbasalspc=0; + + #default output + self.requested_outputs=['default'] return self @@ -68,5 +81,5 @@ md = checkfield(md,'thermal.stabilization','numel',[1],'values',[0,1,2]) md = checkfield(md,'thermal.spctemperature','forcing',1) - if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3: + if EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy and md.mesh.dimension==3: pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))) replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) @@ -74,4 +87,5 @@ md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1]) md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); + md = checkfield(md,'thermal.requested_outputs','stringrow',1) return md @@ -86,3 +100,11 @@ WriteData(fid,'object',self,'fieldname','isenthalpy','format','Boolean') WriteData(fid,'object',self,'fieldname','isdynamicbasalspc','format','Boolean'); + + #process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] + outputs =outputscopy + WriteData(fid,'data',outputs,'enum',ThermalRequestedOutputsEnum(),'format','StringArray') # }}} Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 16469) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 16470) @@ -55,9 +55,5 @@ case ThermalSolutionEnum(), numanalyses=2; - analyses=[ThermalAnalysisEnum();MeltingAnalysisEnum()]; - - case EnthalpySolutionEnum(), - numanalyses=1; - analyses=[EnthalpyAnalysisEnum()]; + analyses=[EnthalpyAnalysisEnum;ThermalAnalysisEnum();MeltingAnalysisEnum()]; case MasstransportSolutionEnum(), Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 16469) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 16470) @@ -20,9 +20,5 @@ elif solutiontype == ThermalSolutionEnum(): numanalyses=2 - analyses=[ThermalAnalysisEnum(),MeltingAnalysisEnum()] - - elif solutiontype == EnthalpySolutionEnum(): - numanalyses=1 - analyses=[EnthalpyAnalysisEnum()] + analyses=[EnthalpyAnalysisEnum(),ThermalAnalysisEnum(),MeltingAnalysisEnum()] elif solutiontype == MasstransportSolutionEnum(): Index: sm/trunk-jpl/src/m/enum/EnthalpySolutionEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/EnthalpySolutionEnum.m (revision 16469) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=EnthalpySolutionEnum() -%ENTHALPYSOLUTIONENUM - Enum of EnthalpySolution -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=EnthalpySolutionEnum() - -macro=StringToEnum('EnthalpySolution'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16469) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16470) @@ -50,8 +50,8 @@ def StressbalanceIsnewtonEnum(): return StringToEnum("StressbalanceIsnewton")[0] def StressbalanceMaxiterEnum(): return StringToEnum("StressbalanceMaxiter")[0] -def StressbalanceNumRequestedOutputsEnum(): return StringToEnum("StressbalanceNumRequestedOutputs")[0] def StressbalancePenaltyFactorEnum(): return StringToEnum("StressbalancePenaltyFactor")[0] def StressbalanceReferentialEnum(): return StringToEnum("StressbalanceReferential")[0] def StressbalanceReltolEnum(): return StringToEnum("StressbalanceReltol")[0] +def StressbalanceNumRequestedOutputsEnum(): return StringToEnum("StressbalanceNumRequestedOutputs")[0] def StressbalanceRequestedOutputsEnum(): return StringToEnum("StressbalanceRequestedOutputs")[0] def StressbalanceRestolEnum(): return StringToEnum("StressbalanceRestol")[0] @@ -228,5 +228,4 @@ def SettingsLowmemEnum(): return StringToEnum("SettingsLowmem")[0] def SettingsOutputFrequencyEnum(): return StringToEnum("SettingsOutputFrequency")[0] -def SettingsResultsAsPatchesEnum(): return StringToEnum("SettingsResultsAsPatches")[0] def SettingsWaitonlockEnum(): return StringToEnum("SettingsWaitonlock")[0] def SurfaceforcingsDelta18oEnum(): return StringToEnum("SurfaceforcingsDelta18o")[0] @@ -265,4 +264,6 @@ def ThermalSpctemperatureEnum(): return StringToEnum("ThermalSpctemperature")[0] def ThermalStabilizationEnum(): return StringToEnum("ThermalStabilization")[0] +def ThermalNumRequestedOutputsEnum(): return StringToEnum("ThermalNumRequestedOutputs")[0] +def ThermalRequestedOutputsEnum(): return StringToEnum("ThermalRequestedOutputs")[0] def GiaMantleViscosityEnum(): return StringToEnum("GiaMantleViscosity")[0] def GiaLithosphereThicknessEnum(): return StringToEnum("GiaLithosphereThickness")[0] @@ -306,5 +307,4 @@ def StressbalanceVerticalAnalysisEnum(): return StringToEnum("StressbalanceVerticalAnalysis")[0] def EnthalpyAnalysisEnum(): return StringToEnum("EnthalpyAnalysis")[0] -def EnthalpySolutionEnum(): return StringToEnum("EnthalpySolution")[0] def FlaimAnalysisEnum(): return StringToEnum("FlaimAnalysis")[0] def FlaimSolutionEnum(): return StringToEnum("FlaimSolution")[0] @@ -521,7 +521,4 @@ def IntExternalResultEnum(): return StringToEnum("IntExternalResult")[0] def JEnum(): return StringToEnum("J")[0] -def PatchEnum(): return StringToEnum("Patch")[0] -def PatchNodesEnum(): return StringToEnum("PatchNodes")[0] -def PatchVerticesEnum(): return StringToEnum("PatchVertices")[0] def PentaP1ElementResultEnum(): return StringToEnum("PentaP1ElementResult")[0] def StringExternalResultEnum(): return StringToEnum("StringExternalResult")[0] Index: sm/trunk-jpl/src/m/enum/PatchEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/PatchEnum.m (revision 16469) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=PatchEnum() -%PATCHENUM - Enum of Patch -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=PatchEnum() - -macro=StringToEnum('Patch'); Index: sm/trunk-jpl/src/m/enum/PatchNodesEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/PatchNodesEnum.m (revision 16469) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=PatchNodesEnum() -%PATCHNODESENUM - Enum of PatchNodes -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=PatchNodesEnum() - -macro=StringToEnum('PatchNodes'); Index: sm/trunk-jpl/src/m/enum/PatchVerticesEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/PatchVerticesEnum.m (revision 16469) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=PatchVerticesEnum() -%PATCHVERTICESENUM - Enum of PatchVertices -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=PatchVerticesEnum() - -macro=StringToEnum('PatchVertices'); Index: sm/trunk-jpl/src/m/enum/SettingsResultsAsPatchesEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SettingsResultsAsPatchesEnum.m (revision 16469) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SettingsResultsAsPatchesEnum() -%SETTINGSRESULTSASPATCHESENUM - Enum of SettingsResultsAsPatches -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=SettingsResultsAsPatchesEnum() - -macro=StringToEnum('SettingsResultsAsPatches'); Index: /issm/trunk-jpl/src/m/enum/ThermalNumRequestedOutputsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/ThermalNumRequestedOutputsEnum.m (revision 16470) +++ /issm/trunk-jpl/src/m/enum/ThermalNumRequestedOutputsEnum.m (revision 16470) @@ -0,0 +1,11 @@ +function macro=ThermalNumRequestedOutputsEnum() +%THERMALNUMREQUESTEDOUTPUTSENUM - Enum of ThermalNumRequestedOutputs +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=ThermalNumRequestedOutputsEnum() + +macro=StringToEnum('ThermalNumRequestedOutputs'); Index: /issm/trunk-jpl/src/m/enum/ThermalRequestedOutputsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/ThermalRequestedOutputsEnum.m (revision 16470) +++ /issm/trunk-jpl/src/m/enum/ThermalRequestedOutputsEnum.m (revision 16470) @@ -0,0 +1,11 @@ +function macro=ThermalRequestedOutputsEnum() +%THERMALREQUESTEDOUTPUTSENUM - Enum of ThermalRequestedOutputs +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=ThermalRequestedOutputsEnum() + +macro=StringToEnum('ThermalRequestedOutputs'); Index: /issm/trunk-jpl/src/m/geometry/NowickiProfile.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/NowickiProfile.m (revision 16469) +++ /issm/trunk-jpl/src/m/geometry/NowickiProfile.m (revision 16470) @@ -12,5 +12,4 @@ delta = 0.1; % ratio of water density and ice density -1 hg = 1; % ice thickness at grounding line -sea = hg/(1+delta); % sea level lamda = 0.1; % ration of deviatoric stress and water pressure beta = 5; % friction coefficient @@ -36,4 +35,4 @@ for i = ceil(numel(x)/2):numel(x) h(i) = (x(i)/(4*(delta+1)*q)+hg^(-2))^(-0.5); % ice thickness for ice shelf from (3.1) - b(i) = sea-h(i)*(1/(1+delta)); + b(i) = -h(i)*(1/(1+delta)); end Index: sm/trunk-jpl/src/m/interp/PatchToVec.m =================================================================== --- /issm/trunk-jpl/src/m/interp/PatchToVec.m (revision 16469) +++ (revision ) @@ -1,21 +1,0 @@ -function vec=PatchToVec(Patch) -%PATCHTOVEC - converts a patch to a vector by averaging over each vertex -% -% Usage: -% vec=PatchToVec(Patch) - -%if the patch is P0: we have element values, return an element vector -switch(size(Patch.value,2)), - case 1, - vec(Patch.element)=Patch.value; - case 3, - connectivity=sparse(Patch.index(:),1,1); - value =sparse(Patch.index(:),1,Patch.value(:)); - vec=full(value./connectivity); - case 6, - connectivity=sparse(Patch.index(:),1,1); - value =sparse(Patch.index(:),1,Patch.value(:)); - vec=full(value./connectivity); - otherwise, - error('interpolation not supported yet'); - end Index: /issm/trunk-jpl/src/m/plot/plot_unit.m =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.m (revision 16469) +++ /issm/trunk-jpl/src/m/plot/plot_unit.m (revision 16470) @@ -57,33 +57,4 @@ end - %Patch plot P1 - case 4, - - if is2d, - patch( 'Faces',elements,'Vertices',[x y],'CData',data,'FaceColor','interp','EdgeColor',edgecolor); - else - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor); - end - - %Patch plot P0 - case 5, - - if is2d, - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', data(:),'FaceColor','flat','EdgeColor',edgecolor); - else - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor); - end - otherwise, error(['case ' num2str(datatype) ' not supported']); Index: /issm/trunk-jpl/src/m/plot/processdata.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processdata.m (revision 16469) +++ /issm/trunk-jpl/src/m/plot/processdata.m (revision 16470) @@ -37,95 +37,72 @@ end -%Process Patch -if isstruct(data) - if (isprop(data,'index') & isprop(data,'value')), - if data.interpolation(1)==P1Enum(), - data=data.value; - data=data'; - data=data(:); - datatype=4; - elseif data.interpolation(1)==P0Enum(), - data=data.value; - datatype=5; - else - error(['interpolation ' data.interpolation(1) ' not supported yet']); - end - else - error('structure other than Patch not supported yet'); - end -else - %initialize datatype - datatype=0; -end +%initialize datatype +datatype=0; %get datasize datasize=size(data); -%non patch processing -if datatype~=4 & datatype~=5, +%transpose data if necessary +if (size(data,2) > size(data,1)), + data=data'; +end +datasize=size(data); - %transpose data if necessary - if (size(data,2) > size(data,1)), - data=data'; - end - datasize=size(data); +%convert to double if necessary +if ~isnumeric(data); + disp('processdata info message: data is not numeric (logical?). Converted to double'); + data=double(data); +end - %convert to double if necessary - if ~isnumeric(data); - disp('processdata info message: data is not numeric (logical?). Converted to double'); - data=double(data); - end +%check length +if datasize(1)~=md.mesh.numberofvertices & datasize(1)~=md.mesh.numberofelements & datasize(1)~=md.mesh.numberofvertices*6 & (md.mesh.dimension==3 & ~(datasize(1)==numberofelements2d | datasize(1)==numberofvertices2d)) + error('plotmodel error message: data not supported yet'); +end - %check length - if datasize(1)~=md.mesh.numberofvertices & datasize(1)~=md.mesh.numberofelements & datasize(1)~=md.mesh.numberofvertices*6 & (md.mesh.dimension==3 & ~(datasize(1)==numberofelements2d | datasize(1)==numberofvertices2d)) - error('plotmodel error message: data not supported yet'); - end +%quiver? +if datasize(2)>1, + datatype=3; - %quiver? - if datasize(2)>1, - datatype=3; - - %check number of columns, add zeros if necessary, - if (md.mesh.dimension==3) - if datasize(2)==2, - data=[data, zeros(datasize(1),1)]; - elseif datasize(2)~=3, - error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot'); - end + %check number of columns, add zeros if necessary, + if (md.mesh.dimension==3) + if datasize(2)==2, + data=[data, zeros(datasize(1),1)]; + elseif datasize(2)~=3, + error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot'); + end %elseif ((md.mesh.dimension==2) & datasize(2)~=2), % error('plotmodel error message: data provided should have 2 columns for quiver plot, and 1 for regular plot'); - end end +end - %treat the case datasize(1)=6*nodes - if datasize(1)==6*md.mesh.numberofvertices - %keep the only norm of data - data1=data(1:6:md.mesh.numberofvertices*6,:); - data2=data(2:6:md.mesh.numberofvertices*6,:); - data=sqrt(data1.^2+data2.^2); - datasize(1)=md.mesh.numberofvertices; - %---> go to node data - end +%treat the case datasize(1)=6*nodes +if datasize(1)==6*md.mesh.numberofvertices + %keep the only norm of data + data1=data(1:6:md.mesh.numberofvertices*6,:); + data2=data(2:6:md.mesh.numberofvertices*6,:); + data=sqrt(data1.^2+data2.^2); + datasize(1)=md.mesh.numberofvertices; + %---> go to node data +end - %treat the case datasize(1)=nodes2d - if (md.mesh.dimension==3 & datasize(1)==numberofvertices2d), - data=project3d(md,'vector',data,'type','node'); - datasize(1)=md.mesh.numberofvertices; - %---> go to node data - end +%treat the case datasize(1)=nodes2d +if (md.mesh.dimension==3 & datasize(1)==numberofvertices2d), + data=project3d(md,'vector',data,'type','node'); + datasize(1)=md.mesh.numberofvertices; + %---> go to node data +end - %treat the case datasize(1)=nodes2d - if (md.mesh.dimension==3 & datasize(1)==numberofelements2d), - data=project3d(md,'vector',data,'type','element'); - datasize(1)=md.mesh.numberofelements; - %---> go to node data - end +%treat the case datasize(1)=nodes2d +if (md.mesh.dimension==3 & datasize(1)==numberofelements2d), + data=project3d(md,'vector',data,'type','element'); + datasize(1)=md.mesh.numberofelements; + %---> go to node data +end - %smoothing? - if exist(options,'smooth') - data=averaging(md,data,getfieldvalue(options,'smooth')); - datasize(1)=md.mesh.numberofvertices; - %---> go to node data - end +%smoothing? +if exist(options,'smooth') + data=averaging(md,data,getfieldvalue(options,'smooth')); + datasize(1)=md.mesh.numberofvertices; + %---> go to node data end Index: /issm/trunk-jpl/src/m/plot/processmesh.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processmesh.m (revision 16469) +++ /issm/trunk-jpl/src/m/plot/processmesh.m (revision 16470) @@ -21,5 +21,5 @@ end -if (isempty(data) | ~isstruct(data)), +if isempty(data) %first load x,y, etc ... to speed up plot @@ -70,19 +70,4 @@ elements=elements2d; end -else - %Process Patch - if (md.mesh.dimension==2), - elements=transpose(reshape(1:3*md.mesh.numberofelements,3,md.mesh.numberofelements)); - x=transpose(reshape(md.mesh.x(data.index)',1,3*md.mesh.numberofelements)); - y=transpose(reshape(md.mesh.y(data.index)',1,3*md.mesh.numberofelements)); - z=zeros(3*md.mesh.numberofelements,1); - is2d=1; - else - elements=transpose(reshape(1:6*md.mesh.numberofelements,6,md.mesh.numberofelements)); - x=transpose(reshape(md.mesh.x(data.index)',1,6*md.mesh.numberofelements)); - y=transpose(reshape(md.mesh.y(data.index)',1,6*md.mesh.numberofelements)); - z=transpose(reshape(md.mesh.z(data.index)',1,6*md.mesh.numberofelements)); - is2d=0; - end end Index: /issm/trunk-jpl/src/m/regional/regionaltransient2d.m =================================================================== --- /issm/trunk-jpl/src/m/regional/regionaltransient2d.m (revision 16469) +++ /issm/trunk-jpl/src/m/regional/regionaltransient2d.m (revision 16470) @@ -102,7 +102,7 @@ for t=find(numElements==1) if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0, - vx=PatchToVec(md1.results.TransientSolution(t).Vx); - vy=PatchToVec(md1.results.TransientSolution(t).Vy); - thickness=PatchToVec(md1.results.TransientSolution(t).Thickness); + vx=md1.results.TransientSolution(t).Vx; + vy=md1.results.TransientSolution(t).Vy; + thickness=md1.results.TransientSolution(t).Thickness; spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)]; spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vy,md2.mesh.x,md2.mesh.y)]; Index: sm/trunk-jpl/src/m/solve/MatlabProcessPatch.m =================================================================== --- /issm/trunk-jpl/src/m/solve/MatlabProcessPatch.m (revision 16469) +++ (revision ) @@ -1,65 +1,0 @@ -function structure=MatlabProcessPatch(structure) -%PROCESSPATCH - create a structure from a patch -% -% Usage: -% Result=ProcessPatch(Result); - -%return if there is no field Patch -if (~isfield(structure,'Patch')), - return; -end - -%loop over steps -for i=1:length(structure), - - %Get Patch for current step - Patch=structure(i).Patch; - numvertices=structure(i).PatchVertices; - - %check that Patch is not empty - if length(Patch)==0 continue; end - - %Get number of fields; - fields=unique(Patch(:,1)); - steps=unique(Patch(:,2)); - - %parse steps - for j=1:length(steps), - - posstep=find(Patch(:,2)==steps(j)); - - %Take all the lines of the Patch for this timestep - temporarypatch=Patch(posstep,:); - time=temporarypatch(1,3); - step=temporarypatch(1,2); - - %parse fields - for i=1:length(fields), - - %get name - fieldname=EnumToString(fields(i)); - - %get line positions - pos=find(temporarypatch(:,1)==fields(i)); - - %Fill Result structure - structure(step).steps=step; - structure(step).time=time; - structure(step).(fieldname).element=temporarypatch(pos,4); - structure(step).(fieldname).interpolation=temporarypatch(pos,5); - structure(step).(fieldname).index=temporarypatch(pos,6:5+numvertices); - if structure(step).(fieldname).interpolation==P1Enum(), - structure(step).(fieldname).value=temporarypatch(pos,6+numvertices:end); - end - if structure(step).(fieldname).interpolation==P0Enum(), - structure(step).(fieldname).value=temporarypatch(pos,6+numvertices); - end - - end - end -end - -%remove fields -structure=rmfield(structure,'Patch'); -structure=rmfield(structure,'PatchVertices'); -structure=rmfield(structure,'PatchNodes'); Index: sm/trunk-jpl/src/m/solve/MatlabProcessPatch.py =================================================================== --- /issm/trunk-jpl/src/m/solve/MatlabProcessPatch.py (revision 16469) +++ (revision ) @@ -1,19 +1,0 @@ -def MatlabProcessPatch(structure): - """ - PROCESSPATCH - create a structure from a patch - - Usage: - Result=ProcessPatch(Result); - """ - - #loop over steps - for structurei in structure: - - #return if there is no field Patch - if not hasattr(structurei,'Patch'): - continue - - raise SystemError("MatlabProcessPatch not implemented in Python.") - - return structure - Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 16469) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 16470) @@ -11,6 +11,4 @@ end -%process patch if necessary -results=MatlabProcessPatch(results); function results=parseresultsfromdiskioserial(filename) % {{{ @@ -76,9 +74,5 @@ %Add result - if strcmpi(result.fieldname,'Patch'), - results(result.step).(result.fieldname)=[0 result.N]; - else - results(result.step).(result.fieldname)=NaN; - end + results(result.step).(result.fieldname)=NaN; %read next result @@ -90,19 +84,6 @@ result=ReadDataDimensions(fid); while ~isempty(result), - - %Add result - if strcmpi(result.fieldname,'Patch'), - patchdimensions=results(result.step).(result.fieldname); - results(result.step).(result.fieldname)=[patchdimensions(1)+result.M result.N]; - end - %read next result result=ReadDataDimensions(fid); -end - -%allocate patches -for i=1:length(results), - results(i).Patch=zeros(results(i).Patch(1),results(i).Patch(2)); - results(i).counter=1; %use to index into the patch end @@ -117,14 +98,5 @@ %Add result - if strcmpi(result.fieldname,'Patch'), - counter=results(result.step).counter; - counter2=counter+size(result.field,1)-1; - results(result.step).(result.fieldname)(counter:counter2,:)=result.field; - - %increment counter: - results(result.step).counter=counter2+1; - else - results(result.step).(result.fieldname)=result.field; - end + results(result.step).(result.fieldname)=result.field; %read next result Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 16469) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 16470) @@ -4,5 +4,4 @@ import results as resultsclass from MatlabFuncs import * -from MatlabProcessPatch import * def parseresultsfromdisk(filename,iosplit): @@ -18,7 +17,4 @@ else: results=parseresultsfromdiskioserial(filename) - - #process patch if necessary - results=MatlabProcessPatch(results) return results @@ -96,8 +92,5 @@ #Add result - if strcmpi(result['fieldname'],'Patch'): - setattr(results[result['step']-1],result['fieldname'],[0,result['N']]) - else: - setattr(results[result['step']-1],result['fieldname'],float('NaN')) + setattr(results[result['step']-1],result['fieldname'],float('NaN')) #read next result @@ -109,17 +102,6 @@ while result: - #Add result - if strcmpi(result['fieldname'],'Patch'): - patchdimensions=getattr(results[result['step']-1],result['fieldname']) - setattr(results[result['step']-1],result['fieldname'],[patchdimensions[0]+result['M'],result['N']]) - #read next result result=ReadDataDimensions(fid) - - #allocate patches - for result in results.itervalues(): - if 'Patch' in result: - setattr(result,'Patch',numpy.zeros((result['Patch'][0],result['Patch'][1]))) - setattr(result,'counter',0) #use to index into the patch #third pass, this time to read the real information @@ -137,12 +119,5 @@ #Add result - if strcmpi(result['fieldname'],'Patch'): - counter=results[result['step']-1].counter - counter2=counter+result['field'].shape[0]-1 - getattr(results[result['step']-1],result['fieldname'])[counter:counter2,:]=result['field'] - #increment counter: - results[result['step']-1].counter=counter2+1 - else: - setattr(results[result['step']-1],result['fieldname'],result['field']) + setattr(results[result['step']-1],result['fieldname'],result['field']) #read next result Index: /issm/trunk-jpl/src/m/solve/process_solve_options.m =================================================================== --- /issm/trunk-jpl/src/m/solve/process_solve_options.m (revision 16469) +++ /issm/trunk-jpl/src/m/solve/process_solve_options.m (revision 16470) @@ -10,5 +10,5 @@ solution_type=getfieldvalue(options,'solution_type'); if ~ismember(solution_type,[StressbalanceSolutionEnum(),MasstransportSolutionEnum(),ThermalSolutionEnum(),... - SteadystateSolutionEnum(),TransientSolutionEnum(),EnthalpySolutionEnum(),... + SteadystateSolutionEnum(),TransientSolutionEnum(),... BalancethicknessSolutionEnum(),BalancethicknessSoftSolutionEnum(),... BalancevelocitySolutionEnum(),BedSlopeSolutionEnum(),... Index: /issm/trunk-jpl/src/m/solve/process_solve_options.py =================================================================== --- /issm/trunk-jpl/src/m/solve/process_solve_options.py (revision 16469) +++ /issm/trunk-jpl/src/m/solve/process_solve_options.py (revision 16470) @@ -19,5 +19,5 @@ solution_type=options.getfieldvalue('solution_type') if solution_type not in (StressbalanceSolutionEnum(),MasstransportSolutionEnum(),ThermalSolutionEnum(),\ - SteadystateSolutionEnum(),TransientSolutionEnum(),EnthalpySolutionEnum(),\ + SteadystateSolutionEnum(),TransientSolutionEnum(),\ BalancethicknessSolutionEnum(),BalancevelocitySolutionEnum(),\ BedSlopeSolutionEnum(),SurfaceSlopeSolutionEnum(),\