Changeset 16470


Ignore:
Timestamp:
10/20/13 20:31:07 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: merging thermal_core with enthalpy_core (removing EnthalpySolution), and replacing InputToResult by OutputRequest

Location:
issm/trunk-jpl/src
Files:
2 added
11 deleted
65 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/CMakeLists.txt

    r16289 r16470  
    9696                                        ./classes/Vertex.cpp
    9797                                        ./classes/Hook.cpp
    98                                         ./classes/Patch.cpp
    9998                                        ./classes/ElementResults/DoubleElementResult.cpp
    10099                                        ./classes/ElementResults/TriaP1ElementResult.cpp
    … …  
    279278                                           ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp
    280279                                           ./analyses/thermal_core.cpp
    281                                            ./analyses/enthalpy_core.cpp
    282280                                           ./solutionsequences/solutionsequence_thermal_nonlinear.cpp)
    283281#}}}

  • issm/trunk-jpl/src/c/Makefile.am

    r16442 r16470  
    6363                                        ./classes/Hook.h\
    6464                                        ./classes/Hook.cpp\
    65                                         ./classes/Patch.h\
    66                                         ./classes/Patch.cpp\
    6765                                        ./classes/ElementResults/ElementResultLocal.h\
    6866                                        ./classes/ElementResults/DoubleElementResult.h\
    … …  
    426424                                           ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\
    427425                                           ./analyses/thermal_core.cpp\
    428                                            ./analyses/enthalpy_core.cpp\
    429426                                           ./solutionsequences/solutionsequence_thermal_nonlinear.cpp
    430427#}}}

  • issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp

    r16442 r16470  
    4949
    5050                case ThermalSolutionEnum:
    51                         numanalyses=2;
     51                        numanalyses=3;
    5252                        analyses=xNew<int>(numanalyses);
    5353                        analyses[0]=ThermalAnalysisEnum;
    5454                        analyses[1]=MeltingAnalysisEnum;
    55                         break;
    56 
    57                 case EnthalpySolutionEnum:
    58                         numanalyses=1;
    59                         analyses=xNew<int>(numanalyses);
    60                         analyses[0]=EnthalpyAnalysisEnum;
     55                        analyses[2]=EnthalpyAnalysisEnum;
    6156                        break;
    6257

  • issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp

    r16181 r16470  
    4040                        #ifdef _HAVE_THERMAL_
    4141                        solutioncore=&thermal_core;
    42                         #else
    43                         _error_("ISSM was not compiled with thermal capabilities. Exiting");
    44                         #endif
    45                         break;
    46                 case EnthalpySolutionEnum:
    47                         #ifdef _HAVE_THERMAL_
    48                         solutioncore=&enthalpy_core;
    4942                        #else
    5043                        _error_("ISSM was not compiled with thermal capabilities. Exiting");

  • issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp

    r15849 r16470  
    3434        if(save_results){
    3535                if(VerboseSolution()) _printf0_("   saving results\n");
    36                 InputToResultx(femmodel,AdjointEnum);
     36                const char* outputs [] = {"Adjoint"};
     37                femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1);
    3738        }
    3839}

  • issm/trunk-jpl/src/c/analyses/adjointstressbalance_core.cpp

    r15849 r16470  
    3535
    3636        /*Save results*/
    37         if(save_results){
     37        if(save_results || true){
    3838                if(VerboseSolution()) _printf0_("   saving results\n");
    39                 InputToResultx(femmodel,AdjointxEnum);
    40                 InputToResultx(femmodel,AdjointyEnum);
    4139                if (isFS){
    42                         InputToResultx(femmodel,AdjointzEnum);
    43                         InputToResultx(femmodel,AdjointpEnum);
     40                        const char* outputs [] = {"Adjointx","Adjointy","Adjointz","Adjointp"};
     41                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],4);
     42                }
     43                else{
     44                        const char* outputs [] = {"Adjointx","Adjointy"};
     45                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2);
    4446                }
    4547        }

  • issm/trunk-jpl/src/c/analyses/analyses.h

    r16442 r16470  
    2323void hydrology_core(FemModel* femmodel);
    2424void thermal_core(FemModel* femmodel);
    25 void enthalpy_core(FemModel* femmodel);
    2625void surfaceslope_core(FemModel* femmodel);
    2726void bedslope_core(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp

    r15849 r16470  
    2626        if(save_results){
    2727                if(VerboseSolution()) _printf0_("   saving results\n");
    28                 InputToResultx(femmodel,ThicknessEnum);
     28                const char* outputs [] = {"Thickness"};
     29                femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1);
    2930        }
    3031

  • issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp

    r16144 r16470  
    3131        if(save_results){
    3232                if(VerboseSolution()) _printf0_("   saving results\n");
    33                 InputToResultx(femmodel,SurfaceSlopeXEnum);
    34                 InputToResultx(femmodel,SurfaceSlopeYEnum);
    35                 InputToResultx(femmodel,VelEnum);
     33                const char* outputs [] = {"SurfaceSlopeX","SurfaceSlopeY","Vel"};
     34                femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3);
    3635        }
    3736

  • issm/trunk-jpl/src/c/analyses/bedslope_core.cpp

    r16382 r16470  
    3232        if(save_results){
    3333                if(VerboseSolution()) _printf0_("   saving results\n");
    34                 InputToResultx(femmodel,BedSlopeXEnum);
    35                 if(meshtype!=Mesh2DverticalEnum) InputToResultx(femmodel,BedSlopeYEnum);
     34                if(meshtype!=Mesh2DverticalEnum){
     35                        const char* outputs [] = {"BedSlopeX"};
     36                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1);
     37                }
     38                else{
     39                        const char* outputs [] = {"BedSlopeX","BedSlopeY"};
     40                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2);
     41                }
    3642        }
    3743

  • issm/trunk-jpl/src/c/analyses/hydrology_core.cpp

    r16219 r16470  
    6868                        if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
    6969                                if(VerboseSolution()) _printf0_("   saving results \n");
    70                                 InputToResultx(femmodel,WatercolumnEnum);
    71                                 InputToResultx(femmodel,HydrologyWaterVxEnum);
    72                                 InputToResultx(femmodel,HydrologyWaterVyEnum);
     70                                const char* outputs [] = {"Watercolumn","HydrologyWaterVx","HydrologyWaterVy"};
     71                                femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3);
    7372
    7473                                /*unload results*/
    … …  
    9291                                InputToResultx(femmodel,SedimentHeadResidualEnum);
    9392                                if(isefficientlayer){
    94                                         InputToResultx(femmodel,EplHeadEnum);
    95                                         InputToResultx(femmodel,HydrologydcMaskEplactiveEnum);
     93                                        const char* outputs [] = {"SedimentHead","SedimentHeadResidual","EplHead","HydrologydcMaskEplactive"};
     94                                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],4);
     95                                }
     96                                else{
     97                                        const char* outputs [] = {"SedimentHead","SedimentHeadResidual"};
     98                                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2);
    9699                                }
    97100                                /*unload results*/

  • issm/trunk-jpl/src/c/analyses/masstransport_core.cpp

    r16464 r16470  
    3232        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
    3333        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     34        femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
    3435        femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
    35         femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
    3636        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
    3737

  • issm/trunk-jpl/src/c/analyses/steadystate_core.cpp

    r16363 r16470  
    3333        IssmDouble  reltol;
    3434        int         numoutputs        = 0;
    35         char**      requested_outputs = NULL;
     35        char** requested_outputs = NULL;
    3636
    3737        /* recover parameters:*/
    3838        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    3939        femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum);
    40         femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
    4140        femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    4241        femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum);
    4342        femmodel->parameters->SetParam(false,SaveResultsEnum);
     43        femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
    4444        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
    4545
    … …  
    5151                if(VerboseSolution()) _printf0_("   computing temperature and velocity for step: " << step << "\n");
    5252                #ifdef _HAVE_THERMAL_
    53                 if(isenthalpy==0){
    54                         thermal_core(femmodel);
    55                         femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
    56                         GetSolutionFromInputsx(&tg,femmodel);
    57                 }
    58                 else{
    59                         enthalpy_core(femmodel);
    60                         GetSolutionFromInputsx(&tg,femmodel);
    61                 }
     53                thermal_core(femmodel);
     54                femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/
     55                GetSolutionFromInputsx(&tg,femmodel);
    6256                #else
    6357                _error_("ISSM was not compiled with thermal capabilities. Exiting");
    … …  
    8579        if(save_results){
    8680                if(VerboseSolution()) _printf0_("   saving results\n");
    87                 InputToResultx(femmodel,VxEnum);
    88                 InputToResultx(femmodel,VyEnum);
    89                 InputToResultx(femmodel,VzEnum);
    90                 InputToResultx(femmodel,VelEnum);
    91                 InputToResultx(femmodel,PressureEnum);
    92                 InputToResultx(femmodel,TemperatureEnum);
    93                 if(isenthalpy)  InputToResultx(femmodel,WaterfractionEnum);
    94                 if(isenthalpy)  InputToResultx(femmodel,EnthalpyEnum);
    95                 if(isenthalpy)  InputToResultx(femmodel,WatercolumnEnum);
    96                 //if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
    97                 InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
    98                 femmodel->RequestedOutputsx(requested_outputs,numoutputs);
     81                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    9982        }
    10083

  • issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp

    r15849 r16470  
    1414        /*parameters: */
    1515        bool save_results;
     16        int  meshtype;
    1617
    1718        /*Recover some parameters: */
    1819        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     20        femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
    1921
    2022        if(VerboseSolution()) _printf0_("computing slope...\n");
    … …  
    2830        if(save_results){
    2931                if(VerboseSolution()) _printf0_("saving results:\n");
    30                 InputToResultx(femmodel,SurfaceSlopeXEnum);
    31                 InputToResultx(femmodel,SurfaceSlopeYEnum);
     32                if(meshtype!=Mesh2DverticalEnum){
     33                        const char* outputs [] = {"SurfaceSlopeX"};
     34                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],1);
     35                }
     36                else{
     37                        const char* outputs [] = {"SurfaceSlopeX","SurfaceSlopeY"};
     38                        femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],2);
     39                }
    3240        }
    3341

  • issm/trunk-jpl/src/c/analyses/thermal_core.cpp

    r16272 r16470  
    1313
    1414        /*intermediary*/
    15         bool   save_results;
    16         bool   dakota_analysis  = false;
    17         int    solution_type;
     15        bool   save_results,isenthalpy;
     16        bool   dakota_analysis;
     17        int    solution_type,numoutputs;
     18        char** requested_outputs = NULL;
    1819
    19         //first recover parameters common to all solutions
     20        /*first recover parameters common to all solutions*/
    2021        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    2122        femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
    2223        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     24        femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
     25        femmodel->parameters->FindParam(&numoutputs,ThermalNumRequestedOutputsEnum);
     26        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,ThermalRequestedOutputsEnum);
    2327
    2428        if(dakota_analysis && solution_type!=TransientSolutionEnum){
    … …  
    3337        }
    3438
    35         if(VerboseSolution()) _printf0_("   setting basal Dirichlet boundary conditions\n");
    36         femmodel->UpdateBasalConstraintsEnthalpyx();
     39        if(isenthalpy){
    3740
    38         if(VerboseSolution()) _printf0_("   computing temperatures\n");
    39         femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
    40         solutionsequence_thermal_nonlinear(femmodel);
     41                if(VerboseSolution()) _printf0_("   computing enthalpy\n");
     42                femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
     43                solutionsequence_nonlinear(femmodel,true);
    4144
    42         if(VerboseSolution()) _printf0_("   computing melting\n");
    43         femmodel->SetCurrentConfiguration(MeltingAnalysisEnum);
    44         solutionsequence_linear(femmodel);
     45                /*transfer enthalpy to enthalpy picard for the next step: */
     46                InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
     47
     48                /*Post process*/
     49                PostprocessingEnthalpyx(femmodel);
     50        }
     51        else{
     52                if(VerboseSolution()) _printf0_("   setting basal Dirichlet boundary conditions\n");
     53                femmodel->UpdateBasalConstraintsEnthalpyx();
     54
     55                if(VerboseSolution()) _printf0_("   computing temperatures\n");
     56                femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
     57                solutionsequence_thermal_nonlinear(femmodel);
     58
     59                if(VerboseSolution()) _printf0_("   computing melting\n");
     60                femmodel->SetCurrentConfiguration(MeltingAnalysisEnum);
     61                solutionsequence_linear(femmodel);
     62        }
    4563
    4664        if(save_results){
    4765                if(VerboseSolution()) _printf0_("   saving results\n");
    48                 InputToResultx(femmodel,TemperatureEnum);
    49                 InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
     66                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    5067        }
    5168}

  • issm/trunk-jpl/src/c/analyses/transient_core.cpp

    r16363 r16470  
    2222        int    i;
    2323        IssmDouble starttime,finaltime,dt,yts;
    24         bool   isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isenthalpy,isdelta18o,isgia;
     24        bool   isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isdelta18o,isgia;
    2525        bool   save_results,dakota_analysis;
    2626        bool   time_adapt=false;
    … …  
    4949        femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
    5050        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    51         femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    5251        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    5352        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    … …  
    116115                        if(VerboseSolution()) _printf0_("   computing temperatures\n");
    117116                        #ifdef _HAVE_THERMAL_
    118                         if(isenthalpy==0){
    119                                 thermal_core(femmodel);
    120                         }
    121                         else{
    122                                 enthalpy_core(femmodel);
    123                                 PostprocessingEnthalpyx(femmodel);
    124                         }
     117                        thermal_core(femmodel);
    125118                        #else
    126119                        _error_("ISSM was not compiled with thermal capabilities. Exiting");
    … …  
    152145                        #endif
    153146                        if(save_results){
    154                                 InputToResultx(femmodel,SurfaceEnum);
    155                                 InputToResultx(femmodel,BedEnum);
    156                                 InputToResultx(femmodel,MaskGroundediceLevelsetEnum);
     147                                const char* outputs [] = {"Surface","Bed","MaskGroundediceLevelset"};
     148                                femmodel->RequestedOutputsx(&femmodel->results,(char**)&outputs[0],3);
    157149                        }
    158150                }

  • issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.cpp

    r15737 r16470  
    8080}
    8181/*}}}*/
    82 /*FUNCTION BoolElementResult::PatchFill{{{*/
    83 void BoolElementResult::PatchFill(int row, Patch* patch){
    84 
    85          /*Here, we fill the result information into the patch object. First, let's remember what is in a row
    86           * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values
    87           * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */
    88         IssmDouble IssmDoublevalue=this->value?1:0;
    89         patch->fillresultinfo(row,this->enum_type,this->step,this->time,P0Enum,&IssmDoublevalue,1);
    90 
    91 }
    92 /*}}}*/
    9382/*FUNCTION BoolElementResult::GetVectorFromResults{{{*/
    9483void BoolElementResult::GetVectorFromResults(Vector<IssmDouble>* vector,int* doflist,int* connectivitylist,int numdofs){

  • issm/trunk-jpl/src/c/classes/ElementResults/BoolElementResult.h

    r15737 r16470  
    3737                int     GetStep(void){return step;};
    3838                int     NumberOfNodalValues(void);
    39                 void    PatchFill(int row, Patch* patch);
    4039                /*}}}*/
    4140                /*BoolElementResult management: {{{*/

  • issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.cpp

    r15737 r16470  
    8080}
    8181/*}}}*/
    82 /*FUNCTION DoubleElementResult::PatchFill{{{*/
    83 void DoubleElementResult::PatchFill(int row, Patch* patch){
    84 
    85          /*Here, we fill the result information into the patch object. First, let's remember what is in a row
    86           * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values
    87           * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */
    88         patch->fillresultinfo(row,this->enum_type,this->step,this->time,P0Enum,&this->value,1);
    89 
    90 }
    91 /*}}}*/
    9282/*FUNCTION DoubleElementResult::GetVectorFromResults{{{1*/
    9383void DoubleElementResult::GetVectorFromResults(Vector<IssmDouble>* vector,int* doflist,int* connectivitylist,int numdofs){

  • issm/trunk-jpl/src/c/classes/ElementResults/DoubleElementResult.h

    r15737 r16470  
    3737                int     GetStep(void){return step;};
    3838                int     NumberOfNodalValues(void);
    39                 void    PatchFill(int row, Patch* patch);
    4039
    4140                /*DoubleElementResult management*/

  • issm/trunk-jpl/src/c/classes/ElementResults/ElementResult.h

    r15737 r16470  
    88/*Headers:*/
    99#include "../../datastructures/datastructures.h"
    10 class Patch;
    1110class Parameters;
    1211
    … …  
    1918                virtual int        GetStep(void) = 0;
    2019                virtual int        NumberOfNodalValues(void) = 0;
    21                 virtual void       PatchFill(int row, Patch *patch)=0;
    2220                virtual int        InstanceEnum() = 0;
    2321                virtual void       GetVectorFromResults(Vector <IssmDouble> *vector,int*doflist,int*connectivitylist,int numdof)=0;

  • issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.cpp

    r15737 r16470  
    8383}
    8484/*}}}*/
    85 /*FUNCTION PentaP1ElementResult::PatchFill{{{*/
    86 void PentaP1ElementResult::PatchFill(int row, Patch* patch){
    87 
    88          /*Here, we fill the result information into the patch object. First, let's remember what is in a row
    89           * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values
    90           * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */
    91         patch->fillresultinfo(row,this->enum_type,this->step,this->time,P1Enum,this->values,6);
    92 
    93 }
    94 /*}}}*/
    9585/*FUNCTION PentaP1ElementResult::GetVectorFromResults{{{*/
    9686void PentaP1ElementResult::GetVectorFromResults(Vector<IssmDouble>* vector,int* doflist,int* connectivitylist,int numdofs){

  • issm/trunk-jpl/src/c/classes/ElementResults/PentaP1ElementResult.h

    r15737 r16470  
    3636                int     GetStep(void){return step;};
    3737                int     NumberOfNodalValues(void);
    38                 void    PatchFill(int row, Patch* patch);
    3938                /*}}}*/
    4039                /*PentaP1ElementResult management: {{{*/

  • issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp

    r15737 r16470  
    8282}
    8383/*}}}*/
    84 /*FUNCTION TriaP1ElementResult::PatchFill{{{*/
    85 void TriaP1ElementResult::PatchFill(int row, Patch* patch){
    86 
    87          /*Here, we fill the result information into the patch object. First, let's remember what is in a row
    88           * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values
    89           * Here, we will supply the enum_type, step, time, interpolation and nodal_values: */
    90         patch->fillresultinfo(row,this->enum_type,this->step,this->time,P1Enum,this->values,3);
    91 
    92 }
    93 /*}}}*/
    9484/*FUNCTION TriaP1ElementResult::GetVectorFromResults{{{*/
    9585void TriaP1ElementResult::GetVectorFromResults(Vector<IssmDouble>* vector,int* doflist,int* connectivitylist,int numdofs){

  • issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.h

    r15737 r16470  
    3535                int     GetStep(void){return step;};
    3636                int     NumberOfNodalValues(void);
    37                 void    PatchFill(int row, Patch* patch);
    3837                /*}}}*/
    3938                /*TriaP1ElementResult management: {{{*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16461 r16470  
    1515class DataSet;
    1616class Parameters;
    17 class Patch;
    1817class Elements;
    1918class Loads;
    … …  
    6059                virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
    6160                virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
    62                 virtual void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes)=0;
    63                 virtual void   PatchFill(int* pcount, Patch* patch)=0;
    6461                virtual void   ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results)=0;
    6562                virtual void   ResultInterpolation(int* pinterpolation,int output_enum)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Elements.cpp

    r16469 r16470  
    1616#include "../ExternalResults/Results.h"
    1717#include "../ExternalResults/GenericExternalResult.h"
    18 #include "../Patch.h"
    1918#include "../../toolkits/toolkits.h"
    2019#include "../../shared/shared.h"
    … …  
    6160}
    6261/*}}}*/
    63 /*FUNCTION Elements::ResultsToPatch{{{*/
    64 Patch* Elements::ResultsToPatch(void){
    65 
    66         /*output: */
    67         Patch* patch=NULL;
    68 
    69         /*intermediary: */
    70         int i;
    71         int count;
    72         int numrows;
    73         int numvertices;
    74         int numnodes;
    75         int max_numvertices;
    76         int max_numnodes;
    77         int element_numvertices;
    78         int element_numrows;
    79         int element_numnodes;
    80 
    81         /*We are going to extract from the results within the elements, the desired results , and create a table
    82          * of patch information, that will hold, for each element that computed the result that
    83          * 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
    84          * at the nodes (could be different from the vertices). This will be used for visualization purposes.
    85          * For example, we could build the following patch table, for velocities:
    86          *
    87          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
    88          VxEnum 1  .5  1 P0  1 2       4.5 NaN  NaN
    89          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
    90          VzEnum 2  .8  2 P1  1 3 4     4.5 3.2  2.5
    91          * ... etc ...
    92          *
    93          * So what do we need to build the table: the maximum number of vertices included in the table,
    94          * and the maximum number of nodal values, as well as the number of rows. Once we have that,
    95          * we ask the elements to fill their own row in the table, by looping on the elememnts.
    96          *
    97          * We will use the Patch object, which will store all of the information needed, and will be able
    98          * to output itself to disk on its own. See object/Patch.h for format of this object.*/
    99 
    100         /*First, determine maximum number of vertices, nodes, and number of results: */
    101         numrows=0;
    102         numvertices=0;
    103         numnodes=0;
    104 
    105         for(i=0;i<this->Size();i++){
    106 
    107                 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
    108                 element->PatchSize(&element_numrows,&element_numvertices,&element_numnodes);
    109 
    110                 numrows+=element_numrows;
    111                 if(element_numvertices>numvertices)numvertices=element_numvertices;
    112                 if(element_numnodes>numnodes)numnodes=element_numnodes;
    113         }
    114 
    115         /*Synchronize across cluster, so as to not end up with different sizes for each patch on each cpu: */
    116         ISSM_MPI_Reduce (&numvertices,&max_numvertices,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm() );
    117         ISSM_MPI_Bcast(&max_numvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    118         numvertices=max_numvertices;
    119 
    120         ISSM_MPI_Reduce (&numnodes,&max_numnodes,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm() );
    121         ISSM_MPI_Bcast(&max_numnodes,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    122         numnodes=max_numnodes;
    123 
    124         /*Ok, initialize Patch object: */
    125         patch=new Patch(numrows,numvertices,numnodes);
    126 
    127         /*Now, go through elements, and fill the Patch object: */
    128         count=0;
    129         for(i=0;i<this->Size();i++){
    130                 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
    131                 element->PatchFill(&count,patch);
    132         }
    133 
    134         return patch;
    135 }
    136 /*}}}*/
    13762/*FUNCTION Elements::SetCurrentConfiguration{{{*/
    13863void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
    … …  
    15681        int num_procs;
    15782
    158         Patch               *patch              = NULL;
    15983        int                 *resultsenums       = NULL;
    16084        int                 *resultssizes       = NULL;
    … …  
    16488        Vector<IssmDouble> *vector = NULL;
    16589        bool                io_gather;
    166         bool                results_as_patches;
    167         int                 numberofvertices     ,numberofelements;
    168         int                 numberofresults      ,vectorsize;
     90        int                 numberofvertices,numberofelements;
     91        int                 numberofresults ,vectorsize;
    16992        int                 rank;
    17093        int                 minrank;
    … …  
    17699        /*Recover parameters: */
    177100        parameters->FindParam(&io_gather,SettingsIoGatherEnum);
    178         parameters->FindParam(&results_as_patches,SettingsResultsAsPatchesEnum);
    179101        parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
    180102        parameters->FindParam(&numberofelements,MeshNumberofelementsEnum);
    181103
    182         if(!results_as_patches){
    183                 /*No patch here, we prepare vectors*/
    184 
    185                 /*Get rank of first cpu that has results*/
    186                 if(this->Size()) rank=my_rank;
    187                 else rank=num_procs;
    188                 ISSM_MPI_Allreduce (&rank,&minrank,1,ISSM_MPI_INT,ISSM_MPI_MIN,IssmComm::GetComm());
    189 
    190                 /*see what the first element of this partition has in stock (this is common to all partitions)*/
    191                 if(my_rank==minrank){
    192                         if(this->Size()==0) _error_("Cannot write results because there is no element??");
    193                         Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(0));
    194                         element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults);
     104        /*Get rank of first cpu that has results*/
     105        if(this->Size()) rank=my_rank;
     106        else rank=num_procs;
     107        ISSM_MPI_Allreduce (&rank,&minrank,1,ISSM_MPI_INT,ISSM_MPI_MIN,IssmComm::GetComm());
     108
     109        /*see what the first element of this partition has in stock (this is common to all partitions)*/
     110        if(my_rank==minrank){
     111                if(this->Size()==0) _error_("Cannot write results because there is no element??");
     112                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(0));
     113                element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults);
     114        }
     115        ISSM_MPI_Bcast(&numberofresults,1,ISSM_MPI_INT,minrank,IssmComm::GetComm());
     116
     117        /*Get out if there is no results. Otherwise broadcast info*/
     118        if(!numberofresults) return;
     119        if(my_rank!=minrank){
     120                resultsenums=xNew<int>(numberofresults);
     121                resultssizes=xNew<int>(numberofresults);
     122                resultstimes=xNew<IssmDouble>(numberofresults);
     123                resultssteps=xNew<int>(numberofresults);
     124        }
     125        ISSM_MPI_Bcast(resultsenums,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
     126        ISSM_MPI_Bcast(resultssizes,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
     127        ISSM_MPI_Bcast(resultstimes,numberofresults,ISSM_MPI_DOUBLE,minrank,IssmComm::GetComm());
     128        ISSM_MPI_Bcast(resultssteps,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
     129
     130        /*Loop over all results and get nodal vector*/
     131        for(int i=0;i<numberofresults;i++){
     132
     133                /*Get vector for result number i*/
     134                if(resultssizes[i]==P1Enum)      vectorsize=numberofvertices;
     135                else if(resultssizes[i]==P0Enum) vectorsize=numberofelements;
     136                else _error_("Unkown result size: " << EnumToStringx(resultssizes[i]));
     137                vector=new Vector<IssmDouble>(vectorsize);
     138
     139                for(int j=0;j<this->Size();j++){
     140                        Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(j));
     141                        element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]);
    195142                }
    196                 ISSM_MPI_Bcast(&numberofresults,1,ISSM_MPI_INT,minrank,IssmComm::GetComm());
    197 
    198                 /*Get out if there is no results. Otherwise broadcast info*/
    199                 if(!numberofresults) return;
    200                 if(my_rank!=minrank){
    201                         resultsenums=xNew<int>(numberofresults);
    202                         resultssizes=xNew<int>(numberofresults);
    203                         resultstimes=xNew<IssmDouble>(numberofresults);
    204                         resultssteps=xNew<int>(numberofresults);
     143                vector->Assemble();
     144
     145                /*Serialize and add to results*/
     146                vector_serial=vector->ToMPISerial();
     147                results->DeleteResult(resultsenums[i],resultssteps[i]);
     148                if(my_rank==0){
     149                        /*No need to add this vector for all cpus*/
     150#ifdef _HAVE_ADOLC_
     151                        IssmPDouble* vector_serial_passive=xNew<IssmPDouble>(vectorsize);
     152                        for(int k=0;k<vectorsize;k++)vector_serial_passive[k]=reCast<IssmPDouble>(vector_serial[k]);
     153                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i]));
     154                        xDelete<IssmPDouble>(vector_serial_passive);
     155#else
     156                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i]));
     157#endif
    205158                }
    206                 ISSM_MPI_Bcast(resultsenums,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
    207                 ISSM_MPI_Bcast(resultssizes,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
    208                 ISSM_MPI_Bcast(resultstimes,numberofresults,ISSM_MPI_DOUBLE,minrank,IssmComm::GetComm());
    209                 ISSM_MPI_Bcast(resultssteps,numberofresults,ISSM_MPI_INT,minrank,IssmComm::GetComm());
    210 
    211                 /*Loop over all results and get nodal vector*/
    212                 for(int i=0;i<numberofresults;i++){
    213 
    214                         /*Get vector for result number i*/
    215                         if(resultssizes[i]==P1Enum)      vectorsize=numberofvertices;
    216                         else if(resultssizes[i]==P0Enum) vectorsize=numberofelements;
    217                         else _error_("Unkown result size: " << EnumToStringx(resultssizes[i]));
    218                         vector=new Vector<IssmDouble>(vectorsize);
    219 
    220                         for(int j=0;j<this->Size();j++){
    221                                 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(j));
    222                                 element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]);
    223                         }
    224                         vector->Assemble();
    225 
    226                         /*Serialize and add to results*/
    227                         vector_serial=vector->ToMPISerial();
    228                         results->DeleteResult(resultsenums[i],resultssteps[i]);
    229                         if(my_rank==0){
    230                                 /*No need to add this vector for all cpus*/
    231                                 #ifdef _HAVE_ADOLC_
    232                                 IssmPDouble* vector_serial_passive=xNew<IssmPDouble>(vectorsize);
    233                                 for(int k=0;k<vectorsize;k++)vector_serial_passive[k]=reCast<IssmPDouble>(vector_serial[k]);
    234                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i]));
    235                                 xDelete<IssmPDouble>(vector_serial_passive);
    236                                 #else
    237                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i]));
    238                                 #endif
    239                         }
    240 
    241                         /*clean up*/
    242                         delete vector;
    243                         xDelete<IssmDouble>(vector_serial);
    244                 }
    245         }
    246         else{
    247                 /*create patch object out of all results in this dataset: */
    248                 patch=this->ResultsToPatch();
    249 
    250                 /*Gather onto master cpu 0, if needed: */
    251                 if(io_gather)patch->Gather();
    252 
    253                 /*create result object and add to results dataset:*/
    254                 if (patch->maxvertices && patch->maxnodes){
    255                         results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchVerticesEnum,patch->maxvertices,1,0));
    256                         results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchNodesEnum,   patch->maxnodes,1,0));
    257                         #ifdef _HAVE_ADOLC_
    258                         IssmPDouble* values_passive=xNew<IssmPDouble>(patch->numrows*patch->numcols);
    259                         for(int k=0;k<(patch->numrows*patch->numcols);k++)values_passive[k]=reCast<IssmPDouble>(patch->values[k]);
    260                         results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum,  values_passive,patch->numrows,patch->numcols,1,0));
    261                         xDelete<IssmPDouble>(values_passive);
    262                         #else
    263                         results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum,  patch->values,patch->numrows,patch->numcols,1,0));
    264                         #endif
    265                 }
     159
     160                /*clean up*/
     161                delete vector;
     162                xDelete<IssmDouble>(vector_serial);
    266163        }
    267164
    … …  
    271168        xDelete<int>(resultssteps);
    272169        xDelete<IssmDouble>(resultstimes);
    273         delete patch;
    274170}
    275171/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Elements.h

    r15128 r16470  
    99class Loads;
    1010class Nodes;
    11 class Patch;
    1211class Results;
    1312
    … …  
    3029                void   SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
    3130                void   ToResults(Results* results,Parameters* parameters);
    32                 Patch* ResultsToPatch(void);
    3331                int    NumberOfElements(void);
    3432                void   InputDuplicate(int input_enum,int output_enum);

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

    r16468 r16470  
    28822882}
    28832883/*}}}*/
    2884 /*FUNCTION Penta::PatchFill{{{*/
    2885 void  Penta::PatchFill(int* pcount, Patch* patch){
    2886 
    2887         int i,count;
    2888         int vertices_ids[6];
    2889 
    2890         /*recover pointer: */
    2891         count=*pcount;
    2892 
    2893         /*will be needed later: */
    2894         for(i=0;i<6;i++) vertices_ids[i]=vertices[i]->Id(); //vertices id start at column 3 of the patch.
    2895 
    2896         for(i=0;i<this->results->Size();i++){
    2897                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2898 
    2899                 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
    2900                  *it to the result object, to fill the rest: */
    2901                 patch->fillelementinfo(count,this->sid+1,vertices_ids,6);
    2902                 elementresult->PatchFill(count,patch);
    2903 
    2904                 /*increment counter: */
    2905                 count++;
    2906         }
    2907 
    2908         /*Assign output pointers:*/
    2909         *pcount=count;
    2910 }/*}}}*/
    2911 /*FUNCTION Penta::PatchSize{{{*/
    2912 void  Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
    2913 
    2914         int     i;
    2915         int     numrows       = 0;
    2916         int     numnodes      = 0;
    2917         int     temp_numnodes = 0;
    2918 
    2919         /*Go through all the results objects, and update the counters: */
    2920         for (i=0;i<this->results->Size();i++){
    2921                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2922                 /*first, we have one more result: */
    2923                 numrows++;
    2924                 /*now, how many vertices and how many nodal values for this result? :*/
    2925                 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
    2926                 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    2927         }
    2928 
    2929         /*Assign output pointers:*/
    2930         *pnumrows=numrows;
    2931         *pnumvertices=NUMVERTICES;
    2932         *pnumnodes=numnodes;
    2933 }
    2934 /*}}}*/
    29352884/*FUNCTION Penta::PositiveDegreeDay{{{*/
    29362885void  Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){
    … …  
    34363385        bool       dakota_analysis;
    34373386        bool       isFS;
    3438         IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
    34393387        int        numnodes;
    34403388        int*       penta_node_ids = NULL;
    … …  
    34443392        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    34453393        iomodel->Constant(&isFS,FlowequationIsFSEnum);
    3446         iomodel->Constant(&beta,MaterialsBetaEnum);
    3447         iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
    3448         iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
    3449         iomodel->Constant(&meltingpoint,MaterialsMeltingpointEnum);
    3450         iomodel->Constant(&latentheat,MaterialsLatentheatEnum);
    34513394
    34523395        /*Checks if debuging*/
    … …  
    36363579                        }
    36373580                        break;
    3638 
    3639 
    3640                 case EnthalpyAnalysisEnum:
    3641                         /*Initialize mesh velocity*/
    3642                         for(i=0;i<6;i++)nodeinputs[i]=0;
    3643                         if (iomodel->Data(TemperatureEnum) && iomodel->Data(WaterfractionEnum) && iomodel->Data(PressureEnum)) {
    3644                                 for(i=0;i<6;i++){
    3645                                         if(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1] < meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]){
    3646                                                 nodeinputs[i]=heatcapacity*(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1]-referencetemperature);
    3647                                         }
    3648                                         else nodeinputs[i]=heatcapacity*
    3649                                          (meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]-referencetemperature)
    3650                                                 +latentheat*iomodel->Data(WaterfractionEnum)[penta_vertex_ids[i]-1];
    3651                                 }
    3652                                 this->inputs->AddInput(new PentaInput(EnthalpyEnum,nodeinputs,P1Enum));
    3653                         }
    3654                         else _error_("temperature and waterfraction required for the enthalpy solution");
    3655                         break;
    3656 
    36573581                default:
    36583582                        /*No update for other solution types*/

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

    r16461 r16470  
    103103                void   ResultInterpolation(int* pinterpolation,int output_enum);
    104104                void   ResultToVector(Vector<IssmPDouble>* vector,int output_enum);
    105                 void   PatchFill(int* pcount, Patch* patch);
    106                 void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
    107105                void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
    108106                void   ResetCoordinateSystem(void);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16461 r16470  
    117117                void        RequestedOutput(int output_enum,int step,IssmDouble time){_error_("not implemented yet");};
    118118                void        ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results){_error_("not implemented yet");};
    119                 void        PatchFill(int* pcount, Patch* patch){_error_("not implemented yet");};
    120                 void        PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){_error_("not implemented yet");};
    121119                void        ResultInterpolation(int* pinterpolation,int output_enum){_error_("not implemented");};
    122120                void        ResultToVector(Vector<IssmPDouble>* vector,int output_enum){_error_("not implemented");};

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

    r16468 r16470  
    24492449        if(found)*pvalue=value;
    24502450        return found;
    2451 }
    2452 /*}}}*/
    2453 /*FUNCTION Tria::PatchFill{{{*/
    2454 void  Tria::PatchFill(int* prow, Patch* patch){
    2455 
    2456         int i,row;
    2457         int vertices_ids[3];
    2458 
    2459         /*recover pointer: */
    2460         row=*prow;
    2461 
    2462         for(i=0;i<3;i++) vertices_ids[i]=vertices[i]->Id(); //vertices id start at column 3 of the patch.
    2463 
    2464         for(i=0;i<this->results->Size();i++){
    2465                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2466 
    2467                 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
    2468                  *it to the result object, to fill the rest: */
    2469                 patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
    2470                 elementresult->PatchFill(row,patch);
    2471 
    2472                 /*increment rower: */
    2473                 row++;
    2474         }
    2475 
    2476         /*Assign output pointers:*/
    2477         *prow=row;
    2478 }
    2479 /*}}}*/
    2480 /*FUNCTION Tria::PatchSize{{{*/
    2481 void  Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
    2482 
    2483         int     i;
    2484         int     numrows       = 0;
    2485         int     numnodes      = 0;
    2486         int     temp_numnodes = 0;
    2487 
    2488         /*Go through all the results objects, and update the counters: */
    2489         for (i=0;i<this->results->Size();i++){
    2490                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2491                 /*first, we have one more result: */
    2492                 numrows++;
    2493                 /*now, how many vertices and how many nodal values for this result? :*/
    2494                 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
    2495                 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    2496         }
    2497 
    2498         /*Assign output pointers:*/
    2499         *pnumrows=numrows;
    2500         *pnumvertices=NUMVERTICES;
    2501         *pnumnodes=numnodes;
    25022451}
    25032452/*}}}*/

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

    r16461 r16470  
    110110                void        ResultInterpolation(int* pinterpolation,int output_enum);
    111111                void        ResultToVector(Vector<IssmPDouble>* vector,int output_enum);
    112                 void        PatchFill(int* pcount, Patch* patch);
    113                 void        PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
    114112                void        ResetCoordinateSystem(void);
    115113                void          SmbGradients();

  • issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h

    r16469 r16470  
    9393                }
    9494                /*}}}*/
    95                 GenericExternalResult(int in_id, char* in_enum_string,ResultType in_value,int in_step, IssmDouble in_time){ /*{{{*/
     95                GenericExternalResult(int in_id,const char* in_enum_string,ResultType in_value,int in_step, IssmDouble in_time){ /*{{{*/
    9696                        id        = in_id;
    9797                        enum_type = -1;

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r16468 r16470  
    533533        IssmDouble  time;
    534534        IssmDouble  double_result;
    535         char      *output_string = NULL;
     535        const char *output_string = NULL;
    536536
    537537        /*recover results*/

  • issm/trunk-jpl/src/c/classes/classes.h

    r16388 r16470  
    117117#include "./DofIndexing.h"
    118118#include "./IoModel.h"
    119 #include "./Patch.h"
    120119#include "./Update.h"
    121120#include "./FemModel.h"

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

    r16388 r16470  
    7575        parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum));
    7676        parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
    77         parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsAsPatchesEnum));
    7877        parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
    7978        parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
    … …  
    152151        iomodel->DeleteData(&requestedoutputs,numoutputs,MasstransportRequestedOutputsEnum);
    153152
     153        iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
     154        parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
     155        if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
     156        iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
     157
    154158        /*Deal with mass flux segments: {{{*/
    155159        iomodel->FetchData(&qmu_mass_flux_present,QmuMassFluxSegmentsPresentEnum);

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

    r16291 r16470  
    1212void    UpdateElementsEnthalpy(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){
    1313
     14        bool dakota_analysis;
     15        bool isenthalpy;
     16
    1417        /*Now, is the model 3d? otherwise, do nothing: */
    1518        if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
     19
     20        /*Is enthalpy requested?*/
     21        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
     22        if(!isenthalpy) return;
    1623
    1724        /*Fetch data needed: */
    … …  
    2835        }
    2936
    30         bool dakota_analysis;
    3137        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    3238
    … …  
    4753        iomodel->FetchDataToInput(elements,TemperatureEnum);
    4854        iomodel->FetchDataToInput(elements,WaterfractionEnum);
     55        iomodel->FetchDataToInput(elements,EnthalpyEnum);
    4956        iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
    5057        iomodel->FetchDataToInput(elements,WatercolumnEnum);

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

    r16291 r16470  
    99#include "../ModelProcessorx.h"
    1010
    11 void    CreateNodesThermal(Nodes** pnodes, IoModel* iomodel){
     11void CreateNodesThermal(Nodes** pnodes, IoModel* iomodel){
    1212
    1313        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);

  • issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp

    r16388 r16470  
    88#include "../../classes/classes.h"
    99
    10 IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,char* output_string){
     10IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,const char* output_string){
    1111
    12         int i;
    13         Definition* definition=NULL;
    14         DataSet* output_definitions=NULL;
    15         IssmDouble return_value;
     12        Definition *definition         = NULL;
     13        DataSet    *output_definitions = NULL;
     14        IssmDouble  return_value;
    1615
    1716        /*Ok, go find the output definitions dataset in the parameters, where our responses are hiding: */
    … …  
    1918
    2019        /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_string: */
    21         for (i=0;i<output_definitions->Size();i++){
     20        for(int i=0;i<output_definitions->Size();i++){
    2221               
    2322                definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));

  • issm/trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.h

    r16388 r16470  
    88
    99/* local prototypes: */
    10 IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,char* output_string);
     10IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,const char* output_string);
    1111
    1212#endif  /* _OUTPUTDEFINITIONSRESPONSEXX_H */

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r16442 r16470  
    5050        StressbalanceIsnewtonEnum,
    5151        StressbalanceMaxiterEnum,
    52         StressbalanceNumRequestedOutputsEnum,
    5352        StressbalancePenaltyFactorEnum,
    5453        StressbalanceReferentialEnum,
    5554        StressbalanceReltolEnum,
     55        StressbalanceNumRequestedOutputsEnum,
    5656        StressbalanceRequestedOutputsEnum,
    5757        StressbalanceRestolEnum,
    … …  
    228228        SettingsLowmemEnum,
    229229        SettingsOutputFrequencyEnum,
    230         SettingsResultsAsPatchesEnum,
    231230        SettingsWaitonlockEnum,
    232231        SurfaceforcingsDelta18oEnum,
    … …  
    258257        SurfaceforcingsBNegEnum,
    259258        ThermalIsenthalpyEnum,
    260         ThermalIsdynamicbasalspcEnum,
     259        ThermalIsdynamicbasalspcEnum,
    261260        ThermalMaxiterEnum,
    262261        ThermalPenaltyFactorEnum,
    … …  
    265264        ThermalSpctemperatureEnum,
    266265        ThermalStabilizationEnum,
     266        ThermalNumRequestedOutputsEnum,
     267        ThermalRequestedOutputsEnum,
    267268        GiaMantleViscosityEnum,
    268269        GiaLithosphereThicknessEnum,
    … …  
    308309        StressbalanceVerticalAnalysisEnum,
    309310        EnthalpyAnalysisEnum,
    310         EnthalpySolutionEnum,
    311311        FlaimAnalysisEnum,
    312312        FlaimSolutionEnum,
    … …  
    539539        IntExternalResultEnum,
    540540        JEnum,
    541         PatchEnum,
    542         PatchNodesEnum,
    543         PatchVerticesEnum,
    544541        PentaP1ElementResultEnum,
    545542        StringExternalResultEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r16442 r16470  
    5858                case StressbalanceIsnewtonEnum : return "StressbalanceIsnewton";
    5959                case StressbalanceMaxiterEnum : return "StressbalanceMaxiter";
    60                 case StressbalanceNumRequestedOutputsEnum : return "StressbalanceNumRequestedOutputs";
    6160                case StressbalancePenaltyFactorEnum : return "StressbalancePenaltyFactor";
    6261                case StressbalanceReferentialEnum : return "StressbalanceReferential";
    6362                case StressbalanceReltolEnum : return "StressbalanceReltol";
     63                case StressbalanceNumRequestedOutputsEnum : return "StressbalanceNumRequestedOutputs";
    6464                case StressbalanceRequestedOutputsEnum : return "StressbalanceRequestedOutputs";
    6565                case StressbalanceRestolEnum : return "StressbalanceRestol";
    … …  
    236236                case SettingsLowmemEnum : return "SettingsLowmem";
    237237                case SettingsOutputFrequencyEnum : return "SettingsOutputFrequency";
    238                 case SettingsResultsAsPatchesEnum : return "SettingsResultsAsPatches";
    239238                case SettingsWaitonlockEnum : return "SettingsWaitonlock";
    240239                case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o";
    … …  
    273272                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
    274273                case ThermalStabilizationEnum : return "ThermalStabilization";
     274                case ThermalNumRequestedOutputsEnum : return "ThermalNumRequestedOutputs";
     275                case ThermalRequestedOutputsEnum : return "ThermalRequestedOutputs";
    275276                case GiaMantleViscosityEnum : return "GiaMantleViscosity";
    276277                case GiaLithosphereThicknessEnum : return "GiaLithosphereThickness";
    … …  
    314315                case StressbalanceVerticalAnalysisEnum : return "StressbalanceVerticalAnalysis";
    315316                case EnthalpyAnalysisEnum : return "EnthalpyAnalysis";
    316                 case EnthalpySolutionEnum : return "EnthalpySolution";
    317317                case FlaimAnalysisEnum : return "FlaimAnalysis";
    318318                case FlaimSolutionEnum : return "FlaimSolution";
    … …  
    529529                case IntExternalResultEnum : return "IntExternalResult";
    530530                case JEnum : return "J";
    531                 case PatchEnum : return "Patch";
    532                 case PatchNodesEnum : return "PatchNodes";
    533                 case PatchVerticesEnum : return "PatchVertices";
    534531                case PentaP1ElementResultEnum : return "PentaP1ElementResult";
    535532                case StringExternalResultEnum : return "StringExternalResult";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r16442 r16470  
    5858              else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
    5959              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    60               else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    6160              else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
    6261              else if (strcmp(name,"StressbalanceReferential")==0) return StressbalanceReferentialEnum;
    6362              else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
     63              else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    6464              else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
    6565              else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
    … …  
    239239              else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
    240240              else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
    241               else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum;
    242241              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    243242              else if (strcmp(name,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum;
    … …  
    260259              else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
    261260              else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
     261              else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
    266               else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
     265              if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
    267266              else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
    268267              else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
    … …  
    279278              else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
    280279              else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
     280              else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
     281              else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum;
    281282              else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum;
    282283              else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum;
    … …  
    320321              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
    321322              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    322               else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum;
    323323              else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
    324324              else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
    … …  
    541541              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    542542              else if (strcmp(name,"J")==0) return JEnum;
    543               else if (strcmp(name,"Patch")==0) return PatchEnum;
    544               else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum;
    545               else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum;
    546543              else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
    547544              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;

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

    r16431 r16470  
    5656                                md = checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
    5757                        end
    58                         if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy) | solution==EnthalpySolutionEnum(),
     58                        if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy)
    5959                                md = checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
    6060                                md = checkfield(md,'initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices 1]);
    … …  
    106106                        WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum);
    107107                        WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);
     108
     109                        if md.thermal.isenthalpy,
     110                                tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
     111                                pos  = find(md.initialization.temperature>tpmp);
     112                                enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
     113                                enthalpy(pos) = md.materials.heatcapacity*tpmp(pos) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction(pos);
     114                                WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum());
     115                        end
    108116                end % }}}
    109117        end

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

    r16431 r16470  
    6868                                md = checkfield(md,'initialization.vz','NaN',1,'size',[md.mesh.numberofvertices])
    6969                        md = checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices])
    70                 if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy) or solution==EnthalpySolutionEnum():
     70                if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy):
    7171                        md = checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
    7272                        md = checkfield(md,'initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices])
    … …  
    8787                WriteData(fid,'data',self.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum())
    8888                WriteData(fid,'data',self.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum())
     89
     90                if md.thermal.isenthalpy:
     91                        tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
     92                        pos  = numpy.nonzero(md.initialization.temperature > tpmp)[0]
     93                        enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
     94                        enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,1) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,1)
     95                        WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum());
     96
    8997        # }}}

  • issm/trunk-jpl/src/m/classes/masstransport.m

    r16464 r16470  
    6060                        end
    6161                end % }}}
    62                 function list=defaultoutputs(self,md) % {{{
     62                function list = defaultoutputs(self,md) % {{{
    6363
    6464                        list = {'Thickness','Surface','Bed'};

  • issm/trunk-jpl/src/m/classes/settings.m

    r15643 r16470  
    88                io_gather           = 0;
    99                lowmem              = 0;
    10                 results_as_patches  = 0;
    1110                output_frequency    = 0;
    1211                waitonlock          = 0;
    … …  
    3837                        obj.output_frequency=1;
    3938
    40                         %do not use patches by default (difficult to plot)
    41                         obj.results_as_patches=0;
    42 
    4339                        %this option can be activated to load automatically the results
    4440                        %onto the model after a parallel run by waiting for the lock file
    … …  
    5551                        md = checkfield(md,'settings.io_gather','numel',[1],'values',[0 1]);
    5652                        md = checkfield(md,'settings.lowmem','numel',[1],'values',[0 1]);
    57                         md = checkfield(md,'settings.results_as_patches','numel',[1],'values',[0 1]);
    5853                        md = checkfield(md,'settings.output_frequency','numel',[1],'>=',1);
    5954                        md = checkfield(md,'settings.waitonlock','numel',[1]);
    … …  
    6560                        fielddisplay(obj,'io_gather','I/O gathering strategy for result outputs (default 1)');
    6661                        fielddisplay(obj,'lowmem','is the memory limited ? (0 or 1)');
    67                         fielddisplay(obj,'results_as_patches','provide results as patches for each element (0 or 1)');
    6862                        fielddisplay(obj,'output_frequency','frequency at which results are saved in all solutions with multiple time_steps');
    6963                        fielddisplay(obj,'waitonlock','maximum number of minutes to wait for batch results (NaN to deactivate)');
    … …  
    7872                        WriteData(fid,'object',obj,'fieldname','io_gather','format','Boolean');
    7973                        WriteData(fid,'object',obj,'fieldname','lowmem','format','Boolean');
    80                         WriteData(fid,'object',obj,'fieldname','results_as_patches','format','Boolean');
    8174                        WriteData(fid,'object',obj,'fieldname','output_frequency','format','Integer');
    8275                        if obj.waitonlock>0,

  • issm/trunk-jpl/src/m/classes/settings.py

    r15131 r16470  
    1515                self.io_gather           = 0
    1616                self.lowmem              = 0
    17                 self.results_as_patches  = 0
    1817                self.output_frequency    = 0
    1918                self.waitonlock          = 0
    … …  
    2827                string="%s\n%s"%(string,fielddisplay(self,"io_gather","I/O gathering strategy for result outputs (default 1)"))
    2928                string="%s\n%s"%(string,fielddisplay(self,"lowmem","is the memory limited ? (0 or 1)"))
    30                 string="%s\n%s"%(string,fielddisplay(self,"results_as_patches","provide results as patches for each element (0 or 1)"))
    3129                string="%s\n%s"%(string,fielddisplay(self,"output_frequency","frequency at which results are saved in all solutions with multiple time_steps"))
    3230                string="%s\n%s"%(string,fielddisplay(self,"waitonlock","maximum number of minutes to wait for batch results, or return 0"))
    … …  
    4442                self.output_frequency=1
    4543
    46                 #do not use patches by default (difficult to plot)
    47                 self.results_as_patches=0
    48 
    4944                #this option can be activated to load automatically the results
    5045                #onto the model after a parallel run by waiting for the lock file
    … …  
    5853                md = checkfield(md,'settings.io_gather','numel',[1],'values',[0,1])
    5954                md = checkfield(md,'settings.lowmem','numel',[1],'values',[0,1])
    60                 md = checkfield(md,'settings.results_as_patches','numel',[1],'values',[0,1])
    6155                md = checkfield(md,'settings.output_frequency','numel',[1],'>=',1)
    6256                md = checkfield(md,'settings.waitonlock','numel',[1])
    … …  
    6761                WriteData(fid,'object',self,'fieldname','io_gather','format','Boolean')
    6862                WriteData(fid,'object',self,'fieldname','lowmem','format','Boolean')
    69                 WriteData(fid,'object',self,'fieldname','results_as_patches','format','Boolean')
    7063                WriteData(fid,'object',self,'fieldname','output_frequency','format','Integer')
    7164                if self.waitonlock>0:

  • issm/trunk-jpl/src/m/classes/steadystate.m

    r16458 r16470  
    2525                        %Relative tolerance for the steadystate convertgence
    2626                        obj.reltol=0.01;
     27
     28                        %default output
     29                        obj.requested_outputs={'default'};
     30                end % }}}
     31                function list=defaultoutputs(self,md) % {{{
     32
     33                        list =  [md.stressbalance.defaultoutputs(md) md.thermal.defaultoutputs(md)];
     34
    2735                end % }}}
    2836                function md = checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5159                        WriteData(fid,'object',obj,'fieldname','reltol','format','Double');
    5260                        WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
    53                         WriteData(fid,'object',obj,'fieldname','requested_outputs','format','StringArray');
     61
     62                        %process requested outputs
     63                        outputs = obj.requested_outputs;
     64                        pos  = find(ismember(outputs,'default'));
     65                        if ~isempty(pos),
     66                                outputs(pos) = [];                         %remove 'default' from outputs
     67                                outputs      = [outputs defaultoutputs(obj,md)]; %add defaults
     68                        end
     69                        WriteData(fid,'data',outputs,'enum',SteadystateRequestedOutputsEnum,'format','StringArray');
    5470                end % }}}
    5571        end

  • issm/trunk-jpl/src/m/classes/steadystate.py

    r16458 r16470  
    2929                return string
    3030                #}}}
     31        def defaultoutputs(self,md): # {{{
     32
     33                return md.stressbalance.defaultoutputs(md)+md.thermal.defaultoutputs(md)
     34
     35        #}}}
    3136        def setdefaultparameters(self): # {{{
    3237               
    … …  
    3742                self.reltol=0.01
    3843
     44                #default output
     45                self.requested_outputs=['default']
    3946                return self
    4047        #}}}
    … …  
    5865                WriteData(fid,'object',self,'fieldname','reltol','format','Double')
    5966                WriteData(fid,'object',self,'fieldname','maxiter','format','Integer')
    60                 WriteData(fid,'object',self,'fieldname','requested_outputs','format','StringArray')
     67
     68                #process requested outputs
     69                outputs = self.requested_outputs
     70                indices = [i for i, x in enumerate(outputs) if x == 'default']
     71                if len(indices) > 0:
     72                        outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
     73                        outputs    =outputscopy
     74                        print outputs
     75                WriteData(fid,'data',outputs,'enum',SteadystateRequestedOutputsEnum(),'format','StringArray')
    6176        # }}}

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

    r16322 r16470  
    1414                isenthalpy        = 0;
    1515                isdynamicbasalspc = 0;
     16                requested_outputs = {};
    1617        end
    1718        methods
    … …  
    2324                                        error('constructor not supported');
    2425                        end
     26                end % }}}
     27                function list = defaultoutputs(self,md) % {{{
     28
     29                        if self.isenthalpy,
     30                                list = {'Enthalpy','Temperature','Waterfraction','Watercolumn','BasalforcingsMeltingRate'};
     31                        else
     32                                list = {'Temperature','BasalforcingsMeltingRate'};
     33                        end
     34
    2535                end % }}}
    2636                function obj = setdefaultparameters(obj) % {{{
    … …  
    4353                        %will basal boundary conditions be set dynamically
    4454                        obj.isdynamicbasalspc=0;
     55
     56                        %default output
     57                        obj.requested_outputs={'default'};
    4558                end % }}}
    4659                function md = checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5164                        md = checkfield(md,'thermal.stabilization','numel',[1],'values',[0 1 2]);
    5265                        md = checkfield(md,'thermal.spctemperature','forcing',1);
    53                         if (ismember(EnthalpyAnalysisEnum(),analyses) & (md.thermal.isenthalpy | solution==EnthalpySolutionEnum()) & md.mesh.dimension==3),
     66                        if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy & md.mesh.dimension==3),
    5467                                pos=find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN);
    5568                                replicate=repmat(md.geometry.surface-md.mesh.z,1,size(md.thermal.spctemperature,2));
    … …  
    6174                                end
    6275            end
     76
     77                 md = checkfield(md,'thermal.requested_outputs','stringrow',1);
    6378    end % }}}
    6479                function disp(obj) % {{{
    … …  
    7388                        fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    7489                        fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
     90                        fielddisplay(obj,'requested_outputs','additional outputs requested');
    7591
    7692                end % }}}
    … …  
    84100                        WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
    85101                        WriteData(fid,'object',obj,'fieldname','isdynamicbasalspc','format','Boolean');
     102
     103                        %process requested outputs
     104                        outputs = obj.requested_outputs;
     105                        pos  = find(ismember(outputs,'default'));
     106                        if ~isempty(pos),
     107                                outputs(pos) = [];                         %remove 'default' from outputs
     108                                outputs      = [outputs defaultoutputs(obj,md)]; %add defaults
     109                        end
     110                        WriteData(fid,'data',outputs,'enum',ThermalRequestedOutputsEnum(),'format','StringArray');
    86111                end % }}}
    87112        end

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

    r16274 r16470  
    2222                self.isenthalpy        = 0
    2323                self.isdynamicbasalspc = 0;
     24                self.requested_outputs = []
    2425
    2526                #set defaults
    … …  
    3637                string="%s\n%s"%(string,fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'))
    3738                string="%s\n%s"%(string,fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'))
     39                string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
    3840                return string
    3941                #}}}
     42        def defaultoutputs(self,md): # {{{
     43
     44                if self.isenthalpy:
     45                        return ['Enthalpy','Temperature','Waterfraction','Watercolumn','BasalforcingsMeltingRate']
     46                else:
     47                        return ['Temperature','BasalforcingsMeltingRate']
     48
     49        #}}}
    4050        def setdefaultparameters(self): # {{{
    4151               
    … …  
    5767                #will basal boundary conditions be set dynamically
    5868                self.isdynamicbasalspc=0;
     69
     70                #default output
     71                self.requested_outputs=['default']
    5972                return self
    6073
    … …  
    6881                md = checkfield(md,'thermal.stabilization','numel',[1],'values',[0,1,2])
    6982                md = checkfield(md,'thermal.spctemperature','forcing',1)
    70                 if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3:
     83                if EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy and md.mesh.dimension==3:
    7184                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])))
    7285                        replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1)))
    … …  
    7487                        md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1])
    7588                        md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0,1]);
     89                        md = checkfield(md,'thermal.requested_outputs','stringrow',1)
    7690
    7791                return md
    … …  
    86100                WriteData(fid,'object',self,'fieldname','isenthalpy','format','Boolean')
    87101                WriteData(fid,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
     102
     103                #process requested outputs
     104                outputs = self.requested_outputs
     105                indices = [i for i, x in enumerate(outputs) if x == 'default']
     106                if len(indices) > 0:
     107                        outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
     108                        outputs    =outputscopy
     109                WriteData(fid,'data',outputs,'enum',ThermalRequestedOutputsEnum(),'format','StringArray')
    88110        # }}}

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m

    r16181 r16470  
    5555        case ThermalSolutionEnum(),
    5656                numanalyses=2;
    57                 analyses=[ThermalAnalysisEnum();MeltingAnalysisEnum()];
    58 
    59         case EnthalpySolutionEnum(),
    60                 numanalyses=1;
    61                 analyses=[EnthalpyAnalysisEnum()];
     57                analyses=[EnthalpyAnalysisEnum;ThermalAnalysisEnum();MeltingAnalysisEnum()];
    6258
    6359        case MasstransportSolutionEnum(),

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py

    r16008 r16470  
    2020        elif solutiontype == ThermalSolutionEnum():
    2121                numanalyses=2
    22                 analyses=[ThermalAnalysisEnum(),MeltingAnalysisEnum()]
    23 
    24         elif solutiontype == EnthalpySolutionEnum():
    25                 numanalyses=1
    26                 analyses=[EnthalpyAnalysisEnum()]
     22                analyses=[EnthalpyAnalysisEnum(),ThermalAnalysisEnum(),MeltingAnalysisEnum()]
    2723
    2824        elif solutiontype == MasstransportSolutionEnum():

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r16443 r16470  
    5050def StressbalanceIsnewtonEnum(): return StringToEnum("StressbalanceIsnewton")[0]
    5151def StressbalanceMaxiterEnum(): return StringToEnum("StressbalanceMaxiter")[0]
    52 def StressbalanceNumRequestedOutputsEnum(): return StringToEnum("StressbalanceNumRequestedOutputs")[0]
    5352def StressbalancePenaltyFactorEnum(): return StringToEnum("StressbalancePenaltyFactor")[0]
    5453def StressbalanceReferentialEnum(): return StringToEnum("StressbalanceReferential")[0]
    5554def StressbalanceReltolEnum(): return StringToEnum("StressbalanceReltol")[0]
     55def StressbalanceNumRequestedOutputsEnum(): return StringToEnum("StressbalanceNumRequestedOutputs")[0]
    5656def StressbalanceRequestedOutputsEnum(): return StringToEnum("StressbalanceRequestedOutputs")[0]
    5757def StressbalanceRestolEnum(): return StringToEnum("StressbalanceRestol")[0]
    … …  
    228228def SettingsLowmemEnum(): return StringToEnum("SettingsLowmem")[0]
    229229def SettingsOutputFrequencyEnum(): return StringToEnum("SettingsOutputFrequency")[0]
    230 def SettingsResultsAsPatchesEnum(): return StringToEnum("SettingsResultsAsPatches")[0]
    231230def SettingsWaitonlockEnum(): return StringToEnum("SettingsWaitonlock")[0]
    232231def SurfaceforcingsDelta18oEnum(): return StringToEnum("SurfaceforcingsDelta18o")[0]
    … …  
    265264def ThermalSpctemperatureEnum(): return StringToEnum("ThermalSpctemperature")[0]
    266265def ThermalStabilizationEnum(): return StringToEnum("ThermalStabilization")[0]
     266def ThermalNumRequestedOutputsEnum(): return StringToEnum("ThermalNumRequestedOutputs")[0]
     267def ThermalRequestedOutputsEnum(): return StringToEnum("ThermalRequestedOutputs")[0]
    267268def GiaMantleViscosityEnum(): return StringToEnum("GiaMantleViscosity")[0]
    268269def GiaLithosphereThicknessEnum(): return StringToEnum("GiaLithosphereThickness")[0]
    … …  
    306307def StressbalanceVerticalAnalysisEnum(): return StringToEnum("StressbalanceVerticalAnalysis")[0]
    307308def EnthalpyAnalysisEnum(): return StringToEnum("EnthalpyAnalysis")[0]
    308 def EnthalpySolutionEnum(): return StringToEnum("EnthalpySolution")[0]
    309309def FlaimAnalysisEnum(): return StringToEnum("FlaimAnalysis")[0]
    310310def FlaimSolutionEnum(): return StringToEnum("FlaimSolution")[0]
    … …  
    521521def IntExternalResultEnum(): return StringToEnum("IntExternalResult")[0]
    522522def JEnum(): return StringToEnum("J")[0]
    523 def PatchEnum(): return StringToEnum("Patch")[0]
    524 def PatchNodesEnum(): return StringToEnum("PatchNodes")[0]
    525 def PatchVerticesEnum(): return StringToEnum("PatchVertices")[0]
    526523def PentaP1ElementResultEnum(): return StringToEnum("PentaP1ElementResult")[0]
    527524def StringExternalResultEnum(): return StringToEnum("StringExternalResult")[0]

  • issm/trunk-jpl/src/m/geometry/NowickiProfile.m

    r16412 r16470  
    1212delta = 0.1;          % ratio of water density and ice density -1
    1313hg    = 1;            % ice thickness at grounding line
    14 sea   = hg/(1+delta); % sea level
    1514lamda = 0.1;          % ration of deviatoric stress and water pressure
    1615beta  = 5;            % friction coefficient
    … …  
    3635for i = ceil(numel(x)/2):numel(x)
    3736        h(i) = (x(i)/(4*(delta+1)*q)+hg^(-2))^(-0.5); % ice thickness for ice shelf from (3.1)
    38         b(i) = sea-h(i)*(1/(1+delta));
     37        b(i) = -h(i)*(1/(1+delta));
    3938end

  • issm/trunk-jpl/src/m/plot/plot_unit.m

    r13009 r16470  
    5757                end
    5858
    59         %Patch plot P1
    60         case 4,
    61 
    62                 if is2d,
    63                         patch( 'Faces',elements,'Vertices',[x y],'CData',data,'FaceColor','interp','EdgeColor',edgecolor);
    64                 else
    65                         A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    66                         patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor);
    67                         patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor);
    68                         patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor);
    69                         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor);
    70                         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', data,'FaceColor','interp','EdgeColor',edgecolor);
    71                 end
    72 
    73         %Patch plot P0
    74         case 5,
    75 
    76                 if is2d,
    77                         A=elements(:,1); B=elements(:,2); C=elements(:,3);
    78                         patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', data(:),'FaceColor','flat','EdgeColor',edgecolor);
    79                 else
    80                         A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    81                         patch( 'Faces', [A B C],  'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor);
    82                         patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor);
    83                         patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor);
    84                         patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor);
    85                         patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', data,'FaceColor','flat','EdgeColor',edgecolor);
    86                 end
    87 
    8859        otherwise,
    8960                error(['case ' num2str(datatype) ' not supported']);

  • issm/trunk-jpl/src/m/plot/processdata.m

    r16325 r16470  
    3737end
    3838
    39 %Process Patch
    40 if isstruct(data)
    41         if (isprop(data,'index') & isprop(data,'value')),
    42                 if data.interpolation(1)==P1Enum(),
    43                         data=data.value;
    44                         data=data';
    45                         data=data(:);
    46                         datatype=4;
    47                 elseif data.interpolation(1)==P0Enum(),
    48                         data=data.value;
    49                         datatype=5;
    50                 else
    51                         error(['interpolation ' data.interpolation(1) ' not supported yet']);
    52                 end
    53         else
    54                 error('structure other than Patch not supported yet');
    55         end
    56 else
    57         %initialize datatype
    58         datatype=0;
    59 end
     39%initialize datatype
     40datatype=0;
    6041
    6142%get datasize
    6243datasize=size(data);
    6344
    64 %non patch processing
    65 if datatype~=4 & datatype~=5,
     45%transpose data if necessary
     46if (size(data,2) > size(data,1)),
     47        data=data';
     48end
     49datasize=size(data);
    6650
    67         %transpose data if necessary
    68         if (size(data,2) > size(data,1)),
    69                 data=data';
    70         end
    71         datasize=size(data);
     51%convert to double if necessary
     52if ~isnumeric(data);
     53        disp('processdata info message: data is not numeric (logical?). Converted to double');
     54        data=double(data);
     55end
    7256
    73         %convert to double if necessary
    74         if ~isnumeric(data);
    75                 disp('processdata info message: data is not numeric (logical?). Converted to double');
    76                 data=double(data);
    77         end
     57%check length
     58if 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))
     59        error('plotmodel error message: data not supported yet');
     60end
    7861
    79         %check length
    80         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))
    81                 error('plotmodel error message: data not supported yet');
    82         end
     62%quiver?
     63if datasize(2)>1,
     64        datatype=3;
    8365
    84         %quiver?
    85         if datasize(2)>1,
    86                 datatype=3;
    87 
    88                 %check number of columns, add zeros if necessary,
    89                 if (md.mesh.dimension==3)
    90                         if datasize(2)==2,
    91                                 data=[data, zeros(datasize(1),1)];
    92                         elseif datasize(2)~=3,
    93                                 error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot');
    94                         end
     66        %check number of columns, add zeros if necessary,
     67        if (md.mesh.dimension==3)
     68                if datasize(2)==2,
     69                        data=[data, zeros(datasize(1),1)];
     70                elseif datasize(2)~=3,
     71                        error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot');
     72                end
    9573                %elseif ((md.mesh.dimension==2) & datasize(2)~=2),
    9674                %       error('plotmodel error message: data provided should have 2 columns for quiver plot, and 1 for regular plot');
    97                 end
    9875        end
     76end
    9977
    100         %treat the case datasize(1)=6*nodes
    101         if datasize(1)==6*md.mesh.numberofvertices
    102                 %keep the only norm of data
    103                 data1=data(1:6:md.mesh.numberofvertices*6,:);
    104                 data2=data(2:6:md.mesh.numberofvertices*6,:);
    105                 data=sqrt(data1.^2+data2.^2);
    106                 datasize(1)=md.mesh.numberofvertices;
    107                 %---> go to node data
    108         end
     78%treat the case datasize(1)=6*nodes
     79if datasize(1)==6*md.mesh.numberofvertices
     80        %keep the only norm of data
     81        data1=data(1:6:md.mesh.numberofvertices*6,:);
     82        data2=data(2:6:md.mesh.numberofvertices*6,:);
     83        data=sqrt(data1.^2+data2.^2);
     84        datasize(1)=md.mesh.numberofvertices;
     85        %---> go to node data
     86end
    10987
    110         %treat the case datasize(1)=nodes2d
    111         if (md.mesh.dimension==3 & datasize(1)==numberofvertices2d),
    112                 data=project3d(md,'vector',data,'type','node');
    113                 datasize(1)=md.mesh.numberofvertices;
    114                 %---> go to node data
    115         end
     88%treat the case datasize(1)=nodes2d
     89if (md.mesh.dimension==3 & datasize(1)==numberofvertices2d),
     90        data=project3d(md,'vector',data,'type','node');
     91        datasize(1)=md.mesh.numberofvertices;
     92        %---> go to node data
     93end
    11694
    117         %treat the case datasize(1)=nodes2d
    118         if (md.mesh.dimension==3 & datasize(1)==numberofelements2d),
    119                 data=project3d(md,'vector',data,'type','element');
    120                 datasize(1)=md.mesh.numberofelements;
    121                 %---> go to node data
    122         end
     95%treat the case datasize(1)=nodes2d
     96if (md.mesh.dimension==3 & datasize(1)==numberofelements2d),
     97        data=project3d(md,'vector',data,'type','element');
     98        datasize(1)=md.mesh.numberofelements;
     99        %---> go to node data
     100end
    123101
    124         %smoothing?
    125         if exist(options,'smooth')
    126                 data=averaging(md,data,getfieldvalue(options,'smooth'));
    127                 datasize(1)=md.mesh.numberofvertices;
    128                 %---> go to node data
    129         end
     102%smoothing?
     103if exist(options,'smooth')
     104        data=averaging(md,data,getfieldvalue(options,'smooth'));
     105        datasize(1)=md.mesh.numberofvertices;
     106        %---> go to node data
    130107end
    131108

  • issm/trunk-jpl/src/m/plot/processmesh.m

    r16325 r16470  
    2121end
    2222
    23 if (isempty(data) | ~isstruct(data)),
     23if isempty(data)
    2424        %first load x,y, etc ... to speed up plot
    2525
    … …  
    7070                elements=elements2d;
    7171        end
    72 else
    73         %Process Patch
    74         if (md.mesh.dimension==2),
    75                 elements=transpose(reshape(1:3*md.mesh.numberofelements,3,md.mesh.numberofelements));
    76                 x=transpose(reshape(md.mesh.x(data.index)',1,3*md.mesh.numberofelements));
    77                 y=transpose(reshape(md.mesh.y(data.index)',1,3*md.mesh.numberofelements));
    78                 z=zeros(3*md.mesh.numberofelements,1);
    79                 is2d=1;
    80         else
    81                 elements=transpose(reshape(1:6*md.mesh.numberofelements,6,md.mesh.numberofelements));
    82                 x=transpose(reshape(md.mesh.x(data.index)',1,6*md.mesh.numberofelements));
    83                 y=transpose(reshape(md.mesh.y(data.index)',1,6*md.mesh.numberofelements));
    84                 z=transpose(reshape(md.mesh.z(data.index)',1,6*md.mesh.numberofelements));
    85                 is2d=0;
    86         end
    8772end
    8873

  • issm/trunk-jpl/src/m/regional/regionaltransient2d.m

    r15771 r16470  
    102102        for t=find(numElements==1)
    103103                if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0,
    104                         vx=PatchToVec(md1.results.TransientSolution(t).Vx);
    105                         vy=PatchToVec(md1.results.TransientSolution(t).Vy);
    106                         thickness=PatchToVec(md1.results.TransientSolution(t).Thickness);
     104                        vx=md1.results.TransientSolution(t).Vx;
     105                        vy=md1.results.TransientSolution(t).Vy;
     106                        thickness=md1.results.TransientSolution(t).Thickness;
    107107                        spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)];
    108108                        spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vy,md2.mesh.x,md2.mesh.y)];

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r15135 r16470  
    1111end
    1212
    13 %process patch if necessary
    14 results=MatlabProcessPatch(results);
    1513
    1614function results=parseresultsfromdiskioserial(filename) % {{{
    … …  
    7674
    7775        %Add result
    78         if strcmpi(result.fieldname,'Patch'),
    79                 results(result.step).(result.fieldname)=[0 result.N];
    80         else
    81                 results(result.step).(result.fieldname)=NaN;
    82         end
     76        results(result.step).(result.fieldname)=NaN;
    8377
    8478        %read next result
    … …  
    9084result=ReadDataDimensions(fid);
    9185while ~isempty(result),
    92 
    93         %Add result
    94         if strcmpi(result.fieldname,'Patch'),
    95                 patchdimensions=results(result.step).(result.fieldname);
    96                 results(result.step).(result.fieldname)=[patchdimensions(1)+result.M result.N];
    97         end
    98 
    9986        %read next result
    10087        result=ReadDataDimensions(fid);
    101 end
    102 
    103 %allocate patches
    104 for i=1:length(results),
    105         results(i).Patch=zeros(results(i).Patch(1),results(i).Patch(2));
    106         results(i).counter=1; %use to index into the patch
    10788end
    10889
    … …  
    11798
    11899        %Add result
    119         if strcmpi(result.fieldname,'Patch'),
    120                 counter=results(result.step).counter;
    121                 counter2=counter+size(result.field,1)-1;
    122                 results(result.step).(result.fieldname)(counter:counter2,:)=result.field;
    123 
    124                 %increment counter:
    125                 results(result.step).counter=counter2+1;
    126         else
    127                 results(result.step).(result.fieldname)=result.field;
    128         end
     100        results(result.step).(result.fieldname)=result.field;
    129101
    130102        %read next result

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

    r15135 r16470  
    44import results as resultsclass
    55from MatlabFuncs import *
    6 from MatlabProcessPatch import *
    76
    87def parseresultsfromdisk(filename,iosplit):
    … …  
    1817        else:
    1918                results=parseresultsfromdiskioserial(filename)
    20 
    21         #process patch if necessary
    22         results=MatlabProcessPatch(results)
    2319
    2420        return results
    … …  
    9692
    9793                #Add result
    98                 if strcmpi(result['fieldname'],'Patch'):
    99                         setattr(results[result['step']-1],result['fieldname'],[0,result['N']])
    100                 else:
    101                         setattr(results[result['step']-1],result['fieldname'],float('NaN'))
     94                setattr(results[result['step']-1],result['fieldname'],float('NaN'))
    10295
    10396                #read next result
    … …  
    109102        while result:
    110103
    111                 #Add result
    112                 if strcmpi(result['fieldname'],'Patch'):
    113                         patchdimensions=getattr(results[result['step']-1],result['fieldname'])
    114                         setattr(results[result['step']-1],result['fieldname'],[patchdimensions[0]+result['M'],result['N']])
    115 
    116104                #read next result
    117105                result=ReadDataDimensions(fid)
    118 
    119         #allocate patches
    120         for result in results.itervalues():
    121                 if 'Patch' in result:
    122                         setattr(result,'Patch',numpy.zeros((result['Patch'][0],result['Patch'][1])))
    123                         setattr(result,'counter',0)    #use to index into the patch
    124106
    125107        #third pass, this time to read the real information
    … …  
    137119
    138120                #Add result
    139                 if strcmpi(result['fieldname'],'Patch'):
    140                         counter=results[result['step']-1].counter
    141                         counter2=counter+result['field'].shape[0]-1
    142                         getattr(results[result['step']-1],result['fieldname'])[counter:counter2,:]=result['field']
    143                         #increment counter:
    144                         results[result['step']-1].counter=counter2+1
    145                 else:
    146                         setattr(results[result['step']-1],result['fieldname'],result['field'])
     121                setattr(results[result['step']-1],result['fieldname'],result['field'])
    147122
    148123                #read next result

  • issm/trunk-jpl/src/m/solve/process_solve_options.m

    r16181 r16470  
    1010solution_type=getfieldvalue(options,'solution_type');
    1111if ~ismember(solution_type,[StressbalanceSolutionEnum(),MasstransportSolutionEnum(),ThermalSolutionEnum(),...
    12                 SteadystateSolutionEnum(),TransientSolutionEnum(),EnthalpySolutionEnum(),...
     12                SteadystateSolutionEnum(),TransientSolutionEnum(),...
    1313                BalancethicknessSolutionEnum(),BalancethicknessSoftSolutionEnum(),...
    1414                BalancevelocitySolutionEnum(),BedSlopeSolutionEnum(),...

  • issm/trunk-jpl/src/m/solve/process_solve_options.py

    r16008 r16470  
    1919        solution_type=options.getfieldvalue('solution_type')
    2020        if solution_type not in (StressbalanceSolutionEnum(),MasstransportSolutionEnum(),ThermalSolutionEnum(),\
    21                         SteadystateSolutionEnum(),TransientSolutionEnum(),EnthalpySolutionEnum(),\
     21                        SteadystateSolutionEnum(),TransientSolutionEnum(),\
    2222                        BalancethicknessSolutionEnum(),BalancevelocitySolutionEnum(),\
    2323                        BedSlopeSolutionEnum(),SurfaceSlopeSolutionEnum(),\
Note: See TracChangeset for help on using the changeset viewer.